Mesh_3: Specific domain for Poisson reconstruction (#8771)

## Summary of Changes

Creating Poisson_mesh_domain including optimizations for Poisson from
Surface_mesher

## Release Management

* Affected package(s): Mesh_3, Poisson_reconstruction_3
* small feature:
https://cgalwiki.geometryfactory.com/CGAL/Members/wiki/Features/Small_Features/Poisson_mesh_domain
-- pre-approved 18 April 2025
* Issue(s) solved (if any): fix #8266
This commit is contained in:
Sebastien Loriot 2025-05-07 09:42:25 +02:00 committed by GitHub
commit ac763ef561
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
15 changed files with 835 additions and 198 deletions

View File

@ -58,6 +58,9 @@
- `initial_points_generator` : enables the user to specify a functor that generates initial points,
- `initial_points` : enables the user to specify a `Range` of initial points.
- Added a new meshing parameter `surface_only`, to improve performances when the user is only interested in surface mesh generation.
### [Poisson Surface Reconstruction](https://doc.cgal.org/6.1/Manual/packages.html#PkgPoissonSurfaceReconstruction3)
- Added a new mesh domain `Poisson_mesh_domain_3` that integrates some optimizations from the deprecated 3D Surface Mesh Generation package.
### [3D Subdivision Methods](https://doc.cgal.org/6.1/Manual/packages.html#PkgSurfaceSubdivisionMethod3)

View File

@ -7,13 +7,8 @@
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Stéphane Tayeb, Aymeric PELLE
//
//******************************************************************************
// File Description :
// class Labeled_mesh_domain_3. See class description.
//******************************************************************************
#ifndef CGAL_LABELED_MESH_DOMAIN_3_H
#define CGAL_LABELED_MESH_DOMAIN_3_H
@ -875,17 +870,15 @@ public:
* function. The domain to be discretized is assumed to be the domain where
* the function has negative values.
*
* The method takes as argument a bounding sphere which is required to
* circumscribe the surface and to have its center inside the domain.
*
* \tparam Function a type compatible with the signature `FT(Point_3)`: it takes a point as argument,
* and returns a scalar value. That object must be model of `CopyConstructible`
* \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
* \tparam Bounding_object either a bounding sphere (of type `Sphere_3`), a bounding box (type `Bbox_3`),
* or a bounding `Iso_cuboid_3`
* or a bounding `Iso_cuboid_3` which is required to circumscribe
* the surface and to have its center inside the domain.
*
* \param function the implicit function
* \param bounding_object object boundint the meshable domain and its center is inside the domain.
* \param bounding_object object bounding the meshable domain and its center is inside the domain.
* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below:
*
* \cgalNamedParamsBegin
@ -936,12 +929,11 @@ public:
/// @}
#ifndef DOXYGEN_RUNNING
template<typename CGAL_NP_TEMPLATE_PARAMETERS>
static Labeled_mesh_domain_3 create_implicit_mesh_domain(const CGAL_NP_CLASS& np = parameters::default_values())
static Labeled_mesh_domain_3 create_implicit_mesh_domain(const CGAL_NP_CLASS &np)
{
static_assert(!parameters::is_default_parameter<CGAL_NP_CLASS, internal_np::function_param_t>::value, "Value for required parameter not found");
static_assert(!parameters::is_default_parameter<CGAL_NP_CLASS, internal_np::bounding_object_param_t>::value, "Value for required parameter not found");
using parameters::get_parameter;
return create_implicit_mesh_domain(parameters::get_parameter(np, internal_np::function_param),
parameters::get_parameter(np, internal_np::bounding_object_param),
np);
@ -1139,11 +1131,7 @@ public:
/*
* Returns a point in the intersection of the primitive `type`
* with some boundary surface.
* `Type1` is either `Segment_3`, `Ray_3` or `Line_3`.
* The integer `dimension` is set to the dimension of the lowest
* dimensional face in the input complex containing the returned point, and
* `index` is set to the index to be stored at a mesh vertex lying
* on this face.
* `Type` is either `Segment_3`, `Ray_3` or `Line_3`.
*/
struct Construct_intersection
{
@ -1170,10 +1158,10 @@ public:
private:
/*
* Returns a point in the intersection of [a,b] with the surface
* `a` must be the source point, and `b` the out point. It's important
* Returns a point in the intersection of `[a,b]` with the surface
* `a` must be the source point, and `b` the out point. It is important
* because it drives bisection cuts.
* Indeed, the returned point is the first intersection from `[a,b]`
* Indeed, the returned point is the first intersection of `[a,b]`
* with a subdomain surface.
*/
Intersection operator()(const Point_3& a, const Point_3& b) const
@ -1196,7 +1184,7 @@ public:
// If both extremities are in the same subdomain,
// there is no intersection.
// This should not happen...
// Should only be able to happen during initial point generation.
if( value_at_p1 == value_at_p2 )
{
return Intersection();

View File

@ -13,8 +13,7 @@ This \cgal component implements a surface reconstruction method which
takes as input point sets with oriented normals and computes an
implicit function. We assume that the input points contain no outliers
and little noise. The output surface mesh is generated by extracting
an isosurface of this function with the \cgal Surface Mesh Generator
\cgalCite{cgal:ry-gsddrm-06} or potentially with any other surface
an isosurface of this function with the \ref PkgMesh3 or potentially with any other surface
contouring algorithm.
\cgalFigureBegin{Poisson_surface_reconstruction_3figintroduction,introduction.jpg}
@ -75,8 +74,8 @@ during refinement is set to zero. It then solves for a scalar
indicator function \f$ f\f$ represented as a piecewise linear function
over the refined triangulation. More specifically, it solves for the
Poisson equation \f$ \Delta f = div(\mathbf{n})\f$ at each vertex of
the triangulation using a sparse linear solver. Eventually, the \cgal
surface mesh generator extracts an isosurface with function value set
the triangulation using a sparse linear solver. Eventually, the \ref PkgMesh3
extracts an isosurface with function value set
by default to be the median value of \f$ f\f$ at all input points.
\section Poisson_surface_reconstruction_3Function Reconstruction Function
@ -86,9 +85,7 @@ provided. It takes points with normals as input and handles the whole
reconstruction pipeline :
- it computes the implicit function
- it reconstructs the surface with a given precision using the \cgal
surface mesh generator based on Delaunay refinement
\cgalCite{cgal:ry-gsddrm-06} \cgalCite{cgal:bo-pgsms-05}
- it reconstructs the surface with a given precision using the \ref PkgMesh3
- it outputs the result in a polygon mesh.
This function aims at providing a quick and user-friendly API for
@ -118,26 +115,19 @@ The following example reads a point set, creates a Poisson implicit function and
\subsection Poisson_surface_reconstruction_3Contouring Contouring
The computed implicit functions can be iso-contoured to reconstruct a
surface by using the \cgal surface mesh generator
\cgalCite{cgal:ry-gsddrm-06} \cgalCite{cgal:bo-pgsms-05} :
surface by using the \ref PkgMesh3 component, and in particular the function
`make_mesh_3()` with the `surface_only()` parameter to only mesh the surface.
`make_surface_mesh()`
The parameter `Tag` affects the behavior of `make_surface_mesh()`:
The following `Tag` parameters affect the behavior of `make_mesh_3()`:
- `Manifold_tag`: the output mesh is guaranteed to be a manifold surface without boundary.
- `Manifold_with_boundary_tag`: the output mesh is guaranteed to be manifold and may have boundaries.
- `Non_manifold_tag`: the output mesh has no guarantee and hence is outputted as a polygon soup.
\subsection Poisson_surface_reconstruction_3Output Output
The surface reconstructed by `make_surface_mesh()` is required to be a
model of the concept `SurfaceMeshComplex_2InTriangulation_3`, a data
structure devised to represent a two dimensional complex embedded into
a three dimensional triangulation.
`SurfaceMeshComplex_2InTriangulation_3` defines the methods to traverse the reconstructed surface, and e.g. convert it to a triangle soup.
The surface reconstructed by `make_mesh_3()` is required to be a model of the concept
`MeshComplex_3InTriangulation_3`, a data structure devised to represent a three dimensional complex embedded into a three dimensional triangulation. The surface facets can then be extracted into a face graph by `facets_in_complex_3_to_triangle_mesh()`.
Other \cgal components provide functions to write the reconstructed
surface mesh to the %Object File Format (OFF) \cgalCite{cgal:p-gmgv16-96}
@ -218,7 +208,7 @@ function over the tetrahedra of a 3D Delaunay triangulation
constructed from the input points then refined through Delaunay
refinement. For this reason, any iso-surface is also piecewise linear
and hence may contain sharp creases. As the contouring algorithm
`make_surface_mesh()` expects a smooth implicit function these
`make_mesh_3()` expects a smooth implicit function these
sharp creases may create spurious clusters of vertices in the final
reconstructed surface mesh when setting a small mesh sizing or surface
approximation error parameter (see
@ -339,136 +329,238 @@ Left: 5K points sampled on a mechanical piece with sharp features (creases, dart
We provide some performance numbers for scanning data. We measure the Poisson implicit function computation time,
the contouring time for a range of approximation distances, the memory occupancy as well as the influence of
the point set simplification. The machine used is a PC running Windows 7 64 bits with an Intel CPU Core 2 Duo
processor clocked at 2.81 GHz and with 8 GB of RAM. The software is compiled with Visual \CC 2010 (VC9) compiler
the point set simplification. The machine used is a PC running Windows 10 64 bits with an Intel CPU Core i7-11850H
processor with 8 cores and 32 GB of RAM. The software is compiled with Visual \CC 2022 compiler
with the 03 option which maximizes speed. All measurements were done using the \ref thirdpartyEigen "Eigen" library.
\subsection SurfReconstPerfPIF Poisson Implicit Function
The point set chosen for benchmarking the Poisson implicit function is the Bimba con Nastrino point set
(1.6 million points) depicted by \cgalFigureRef{Poisson_surface_reconstruction_3-fig-contouring_bench}.
The point set chosen for benchmarking the Poisson implicit function is the Lucy statue point set from the
<a href="https://graphics.stanford.edu/data/3Dscanrep/">The Stanford 3D Scanning Repository</a>
(originally 14 million points, here downsampled to 2.9 million points) depicted by \cgalFigureRef{Poisson_surface_reconstruction_3-fig-contouring_bench}.
We measure the Poisson implicit function computation (i.e., the call to
`Poisson_reconstruction_function::compute_implicit_function()` denoted by Poisson solve hereafter)
for this point set as well as for simplified versions obtained through random simplification.
The following table provides Poisson solve computation times in seconds for an increasing number of points.
<CENTER>
<TABLE CELLSPACING=5 >
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
Number of points (x1000)
<TD class="math" ALIGN=CENTER NOWRAP>
Poisson solve duration (in s)
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
30
<TD class="math" ALIGN=CENTER NOWRAP>
3.3
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
60
<TD class="math" ALIGN=CENTER NOWRAP>
15
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
100
<TD class="math" ALIGN=CENTER NOWRAP>
25
7.7
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
250
120
<TD class="math" ALIGN=CENTER NOWRAP>
96
18.1
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
500
237.5
<TD class="math" ALIGN=CENTER NOWRAP>
150
35.1
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,000
375
<TD class="math" ALIGN=CENTER NOWRAP>
249
64
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,800
750
<TD class="math" ALIGN=CENTER NOWRAP>
478
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
129
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,500
<TD class="math" ALIGN=CENTER NOWRAP>
303
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
2,900
<TD class="math" ALIGN=CENTER NOWRAP>
486
<TR><TD ALIGN=CENTER NOWRAP COLSPAN=2><HR>
</TABLE>
</CENTER>
\subsection SurfReconstPerfCont Contouring
The point set chosen for benchmarking the contouring stage is the Bimba con Nastrino point
set simplified to 100k points. We measure the contouring (i.e.\ the call to `make_surface_mesh()`)
The point set chosen for benchmarking the contouring stage is the Lucy point
set simplified to 2.9M points. We measure the contouring (i.e., the calls to `make_mesh_3()` and `facets_in_complex_3_to_triangle_mesh()`)
duration and the reconstruction error for a range of approximation distances.
The reconstruction error is expressed as the average distance from input points to the reconstructed surface
in mm (the Bimba con Nastrino statue is 324 mm tall).
The reconstruction error is expressed as the average distance from input points to the reconstructed surface in mm (the Lucy statue is 1597 mm tall).
<CENTER>
<TABLE CELLSPACING=5 >
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=3><HR>
<TR><TD ALIGN=CENTER NOWRAP COLSPAN=4><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
Approx. distance (*average spacing)
<TD class="math" ALIGN=CENTER NOWRAP>
Contouring duration (in s)
Contouring duration single-thread (in s)
<TD class="math" ALIGN=CENTER NOWRAP>
Contouring duration parallel (in s)
<TD class="math" ALIGN=CENTER NOWRAP>
Reconstruction error (mm)
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=3><HR>
<TR><TD ALIGN=CENTER NOWRAP COLSPAN=4><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
0.05
<TD class="math" ALIGN=CENTER NOWRAP>
582
<TD class="math" ALIGN=CENTER NOWRAP>
112
<TD class="math" ALIGN=CENTER NOWRAP>
0.114
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
0.1
<TD class="math" ALIGN=CENTER NOWRAP>
19.2
221
<TD class="math" ALIGN=CENTER NOWRAP>
0.055
26
<TD class="math" ALIGN=CENTER NOWRAP>
0.119
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
0.15
<TD class="math" ALIGN=CENTER NOWRAP>
104
<TD class="math" ALIGN=CENTER NOWRAP>
25
<TD class="math" ALIGN=CENTER NOWRAP>
0.129
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
0.2
<TD class="math" ALIGN=CENTER NOWRAP>
69.4
<TD class="math" ALIGN=CENTER NOWRAP>
9.2
<TD class="math" ALIGN=CENTER NOWRAP>
0.14
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
0.25
<TD class="math" ALIGN=CENTER NOWRAP>
6.9
53.6
<TD class="math" ALIGN=CENTER NOWRAP>
0.106
7.0
<TD class="math" ALIGN=CENTER NOWRAP>
0.151
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
0.5
<TD class="math" ALIGN=CENTER NOWRAP>
3.2
25.2
<TD class="math" ALIGN=CENTER NOWRAP>
0.18
3.6
<TD class="math" ALIGN=CENTER NOWRAP>
0.209
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
0.75
<TD class="math" ALIGN=CENTER NOWRAP>
16.4
<TD class="math" ALIGN=CENTER NOWRAP>
4.9
<TD class="math" ALIGN=CENTER NOWRAP>
0.209
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1
<TD class="math" ALIGN=CENTER NOWRAP>
1.65
12.4
<TD class="math" ALIGN=CENTER NOWRAP>
0.36
2.3
<TD class="math" ALIGN=CENTER NOWRAP>
0.33
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1.5
<TD class="math" ALIGN=CENTER NOWRAP>
8.2
<TD class="math" ALIGN=CENTER NOWRAP>
1.4
<TD class="math" ALIGN=CENTER NOWRAP>
0.455
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
2
<TD class="math" ALIGN=CENTER NOWRAP>
6.1
<TD class="math" ALIGN=CENTER NOWRAP>
1.1
<TD class="math" ALIGN=CENTER NOWRAP>
0.59
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
3
<TD class="math" ALIGN=CENTER NOWRAP>
4.0
<TD class="math" ALIGN=CENTER NOWRAP>
0.8
<TD class="math" ALIGN=CENTER NOWRAP>
0.76
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=3><HR>
0.87
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
5
<TD class="math" ALIGN=CENTER NOWRAP>
2.3
<TD class="math" ALIGN=CENTER NOWRAP>
0.7
<TD class="math" ALIGN=CENTER NOWRAP>
1.50
<TR>
<TD ALIGN=CENTER NOWRAP COLSPAN=4><HR>
</TABLE>
</CENTER>
\cgalFigureBegin{Poisson_surface_reconstruction_3-fig-contouring_bench,contouring_bench.jpg}
\cgalFigureAnchor{Poisson_surface_reconstruction_3-fig-contouring_bench}
<center>
<img src="contouring_bench.jpg" style="max-width:80%;"/>
</center>
\cgalFigureCaptionBegin{Poisson_surface_reconstruction_3-fig-contouring_bench}
Contouring duration (in s) and reconstruction error (mm)
against several approximation distance parameters
for the Bimba con Nastrino point set simplified to 100k points.
\cgalFigureEnd
for the Lucy point set simplified to 100k points.
\cgalFigureCaptionEnd
\subsection SurfReconstPerfMem Memory
We measure the memory occupancy for the reconstruction of the full Bimba con Nastrino point
set (1.8 millions points) as well as for simplified versions.\n
We measure the memory occupancy for the reconstruction of the Lucy point
set (2.9 millions points) as well as for further simplified versions.\n
The Poisson implicit function computation has a memory peak when solving the Poisson linear
system using the sparse linear solver.
<CENTER>
<TABLE CELLSPACING=5 >
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
@ -479,47 +571,59 @@ Memory occupancy (MBytes)
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
30
<TD class="math" ALIGN=CENTER NOWRAP>
128
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
60
<TD class="math" ALIGN=CENTER NOWRAP>
180
226
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
100
120
<TD class="math" ALIGN=CENTER NOWRAP>
270
431
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
250
237.5
<TD class="math" ALIGN=CENTER NOWRAP>
790
813
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
500
375
<TD class="math" ALIGN=CENTER NOWRAP>
1300
1,232
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,000
750
<TD class="math" ALIGN=CENTER NOWRAP>
2200
2,283
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,800
1,500
<TD class="math" ALIGN=CENTER NOWRAP>
3800
4,042
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
2,900
<TD class="math" ALIGN=CENTER NOWRAP>
6,868
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
</TABLE>
</CENTER>
\subsection SurfReconstPerfPSS Point Set Simplification
Due to the memory limitations described above, we recommend to simplify the point sets captured by laser scanners.\n
We measure the reconstruction error for the Bimba con Nastrino point set (1.6M points) as well as for
We measure the reconstruction error for the Lucy point set (2.9M points) as well as for
simplified versions. All reconstructions use the recommended contouring parameter
`approximation distance = 0.25 * the input point` set's average spacing.
The reconstruction error is expressed as the average distance from input points to the reconstructed surface in mm
(the Bimba con Nastrino statue is 324 mm tall).
(the Lucy statue is 1597 mm tall).
<CENTER>
<TABLE CELLSPACING=5 >
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
@ -530,46 +634,67 @@ Reconstruction error (mm)
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
3.75
<TD class="math" ALIGN=CENTER NOWRAP>
9.88395
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
7.5
<TD class="math" ALIGN=CENTER NOWRAP>
5.81843
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
15
<TD class="math" ALIGN=CENTER NOWRAP>
3.13479
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
30
<TD class="math" ALIGN=CENTER NOWRAP>
2.25391
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
60
<TD class="math" ALIGN=CENTER NOWRAP>
0.27
1.42965
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
120
<TD class="math" ALIGN=CENTER NOWRAP>
0.15
1.17589
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
250
237.5
<TD class="math" ALIGN=CENTER NOWRAP>
0.11
0.99509
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
500
375
<TD class="math" ALIGN=CENTER NOWRAP>
0.079
0.75215
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,000
750
<TD class="math" ALIGN=CENTER NOWRAP>
0.066
0.344654
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,500
<TD class="math" ALIGN=CENTER NOWRAP>
0.061
0.225341
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,600
2,900
<TD class="math" ALIGN=CENTER NOWRAP>
0.06
0.150947
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
</TABLE>
</CENTER>
\cgalFigureBegin{Poisson_surface_reconstruction_3-fig-simplification_bench,simplification_bench.jpg}
Reconstruction error (mm) against number of points
for the Bimba con Nastrino point set with 1.6M points
for the Lucy point set with 2.9M points
as well as for simplified versions.
\cgalFigureEnd

View File

@ -1,8 +1,11 @@
Manual
Kernel_23
BGL
STL_Extension
Mesh_3
Algebraic_foundations
Circulator
Stream_support
Point_set_processing_3
Solver_interface
SMDS_3

Binary file not shown.

Before

Width:  |  Height:  |  Size: 66 KiB

After

Width:  |  Height:  |  Size: 133 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 65 KiB

After

Width:  |  Height:  |  Size: 53 KiB

View File

@ -7,7 +7,7 @@
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Poisson_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/facets_in_complex_3_to_triangle_mesh.h>
@ -44,7 +44,7 @@ namespace params = CGAL::parameters;
template<typename Concurrency_tag, typename PointSet>
void poisson_reconstruction(const PointSet& points, const char* output)
{
typedef CGAL::Labeled_mesh_domain_3<Kernel> Mesh_domain;
typedef CGAL::Poisson_mesh_domain_3<Kernel> Mesh_domain;
typedef typename CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
@ -111,15 +111,14 @@ void poisson_reconstruction(const PointSet& points, const char* output)
params::facet_size = sm_radius * average_spacing,
params::facet_distance = sm_distance * average_spacing);
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(function, bsphere,
Mesh_domain domain = Mesh_domain::create_Poisson_mesh_domain(function, bsphere,
params::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
// Generates surface mesh with manifold option
std::cout << "Start meshing...";
std::cout.flush();
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::no_exude()
.no_perturb()
params::surface_only()
.manifold_with_boundary());
time.stop();

View File

@ -18,7 +18,7 @@
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Poisson_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/facets_in_complex_3_to_triangle_mesh.h>
@ -56,7 +56,7 @@ typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef CGAL::Poisson_reconstruction_function<Kernel> Poisson_reconstruction_function;
// Mesh_3
typedef CGAL::Labeled_mesh_domain_3<Kernel> Mesh_domain;
typedef CGAL::Poisson_mesh_domain_3<Kernel> Mesh_domain;
typedef typename CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
@ -327,12 +327,12 @@ int main(int argc, char * argv[])
<< " manifold_with_boundary()\n";
// Defines mesh domain
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(function, bsphere,
Mesh_domain domain = Mesh_domain::create_Poisson_mesh_domain(function, bsphere,
CGAL::parameters::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
// Generates mesh with manifold option
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
CGAL::parameters::no_exude().no_perturb()
CGAL::parameters::surface_only()
.manifold_with_boundary());
// Prints status

View File

@ -5,7 +5,7 @@
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Poisson_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/facets_in_complex_3_to_triangle_mesh.h>
@ -32,7 +32,7 @@ typedef Kernel::Sphere_3 Sphere;
typedef std::vector<Point_with_normal> PointList;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef CGAL::Poisson_reconstruction_function<Kernel> Poisson_reconstruction_function;
typedef CGAL::Labeled_mesh_domain_3<Kernel> Mesh_domain;
typedef CGAL::Poisson_mesh_domain_3<Kernel> Mesh_domain;
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
@ -85,12 +85,12 @@ int main(void)
CGAL::parameters::facet_distance = sm_distance*average_spacing);
// Defines mesh domain
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(function, bsphere,
Mesh_domain domain = Mesh_domain::create_Poisson_mesh_domain(function, bsphere,
CGAL::parameters::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
// Generates mesh with manifold option
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
CGAL::parameters::no_exude().no_perturb()
CGAL::parameters::surface_only()
.manifold_with_boundary());
const Tr& tr = c3t3.triangulation();

View File

@ -0,0 +1,409 @@
// Copyright (c) 2025 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Sven Oesau
#ifndef CGAL_POISSON_MESH_DOMAIN_3_H
#define CGAL_POISSON_MESH_DOMAIN_3_H
#include <CGAL/license/Poisson_surface_reconstruction_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Poisson_reconstruction_function.h>
namespace CGAL {
/*!
\ingroup PkgPoissonSurfaceReconstruction3Ref
\brief The class `Poisson_mesh_domain_3` derives from `Labeled_mesh_domain_3` for the handling of `Poisson_reconstruction_function`.
This class has a constructor taking a labeling function. It has also a static template member
function that acts as named constructor:
<ul><li>`create_Poisson_mesh_domain()`</li>, to create a domain from a `Poisson_reconstruction_function`</ul>
\tparam BGT is a geometric traits class that provides
the basic operations to implement intersection tests and intersection computations through a bisection
method. This parameter must be instantiated with a model of the concept `BisectionGeometricTraits_3`.
\cgalModels{MeshDomain_3}
\sa `CGAL::Labeled_mesh_domain_3`
\sa `CGAL::make_mesh_3()`
*/
template<class BGT>
class Poisson_mesh_domain_3
#ifndef DOXYGEN_RUNNING
: public Labeled_mesh_domain_3<BGT>
#endif
{
public:
using Base = Labeled_mesh_domain_3<BGT>;
typedef typename Base::Subdomain Subdomain;
typedef typename Base::Subdomain_index Subdomain_index;
typedef typename Base::Surface_patch_index Surface_patch_index;
typedef typename Base::Intersection Intersection;
// Type of indexes for cells of the input complex
typedef std::optional<Surface_patch_index> Surface_patch;
// Type of indexes to characterize the lowest dimensional face of the input
// complex on which a vertex lie
typedef typename CGAL::Mesh_3::internal::Index_generator<Subdomain_index, Surface_patch_index>::Index Index;
// Geometric object types
#ifdef DOXYGEN_RUNNING
/// \name Types imported from the geometric traits class
///@{
/// The point type of the geometric traits class
typedef typename Geom_traits::Point_3 Point_3;
/// The sphere type of the geometric traits class
typedef typename Geom_traits::Sphere_3 Sphere_3;
/// The iso-cuboid type of the geometric traits class
typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid_3;
/// The bounding box type
typedef CGAL::Bbox_3 Bbox_3;
/// The number type (a field type) of the geometric traits class
typedef typename Geom_traits::FT FT;
/// The ray type of the geometric traits class
typedef typename Geom_traits::Ray_3 Ray_3;
/// The line type of the geometric traits class
typedef typename Geom_traits::Line_3 Line_3;
/// The segment type of the geometric traits class
typedef typename Geom_traits::Segment_3 Segment_3;
/// The Poisson function type
typedef CGAL::Poisson_reconstruction_function<Geom_traits> Function;
///@}
#else
/// The point type of the geometric traits class
typedef typename BGT::Point_3 Point_3;
/// The sphere type of the geometric traits class
typedef typename BGT::Sphere_3 Sphere_3;
/// The iso-cuboid type of the geometric traits class
typedef typename BGT::Iso_cuboid_3 Iso_cuboid_3;
/// The bounding box type
typedef CGAL::Bbox_3 Bbox_3;
/// The number type (a field type) of the geometric traits class
typedef typename BGT::FT FT;
/// The ray type of the geometric traits class
typedef typename BGT::Ray_3 Ray_3;
/// The line type of the geometric traits class
typedef typename BGT::Line_3 Line_3;
/// The segment type of the geometric traits class
typedef typename BGT::Segment_3 Segment_3;
/// The Poisson function type
typedef CGAL::Poisson_reconstruction_function<BGT> Function;
#endif
Function poisson_function;
/// \name Creation
/// @{
/*! \brief Construction from a function, a bounding object and a relative error bound.
*
* \tparam Bounding_object either a bounding sphere (of type `Sphere_3`), a bounding box (type `Bbox_3`),
* or a bounding `Iso_cuboid_3`
* \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
*
* \param function the Poisson reconstruction function
* \param bounding_object the bounding object bounding the meshable space.
* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below:
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{relative_error_bound}
* \cgalParamDescription{the relative error bound used to compute intersection points between the implicit surface and query segments.
* The bisection is stopped when the length of the intersected segment is less than the product
* of `relative_error_bound` by the diameter of the bounding object.}
* \cgalParamDefault{FT(1e-3)}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*/
template<typename Bounding_object, typename CGAL_NP_TEMPLATE_PARAMETERS>
Poisson_mesh_domain_3(const Function& function,
const Bounding_object& bounding_object,
const CGAL_NP_CLASS& np = parameters::default_values()
#ifndef DOXYGEN_RUNNING
, typename std::enable_if<!is_named_function_parameter<Function>>::type* = nullptr
#endif // DOXYGEN_RUNNING
)
: Base(make_implicit_to_labeling_function_wrapper<BGT>(function), bounding_object, np),
poisson_function(function)
{}
/*! \brief Construction from a function, a bounding object and a relative error bound.
*
* \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
*
* \param function the Poisson reconstruction function
* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below:
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{relative_error_bound}
* \cgalParamDescription{the relative error bound used to compute intersection points between the implicit surface and query segments.
* The bisection is stopped when the length of the intersected segment is less than the product
* of `relative_error_bound` by the diameter of the bounding object.}
* \cgalParamDefault{FT(1e-3)}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*/
template<typename CGAL_NP_TEMPLATE_PARAMETERS>
Poisson_mesh_domain_3(const Function & function,
const CGAL_NP_CLASS& np = parameters::default_values()
#ifndef DOXYGEN_RUNNING
, typename std::enable_if<!is_named_function_parameter<Function>>::type * = nullptr
#endif // DOXYGEN_RUNNING
)
: Base(make_implicit_to_labeling_function_wrapper<BGT>(function), function.bounding_sphere(), np),
poisson_function(function)
{}
///@}
#ifndef DOXYGEN_RUNNING
template<typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT>
Poisson_mesh_domain_3(const CGAL_NP_CLASS& np)
: Base(np)
{}
// Overload handling parameters passed with operator=
template<typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_1,
typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_2,
typename ... NP>
Poisson_mesh_domain_3(const Function& function,
const CGAL_NP_CLASS_1& np1,
const CGAL_NP_CLASS_2& np2,
const NP& ... nps)
: Base(internal_np::combine_named_parameters(
CGAL::parameters::function(make_implicit_to_labeling_function_wrapper<BGT>(function)), np1, np2, nps...)),
poisson_function(function)
{}
#endif
/// \name Creation of domains from Poisson implicit functions
/// @{
/*! \brief Construction from a Poisson implicit function
*
* This static method is a <em>named constructor</em>. It constructs a domain
* whose bounding surface is described implicitly as the zero level set of a
* function. The domain to be discretized is assumed to be the domain where
* the function has negative values.
*
* \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
* \tparam Bounding_object either a bounding sphere (of type `Sphere_3`), a bounding box (type `Bbox_3`),
* or a bounding `Iso_cuboid_3` which is required to circumscribe
* the surface and to have its center inside the domain.
*
* \param function the Poisson reconstruction function
* \param bounding_object object bounding the meshable domain and its center is inside the domain.
* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below:
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{relative_error_bound}
* \cgalParamDescription{ is the relative error
* bound, relative to the diameter of the box of the image.}
* \cgalParamDefault{FT(1e-3)}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
*/
template<typename Bounding_object, typename CGAL_NP_TEMPLATE_PARAMETERS>
static Poisson_mesh_domain_3 create_Poisson_mesh_domain(const Function& function,
const Bounding_object& bounding_object,
const CGAL_NP_CLASS& np = parameters::default_values()
#ifndef DOXYGEN_RUNNING
,typename std::enable_if<!is_named_function_parameter<Function>>::type* = nullptr
#endif
)
{
using parameters::get_parameter;
using parameters::choose_parameter;
FT relative_error_bound_ = choose_parameter(get_parameter(np, internal_np::error_bound), FT(1e-3));
CGAL::Random* p_rng_ = choose_parameter(get_parameter(np, internal_np::rng), nullptr);
auto null_subdomain_index_ = choose_parameter(get_parameter(np, internal_np::null_subdomain_index_param), Null_functor());
auto construct_surface_patch_index_ = choose_parameter(get_parameter(np, internal_np::surface_patch_index), Null_functor());
return Poisson_mesh_domain_3(function,
bounding_object,
CGAL::parameters::relative_error_bound(relative_error_bound_)
.function(make_implicit_to_labeling_function_wrapper<BGT>(function))
.p_rng(p_rng_)
.null_subdomain_index(Base::create_null_subdomain_index(null_subdomain_index_))
.construct_surface_patch_index(Base::create_construct_surface_patch_index(construct_surface_patch_index_)));
}
/// @}
#ifndef DOXYGEN_RUNNING
template<typename CGAL_NP_TEMPLATE_PARAMETERS>
static Poisson_mesh_domain_3 create_Poisson_mesh_domain(const CGAL_NP_CLASS& np) {
using parameters::get_parameter;
static_assert(!parameters::is_default_parameter<CGAL_NP_CLASS, internal_np::function_param_t>::value, "Value for required parameter not found");
static_assert(!parameters::is_default_parameter<CGAL_NP_CLASS, internal_np::bounding_object_param_t>::value, "Value for required parameter not found");
return create_Poisson_mesh_domain(get_parameter(np, internal_np::function_param),
get_parameter(np, internal_np::bounding_object_param),
np);
}
// Overload handling parameters passed with operator=
template<typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_1,
typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_2,
typename ... NP>
static Poisson_mesh_domain_3 create_Poisson_mesh_domain(const CGAL_NP_CLASS_1& np1,
const CGAL_NP_CLASS_2& np2,
const NP& ... nps)
{
return create_Poisson_mesh_domain(internal_np::combine_named_parameters(np1, np2, nps...));
}
#endif
/*
* Returns a point in the intersection of the primitive `type`
* with some boundary surface.
* `Type` is either `Segment_3`, `Ray_3` or `Line_3`.
*/
struct Construct_intersection
{
Construct_intersection(const Poisson_mesh_domain_3& domain) : domain_(domain) {}
Intersection operator()(const Segment_3& s) const
{
#ifndef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
CGAL_precondition(r_domain_.do_intersect_surface_object()(s) != std::nullopt);
#endif // NOT CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
return this->operator()(s.source(), s.target());
}
Intersection operator()(const Ray_3& r) const
{
return clip_to_segment(r);
}
Intersection operator()(const Line_3& l) const
{
return clip_to_segment(l);
}
private:
/*
* Returns a point in the intersection of `[a,b]` with the surface
* `a` must be the source point, and `b` the out point. It is important
* because it drives bisection cuts.
* Indeed, the returned point is the first intersection of `[a,b]`
* with a subdomain surface.
*
* The difference from the Labeled_mesh_domain_3::Construct_intersection is that
* the underlying Delaunay triangulation in the Poisson function is used for the bisection.
*/
Intersection operator()(const Point_3& a, const Point_3& b) const
{
// Functors
typename BGT::Compute_squared_distance_3 squared_distance =
BGT().compute_squared_distance_3_object();
typename BGT::Construct_midpoint_3 midpoint =
BGT().construct_midpoint_3_object();
// Non const points
Point_3 p1 = a;
Point_3 p2 = b;
Point_3 mid = midpoint(p1, p2);
FT value_at_p1, value_at_p2;
typename Function::Cell_handle c1, c2;
bool c1_is_inf, c2_is_inf;
std::tie(value_at_p1, c1, c1_is_inf) = domain_.poisson_function.special_func(p1);
std::tie(value_at_p2, c2, c2_is_inf) = domain_.poisson_function.special_func(p2);
Subdomain_index label_at_p1 = (value_at_p1 < 0) ? 1 : 0;
Subdomain_index label_at_p2 = (value_at_p2 < 0) ? 1 : 0;
// If both extremities are in the same subdomain,
// there is no intersection.
// Should only be able to happen during initial point generation.
if(label_at_p1 == label_at_p2)
return Intersection();
// Else lets find a point (by bisection)
// Bisection ends when the point is nearer from surface than the error bound
while(true) {
if(c1 == c2) {
if(c1_is_inf) {
std::cout << "Intersection(): c1 == c2 and inf!" << std::endl;
return Intersection();
} else {
const Surface_patch_index sp_index = domain_.make_surface_index(label_at_p1, label_at_p2);
const Index index = domain_.index_from_surface_patch_index(sp_index);
return Intersection(Point_3(ORIGIN + ((value_at_p2 * (p1 - ORIGIN)) - (value_at_p1 * (p2 - ORIGIN))) /
(value_at_p2 - value_at_p1)), index, 2);
}
}
mid = midpoint(p1, p2);
// If the two points are enough close, then we return midpoint
if ( squared_distance(p1, p2) < domain_.squared_error_bound_ )
{
CGAL_assertion(value_at_p1 * value_at_p2 <= 0);
const Surface_patch_index sp_index = domain_.make_surface_index(label_at_p1, label_at_p2);
const Index index = domain_.index_from_surface_patch_index(sp_index);
return Intersection(mid, index, 2);
}
// Cannot be const: those values are modified below.
FT value_at_mid;
typename Function::Cell_handle c_at_mid;
bool c_is_inf;
std::tie(value_at_mid, c_at_mid, c_is_inf) = domain_.poisson_function.special_func(mid);
Subdomain_index label_at_mid = (value_at_mid < 0) ? 1 : 0;
// Else we must go on
// Here we consider that p1(a) is the source point. Thus, we keep p1 and
// change p2 if f(p1)!=f(p2).
// That allows us to find the first intersection from a of [a,b] with
// a surface.
if(label_at_p1 != label_at_mid && !(domain_.null(label_at_p1) && domain_.null(label_at_mid)))
{
p2 = mid;
value_at_p2 = value_at_mid;
label_at_p2 = label_at_mid;
}
else
{
p1 = mid;
value_at_p1 = value_at_mid;
label_at_p1 = label_at_mid;
}
}
}
// Clips `query` to a segment `s`, and call `operator()(s)`
template<typename Query>
Intersection clip_to_segment(const Query& query) const
{
const auto clipped = CGAL::intersection(query, domain_.bbox_);
if (clipped)
if (const Segment_3* s = std::get_if<Segment_3>(&*clipped))
return this->operator()(*s);
return Intersection();
}
const Poisson_mesh_domain_3& domain_;
};
// Returns Construct_intersection object
Construct_intersection construct_intersection_object() const
{
return Construct_intersection(*this);
}
}; // end class Poisson_mesh_domain_3
} // end namespace CGAL
#endif // CGAL_LABELED_MESH_DOMAIN_3_H

View File

@ -460,7 +460,7 @@ public:
// Add a pass of Delaunay refinement.
//
// In that pass, the sizing field, of the refinement process of the
// triangulation, is based on the result of a poisson function with a
// triangulation, is based on the result of a Poisson function with a
// sample of the input points. The ratio is 'approximation_ratio'.
//
// For optimization reasons, the cell criteria of the refinement

View File

@ -17,7 +17,7 @@
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Poisson_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/facets_in_complex_3_to_triangle_mesh.h>
#include <CGAL/Poisson_reconstruction_function.h>
@ -99,7 +99,7 @@ namespace CGAL {
typedef typename Kernel::FT FT;
typedef CGAL::Poisson_reconstruction_function<Kernel> Poisson_reconstruction_function;
typedef CGAL::Labeled_mesh_domain_3<Kernel> Mesh_domain;
typedef CGAL::Poisson_mesh_domain_3<Kernel> Mesh_domain;
typedef typename CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Sequential_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
@ -115,7 +115,7 @@ namespace CGAL {
FT sm_sphere_radius = 5.0 * radius;
FT sm_dichotomy_error = sm_distance * spacing / 1000.0;
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(function, Sphere(inner_point, sm_sphere_radius),
Mesh_domain domain = Mesh_domain::create_Poisson_mesh_domain(function, Sphere(inner_point, sm_sphere_radius),
CGAL::parameters::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
Mesh_criteria criteria(CGAL::parameters::facet_angle = sm_angle,
@ -134,7 +134,7 @@ namespace CGAL {
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
turn_tag_into_mesh_3_manifold_option(tag)
.no_exude().no_perturb()
.surface_only()
.manifold_with_boundary());
const auto& tr = c3t3.triangulation();

View File

@ -5,6 +5,8 @@ project(Poisson_surface_reconstruction_3_Tests)
# Find CGAL
find_package(CGAL REQUIRED)
find_package(TBB QUIET)
include(CGAL_TBB_support)
# VisualC++ optimization for applications dealing with large data
if(MSVC)
@ -22,15 +24,25 @@ find_package(Eigen3 3.1.0 QUIET) #(requires 3.1.0 or greater)
include(CGAL_Eigen3_support)
if(TARGET CGAL::Eigen3_support)
# Executables that require Eigen 3.1
create_single_source_cgal_program("poisson_reconstruction_test_surface_mesher.cpp")
target_link_libraries(poisson_reconstruction_test_surface_mesher PRIVATE CGAL::Eigen3_support)
if(TARGET CGAL::TBB_support)
create_single_source_cgal_program("poisson_reconstruction_test_surface_mesher.cpp")
target_link_libraries(poisson_reconstruction_test_surface_mesher PRIVATE CGAL::Eigen3_support CGAL::TBB_support)
create_single_source_cgal_program("poisson_reconstruction_test_mesh_3.cpp")
target_link_libraries(poisson_reconstruction_test_mesh_3 PRIVATE CGAL::Eigen3_support)
create_single_source_cgal_program("poisson_reconstruction_test_mesh_3.cpp")
target_link_libraries(poisson_reconstruction_test_mesh_3 PRIVATE CGAL::Eigen3_support CGAL::TBB_support)
create_single_source_cgal_program("compare_mesh_3_vs_Poisson_implicit_surface_3.cpp")
target_link_libraries(compare_mesh_3_vs_Poisson_implicit_surface_3 PRIVATE CGAL::Eigen3_support)
create_single_source_cgal_program("compare_mesh_3_vs_Poisson_implicit_surface_3.cpp")
target_link_libraries(compare_mesh_3_vs_Poisson_implicit_surface_3 PRIVATE CGAL::Eigen3_support CGAL::TBB_support)
else()
create_single_source_cgal_program("poisson_reconstruction_test_surface_mesher.cpp")
target_link_libraries(poisson_reconstruction_test_surface_mesher PRIVATE CGAL::Eigen3_support)
create_single_source_cgal_program("poisson_reconstruction_test_mesh_3.cpp")
target_link_libraries(poisson_reconstruction_test_mesh_3 PRIVATE CGAL::Eigen3_support)
create_single_source_cgal_program("compare_mesh_3_vs_Poisson_implicit_surface_3.cpp")
target_link_libraries(compare_mesh_3_vs_Poisson_implicit_surface_3 PRIVATE CGAL::Eigen3_support)
endif()
else()
message("NOTICE: Tests in this directory require Eigen 3.1 (or greater), and will not be compiled.")
endif()

View File

@ -14,6 +14,7 @@
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Poisson_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/facets_in_complex_3_to_triangle_mesh.h>
@ -61,7 +62,12 @@ typedef CGAL::Surface_mesh_default_triangulation_3 STr;
typedef CGAL::Surface_mesh_complex_2_in_triangulation_3<STr> C2t3;
// Mesh_3
typedef CGAL::Labeled_mesh_domain_3<Kernel> Mesh_domain;
typedef CGAL::Labeled_mesh_domain_3<Kernel> LMesh_domain;
typedef typename CGAL::Mesh_triangulation_3<LMesh_domain>::type LTr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<LTr> LC3t3;
typedef CGAL::Mesh_criteria_3<LTr> LMesh_criteria;
typedef CGAL::Poisson_mesh_domain_3<Kernel> Mesh_domain;
typedef typename CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
@ -290,10 +296,10 @@ int main(int argc, char * argv[])
}
#endif
task_timer.stop();
// Prints status
std::cerr << "Total implicit function (triangulation+refinement+solver): " << task_timer.time() << " seconds\n";
task_timer.reset();
//***************************************
// Surface mesh generation
@ -331,47 +337,88 @@ int main(int argc, char * argv[])
const double implicit_function_time = reconstruction_timer.time();
reconstruction_timer.reset();
// MESH_3
// MESH_3 labeled
{
// Defines generation criteria
LMesh_criteria criteria(CGAL::parameters::facet_angle = fangle, CGAL::parameters::facet_size = fsize,
CGAL::parameters::facet_distance = fdist);
std::cout << "* Use Mesh_3 implicit *" << std::endl;
CGAL::Real_timer meshing_timer;
meshing_timer.start();
std::cout << "* Use Mesh_3 *" << std::endl;
// Defines generation criteria
Mesh_criteria criteria(CGAL::parameters::facet_angle = fangle,
CGAL::parameters::facet_size = fsize,
CGAL::parameters::facet_distance = fdist);
// Defines mesh domain
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(function, bsphere,
CGAL::parameters::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
LMesh_domain domain = LMesh_domain::create_implicit_mesh_domain(
function, bsphere, CGAL::parameters::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
// Generates mesh with manifold option
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
CGAL::parameters::no_exude().no_perturb()
.manifold_with_boundary());
LC3t3 c3t3 =
CGAL::make_mesh_3<LC3t3>(domain, criteria, CGAL::parameters::surface_only().manifold_with_boundary());
meshing_timer.stop();
const Tr& tr = c3t3.triangulation();
const LTr& tr = c3t3.triangulation();
// Prints status
std::cerr << "Mesh_3 meshing: " << meshing_timer.time() << " seconds, "
<< tr.number_of_vertices() << " output vertices"
<< std::endl;
std::cerr << "Mesh_3 meshing: " << meshing_timer.time() << " seconds, " << tr.number_of_vertices()
<< " output vertices" << std::endl;
if (tr.number_of_vertices() == 0)
if(tr.number_of_vertices() == 0)
return EXIT_FAILURE;
// Prints total reconstruction duration
reconstruction_timer.stop();
std::cerr << "Total reconstruction (implicit function + meshing): "
<< (implicit_function_time + reconstruction_timer.time()) << " seconds\n";
reconstruction_timer.reset();
<< (implicit_function_time + meshing_timer.time()) << " seconds\n";
// Converts to polyhedron
Polyhedron output_mesh;
CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, output_mesh);
std::ofstream out(output_basename + "_mesh_3.off");
std::ofstream out(output_basename + "_implicit_mesh_3.off");
out << output_mesh;
out.close();
}
// MESH_3 poisson
{
// Defines generation criteria
Mesh_criteria criteria(CGAL::parameters::facet_angle = fangle, CGAL::parameters::facet_size = fsize,
CGAL::parameters::facet_distance = fdist);
std::cout << "* Use Mesh_3 poisson *" << std::endl;
CGAL::Real_timer meshing_timer;
meshing_timer.start();
// Defines generation criteria
// Defines mesh domain
Mesh_domain domain = Mesh_domain::create_Poisson_mesh_domain(
function, bsphere, CGAL::parameters::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
// Generates mesh with manifold option
C3t3 c3t3 =
CGAL::make_mesh_3<C3t3>(domain, criteria, CGAL::parameters::surface_only().manifold_with_boundary());
meshing_timer.stop();
const Tr& tr = c3t3.triangulation();
// Prints status
std::cerr << "Mesh_3 meshing: " << meshing_timer.time() << " seconds, " << tr.number_of_vertices()
<< " output vertices" << std::endl;
if(tr.number_of_vertices() == 0)
return EXIT_FAILURE;
// Prints total reconstruction duration
std::cerr << "Total reconstruction (poisson function + meshing): "
<< (implicit_function_time + meshing_timer.time()) << " seconds\n";
// Converts to polyhedron
Polyhedron output_mesh;
CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, output_mesh);
std::ofstream out(output_basename + "_poisson_mesh_3.off");
out << output_mesh;
out.close();
}
@ -380,7 +427,6 @@ int main(int argc, char * argv[])
{
CGAL::Real_timer meshing_timer;
meshing_timer.start();
reconstruction_timer.start();
std::cout << "\n\n* Use Surface_mesher with Poisson_implicit_surface_3 *" << std::endl;
Surface_3 surface(function,
@ -408,9 +454,8 @@ int main(int argc, char * argv[])
return EXIT_FAILURE;
// Prints total reconstruction duration
reconstruction_timer.stop();
std::cerr << "Total reconstruction (implicit function + meshing): "
<< (implicit_function_time + reconstruction_timer.time()) << " seconds\n";
<< (implicit_function_time + meshing_timer.time()) << " seconds\n";
Polyhedron output_mesh;
CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, output_mesh);

View File

@ -18,6 +18,7 @@
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Poisson_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/facets_in_complex_3_to_triangle_mesh.h>
@ -57,12 +58,18 @@ typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
// Poisson implicit function
typedef CGAL::Poisson_reconstruction_function<Kernel> Poisson_reconstruction_function;
// Mesh_3
// Mesh_3 using Labeled_mesh_domain_3
typedef CGAL::Labeled_mesh_domain_3<Kernel> Mesh_domain;
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
// Mesh_3 using Poisson_mesh_domain_3
typedef CGAL::Poisson_mesh_domain_3<Kernel> Poisson_mesh_domain;
typedef CGAL::Mesh_triangulation_3<Poisson_mesh_domain>::type Poisson_Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Poisson_Tr> Poisson_C3t3;
typedef CGAL::Mesh_criteria_3<Poisson_Tr> Poisson_mesh_criteria;
namespace params = CGAL::parameters;
// ----------------------------------------------------------------------------
@ -171,7 +178,6 @@ int main(int argc, char * argv[])
<< task_timer.time() << " seconds, "
<< (memory>>20) << " Mb allocated"
<< std::endl;
task_timer.reset();
//***************************************
// Checks requirements
@ -192,14 +198,14 @@ int main(int argc, char * argv[])
continue;
}
CGAL::Timer reconstruction_timer; reconstruction_timer.start();
//***************************************
// Computes implicit function
//***************************************
std::cerr << "Computes Poisson implicit function...\n";
task_timer.reset();
// Creates implicit function from the read points.
// Note: this method requires an iterator over points
// + property maps to access each point's position and normal.
@ -218,9 +224,9 @@ int main(int argc, char * argv[])
continue;
}
double poisson_time = task_timer.time();
// Prints status
std::cerr << "Total implicit function (triangulation+refinement+solver): " << task_timer.time() << " seconds\n";
task_timer.reset();
std::cerr << "Total implicit function (triangulation+refinement+solver): " << poisson_time << " seconds\n";
//***************************************
// Surface mesh generation
@ -250,45 +256,92 @@ int main(int argc, char * argv[])
FT sm_sphere_radius = 5.0 * radius;
FT sm_dichotomy_error = sm_distance*average_spacing/1000.0; // Dichotomy error must be << sm_distance
// Defines surface mesh generation criteria
Mesh_criteria criteria(params::facet_angle = sm_angle,
params::facet_size = sm_radius*average_spacing,
params::facet_distance = sm_distance*average_spacing);
{
task_timer.reset();
// Defines surface mesh generation criteria
Mesh_criteria criteria(params::facet_angle = sm_angle, params::facet_size = sm_radius * average_spacing,
params::facet_distance = sm_distance * average_spacing);
std::cerr << " make_mesh_3 with sphere center=("<<inner_point << "),\n"
<< " sphere radius="<<sm_sphere_radius<<",\n"
<< " angle="<<sm_angle << " degrees,\n"
<< " triangle size="<<sm_radius<<" * average spacing="<<sm_radius*average_spacing<<",\n"
<< " distance="<<sm_distance<<" * average spacing="<<sm_distance*average_spacing<<",\n"
<< " dichotomy = distance/"<<sm_distance*average_spacing/sm_dichotomy_error<<",\n"
<< " manifold_with_boundary()\n";
std::cerr << " make_mesh_3 with sphere center=(" << inner_point << "),\n"
<< " sphere radius=" << sm_sphere_radius << ",\n"
<< " angle=" << sm_angle << " degrees,\n"
<< " triangle size=" << sm_radius
<< " * average spacing=" << sm_radius * average_spacing << ",\n"
<< " distance=" << sm_distance
<< " * average spacing=" << sm_distance * average_spacing << ",\n"
<< " dichotomy = distance/" << sm_distance * average_spacing / sm_dichotomy_error
<< ",\n"
<< " manifold_with_boundary()\n";
// Generates surface mesh with manifold option
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(function, bsphere,
params::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::no_exude().no_perturb().manifold_with_boundary());
// Generates surface mesh with manifold option
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(
function, bsphere, params::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::surface_only().manifold_with_boundary());
// Prints status
/*long*/ memory = CGAL::Memory_sizer().virtual_size();
const Tr& tr = c3t3.triangulation();
std::cerr << "Surface meshing: " << task_timer.time() << " seconds, "
<< tr.number_of_vertices() << " output vertices, "
<< (memory>>20) << " Mb allocated"
<< std::endl;
task_timer.reset();
// Prints status
/*long*/ memory = CGAL::Memory_sizer().virtual_size();
const Tr& tr = c3t3.triangulation();
std::cerr << "Surface meshing: " << task_timer.time() << " seconds, " << tr.number_of_vertices()
<< " output vertices, " << (memory >> 20) << " Mb allocated" << std::endl;
task_timer.reset();
if(tr.number_of_vertices() == 0) {
accumulated_fatal_err = EXIT_FAILURE;
continue;
if(tr.number_of_vertices() == 0) {
accumulated_fatal_err = EXIT_FAILURE;
continue;
}
// Converts to polyhedron
Polyhedron output_mesh;
CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, output_mesh);
// Prints total reconstruction duration
std::cerr << "Total reconstruction using Labeled_mesh_domain_3 (implicit function + meshing): "
<< poisson_time + task_timer.time() << " seconds\n";
}
// Converts to polyhedron
Polyhedron output_mesh;
CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, output_mesh);
{
task_timer.reset();
// Defines surface mesh generation criteria
Poisson_mesh_criteria criteria(params::facet_angle = sm_angle, params::facet_size = sm_radius * average_spacing,
params::facet_distance = sm_distance * average_spacing);
// Prints total reconstruction duration
std::cerr << "Total reconstruction (implicit function + meshing): " << reconstruction_timer.time() << " seconds\n";
std::cerr << " make_mesh_3 with sphere center=(" << inner_point << "),\n"
<< " sphere radius=" << sm_sphere_radius << ",\n"
<< " angle=" << sm_angle << " degrees,\n"
<< " triangle size=" << sm_radius
<< " * average spacing=" << sm_radius * average_spacing << ",\n"
<< " distance=" << sm_distance
<< " * average spacing=" << sm_distance * average_spacing << ",\n"
<< " dichotomy = distance/" << sm_distance * average_spacing / sm_dichotomy_error
<< ",\n"
<< " manifold_with_boundary()\n";
// Generates surface mesh with manifold option
Poisson_mesh_domain domain = Poisson_mesh_domain::create_Poisson_mesh_domain(
function, bsphere, params::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));
Poisson_C3t3 c3t3 =
CGAL::make_mesh_3<Poisson_C3t3>(domain, criteria, params::surface_only().manifold_with_boundary());
// Prints status
/* long */ memory = CGAL::Memory_sizer().virtual_size();
const Poisson_Tr& tr = c3t3.triangulation();
std::cerr << "Surface meshing: " << task_timer.time() << " seconds, " << tr.number_of_vertices()
<< " output vertices, " << (memory >> 20) << " Mb allocated" << std::endl;
task_timer.reset();
if(tr.number_of_vertices() == 0) {
accumulated_fatal_err = EXIT_FAILURE;
continue;
}
// Converts to polyhedron
Polyhedron output_mesh;
CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, output_mesh);
// Prints total reconstruction duration
std::cerr << "Total reconstruction using Poisson_mesh_domain_3 (implicit function + meshing): "
<< poisson_time + task_timer.time() << " seconds\n";
}
} // for each input file