Pierre Alliez, Laurent Rineau, Stéphane Tayeb, Jane Tournois, Mariette Yvinec
This package is devoted to the generation of isotropic simplicial meshes discretizing 3D domains. The domain to be meshed is a region of 3D space that has to be bounded. The region may be connected or composed of multiple components and/or subdivided in several subdomains.
Boundary and subdivision surfaces are either smooth or piecewise smooth surfaces, formed with planar or curved surface patches. Surfaces may exhibit 1dimensional features (e.g. crease edges) and 0dimensional features (e.g. singular points as corners tips, cusps or darts), that have to be fairly approximated in the mesh.
The output mesh is a 3dimensional triangulation, including subcomplexes that approximate each input domain feature: subdomain, boundary surface patch or input domain feature with dimension 0 or 1. Thus, the output mesh includes a 3D submesh covering each subdomain, a surface mesh approximating each boundary or subdividing surface patch, a polyline approximation for each 1dimensional feature and of course a vertex on each corner.
The main entry points of the package are two global functions that respectively generate and refine such meshes. The mesh generator is customized to output a mesh that fits as much as possible the user needs, for instance in terms of sizing field or with respect to some user customized quality criteria.
The meshing engine used in this mesh generator is based on Delaunay refinement [Che93, Rup95, She98b]. It uses the notion of restricted Delaunay triangulation to approximate 1dimensional curve segments and surface patches [BO05]. Before the refinement, a mechanism of protecting balls is set up on 1dimensional features, if any, to ensure a fair representation of those features in the mesh, and also to guarantee the termination of the refinement process, whatever may be the input geometry, in particular whatever small angles the boundary and subdivision surface patches may form [CDL07, CDR07]. The Delaunay refinement is followed by a mesh optimization phase to remove slivers and provide a good quality mesh.
The domain to be meshed is assumed to be bounded and representable as a pure 3D complex. A 3D complex is a set of faces with dimension 0, 1, 2 and 3 such that all faces are pairwise interior disjoint, and the boundary of each face of the complex is the union of faces of the complex. The 3D complex is pure, meaning that each face is included in a face of dimension 3, so that the complex is entirely described by the set of its 3D faces and their subfaces. However the 3D complex needs not be connected. The set of faces with dimension lower or equal than 2 forms a 2D subcomplex which needs not be manifold, neither pure, nor connected: some 3D faces may have dangling 2D or 1D faces in their boundary faces.
In the rest of the documentation, we will refer to the input 3D complex as the input domain. The faces of the input domain with dimension 0, 1, 2 and 3 are called respectively corners, curve segments, surface patches and subdomains to clearly distinguish them from the faces of the mesh that are called vertices, edges, facets and cells.
Note that the input complex faces are not required to be linear nor smooth. Surface patches, for instance, may be smooth surface patches, or portions of surface meshes with boundaries. Curve segments may be for instance straight segments, curved segments or polylines. Each of those features will be accurately represented in the final mesh.
The 0 and 1dimensional features of the input domain are usually singular points of the subdomain boundaries, however this is not required. Furthermore those features are not required to cover all the subdomains boundaries singularities but only those that need to be accurately represented in the final mesh. In the following, we say that a domain has features when it has 0 and 1dimensional features that need to be accurately represented in the mesh, and we call those features exposed features. Therefore, a domain may be without features either because all boundary surface patches are smooth closed surfaces, or simply because the curves joining different surface patches and the singularities of those patches need not be accurately approximated in the final mesh.
Note also that input complex faces are not required to be connected. Faces of the input domain are identified by indexes. If a subdomain is not connected, its different components receive the same index. Likewise different surface patches, segment curves or corners may share the same index. Each connected component of a feature will be accurately represented in the final mesh. Note however that the occurrence of multiply connected faces in the input complex may affect the relevance of internal topological checks performed by the mesh generator. Also the mesh generator will not be able to apply different meshing criteria, e.g. different sizing field, for the different connected components of a single feature.
The domain is input to the mesh generation function, as a domain class, often called the oracle, that provides predicates and constructors related to the domain, the subdomains, the boundary surface patches and also the 0 and 1dimensional exposed features, if any. Mainly, the oracle provides a predicate to test if a given query point belongs to the domain or not and to find in which subdomain it lies in the affirmative case. The domain class also provides predicates and constructors to test the intersection of a query line segment with the boundary surface patches and to build some intersection points if any. Lastly, if the input domain includes 1dimensional exposed features, the domain class provides a way to construct sample points on these features.
The current implementation provides classes to represent domains bounded by isosurfaces of implicit functions, polyhedral domains and domains defined through 3D labeled images. Currently, 1dimensional features may be defined as segments and polyline segments.
The resulting mesh is output as a subcomplex of a 3D Delaunay triangulation, in a class providing various iterators on mesh elements.
The 3D triangulation provides approximations of the subdomains, surface patches and curve segments and corners, according to the restricted Delaunay triangulation paradigm. This means that each subdomain is approximated by the union of the tetrahedral cells whose circumcenters are located inside the domain (or subdomain). Each surface patch is approximated by the union of the Delaunay mesh facets whose dual Voronoi edges intersect the surface patch. Such mesh facets are called surface facets in the following. The 1dimensional exposed features are approximated by sequences of mesh edges and the 0dimensional exposed features are represented by mesh vertices.
The mesh generation algorithm is mainly a Delaunay refinement process. The Delaunay refinement is preceded by a protecting phase to ensure an accurate representation of 1dimensional features if any, and followed by an optimization phase to achieve a good quality mesh.
The Delaunay refinement process is driven by criteria concerning either the size and shape of mesh cells and surface facets. The refinement process terminates when there are no more mesh cells or surface facets violating the criteria.
The criteria are designed to achieve a nice spread of the mesh vertices while ensuring the termination of the refinement process. Those criteria may be somehow tuned to the user needs to achieve for instance the respect of a sizing field by mesh elements, some topological conditions on the representation of boundary surfaces in the mesh, and/or some error bound for the approximation of boundary surfaces. To some extend, the user may tune the Delaunay refinement to a prescribed tradeoff between mesh quality and mesh density. The mesh density refers to the number of mesh vertices and cells, i.e. to the complexity of the mesh. The mesh quality referred to here is measured by the radius edge ratio of surface facets end mesh cells, where the radius edge ratio of a simplex (triangle or tetrahedron) is the the ratio between its circumradius and its shortest edge.
If the domain description includes 0 dimensional features, the corresponding points are inserted into the Delaunay triangulation from the start.
If the domain has 1dimensional exposed features, the method of protecting balls [CDR07, CDL07] is used to achieve an accurate representation of those features in the mesh and to guarantee that the refinement process terminates whatever may be the dihedral angles formed by input surface patches incident to a given 1feature or the angles formed by two 1features incident to a 0feature.
According to this method, the 1dimensional features are
sampled with points and covered by protecting balls centered on the sample points,
in such a way that :
no three balls intersect
no pair of balls centered on different 1features intersect.
The triangulation embedding the mesh is in fact a weighted Delaunay triangulation, and the triangulation is initialized by the insertion of all the protecting balls, regarded as weighted points. The Delaunay refinement process is then launched as before except that refinement points are no longer circumcenters but are weighted circumcenters. All Steiner vertices inserted by the refinement process are given a zero weight.
The method guarantees:
1) that each segment joining
two successive centers on a 1dimensional feature will stay in the triangulation,
thus ensuring an accurate approximation of the 1dimensional features.
2) that the refinement process will never try to insert a refinement point in the union of the
protecting balls, which ensures the termination of the refinement process.
Any tetrahedron that is quasi degenerate has a big radius edge ratio except those belonging to the family of slivers. A sliver is easily obtained as the convex hull of 4 points close to the equatorial circle of a 3D ball and roughly equally spread along this circle. The Delaunay refinement tracks tetrahedra with big radius edge ratio and therefore eliminates all kinds of badly shaped tetrahedra except slivers.
Therefore, at the end of the refinement process, some sliver shaped tetrahedra may occur in the mesh. The optimization phase aims at eliminating slivers.
The optimization phase is a succession of optimization processes, including possibly a Lloyd smoothing, an odtsmoothing, a perturber and an exuder.
The Lloyd and odtsmoother are global optimizers moving the mesh vertices to minimize a mesh energy. Those optimizers are described respectively in [DFG99, DW02] and in [Che04, ACSYD05]. In both cases the mesh energy is the L1 error resulting from the interpolation of the function f(x) =x^{2} by a piecewise linear function. In the case of the Lloyd smoother, the interpolation is linear in each Voronoi cell of the set of mesh vertices. In the case of the odtsmoother, the interpolation is linear in each cell of the Delaunay triangulation of the mesh vertices, hence the name odt which is an abbreviation for ``optimal Delaunay triangulation''.
The Lloyd optimizer is known to be blind to the occurrence of slivers in the mesh while the odtsmoother tends to chase them out. Both of them are global optimizers, meaning that they try to improve the whole mesh rather than focusing on the worst elements. However, both are empirically known to be very efficient as a preliminary step of optimization, as they tend to enhance the efficiency of the perturber and/or exuder applied next, see Figure 50.3
The perturber and the exuder focus on improving the worst mesh elements. The perturber [TSA09] improves the meshes by local changes in the vertices positions aiming to make sliver disappear. The exuder [CDE^{+}00] chases the remaining slivers by reweighting mesh vertices with optimal weights.
Each optimization process can be activated or not, and tuned according to the user requirements and the available time. By default, only the perturber and the exuder are activated.
Optimization processes are designed to improve mesh quality. However, beware that such an improvement is obtained by perturbing mesh vertices and modifying the mesh connectivity which has an impact on the strict compliance to the refinement criteria. Though a strict compliance to mesh criteria is granted at the end of the Delaunay refinement, this may no longer be true after some optimization processes. Also beware that the default behavior does involve some optimization processes.
 

 
 


The function make_mesh_3 generates from scratch a mesh of the input domain, while the function refine_mesh_3 refines an existing mesh of the input domain. Note that as the protection of 0 and 1dimensional features does not rely on Delaunay refinement, the function refine_mesh_3 has no parameter to preserve features.
This concept provides, among others, member functions to test whether or not a query segment intersects boundary surfaces, and to compute an intersection point in the affirmative. The MeshDomain_3 concept adds member functions which given a query point tell whether the point lies inside or outside the domain and in which subdomain the point lies if inside.
If the domain description includes 0 and 1dimensional features that have to be accurately represented in the final mesh, the template parameter MeshDomain_3 is required to be of a model of the refined concept MeshDomainWithFeatures_3. Mainly the concept MeshDomainWithFeatures_3 provides the incidence graph of 0, 1 and 2dimensional features, and a member function to construct sample points on curve segments.
Using the parameter of type Features, the user whose domain is a model of MeshDomainWithFeatures_3 can choose to have the corners and curve segments of the domain represented in the mesh or not. The type Features of this parameter is an internal undescribed type. The library provides functions to construct appropriate values of that type.
The criteria for surface facets are governed by the four following parameters:
The criteria for mesh cells are governed by two parameters:
Figure 50.4 shows how the mesh generation process behaves with respect to these parameters.
If the domain has 1dimensional exposed features, the criteria includes a sizing field to guide the sampling of 1dimensional features with protecting balls centers.
The four additional parameters are optimization parameters. They control which optimization processes are performed and allow the user to tune the parameters of the activated optimization processes. These parameters have internal types which are not described but the library provides global functions to generate appropriate values of these types:
These parameters are optional and can be passed in any order. If one parameter is not passed the default value is used. By default, only the perturber and the exuder are activated. Note that whatever may be the optimization processes activated by make_mesh_3 or refine_mesh_3, they are always launched in the order that is a suborder of the following: odt smoother, Lloyd smoother, perturber and exuder.
The package also provides four global functions to launch each optimization process independently. These functions are useful for advanced experimentation on the efficiency of each optimization method. Note however that the exuder adds on mesh vertices weights that are conditioned by vertices positions. Therefore an exudation process should never be run before a smoother or a perturber. For a maximum efficiency, whatever may be the optimization processes activated, they should be launched in the order that is a suborder of the following: odtsmoother, Lloydsmoother, perturber, exuder.
 

 
 

 
 

 
 


Note that the global functions activating the optimization processes or launching those processes have themselves parameters (see details in reference pages) to tune the optimization process.
Note the use of named parameters (from Boost library) in the constructor of the Mesh_criteria instance.
File: examples/Mesh_3/mesh_implicit_sphere.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Implicit_mesh_domain_3.h> #include <CGAL/make_mesh_3.h> // Domain typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef K::FT FT; typedef K::Point_3 Point; typedef FT (Function)(const Point&); typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain; // Triangulation typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3; // Criteria typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria; // To avoid verbose function and named parameters call using namespace CGAL::parameters; // Function FT sphere_function (const Point& p) { return CGAL::squared_distance(p, Point(CGAL::ORIGIN))1; } int main() { // Domain (Warning: Sphere_3 constructor uses squared radius !) Mesh_domain domain(sphere_function, K::Sphere_3(CGAL::ORIGIN, 2.)); // Mesh criteria Mesh_criteria criteria(facet_angle=30, facet_size=0.1, facet_distance=0.025, cell_radius_edge=2, cell_size=0.1); // Mesh generation C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria); // Output std::ofstream medit_file("out.mesh"); c3t3.output_to_medit(medit_file); return 0; }
File: examples/Mesh_3/mesh_polyhedral_domain.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Polyhedral_mesh_domain_3.h> #include <CGAL/make_mesh_3.h> #include <CGAL/refine_mesh_3.h> // IO #include <CGAL/IO/Polyhedron_iostream.h> // Domain typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Polyhedron_3<K> Polyhedron; typedef CGAL::Polyhedral_mesh_domain_3<Polyhedron, K> Mesh_domain; // Triangulation typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3; // Criteria typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria; // To avoid verbose function and named parameters call using namespace CGAL::parameters; int main() { // Create input polyhedron Polyhedron polyhedron; std::ifstream input("data/elephant.off"); input >> polyhedron; // Create domain Mesh_domain domain(polyhedron); // Mesh criteria (no cell_size set) Mesh_criteria criteria(facet_angle=25, facet_size=0.15, facet_distance=0.008, cell_radius_edge=3); // Mesh generation C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), no_exude()); // Output std::ofstream medit_file("out_1.mesh"); c3t3.output_to_medit(medit_file); medit_file.close(); // Set tetrahedron size (keep cell_radius_edge), ignore facets Mesh_criteria new_criteria(cell_radius_edge=3, cell_size=0.03); // Mesh refinement CGAL::refine_mesh_3(c3t3, domain, new_criteria); // Output medit_file.open("out_2.mesh"); c3t3.output_to_medit(medit_file); return 0; }
In the following example, the image is read from the file liner.inr.gz which is encoded in the format of the library Inrimage http://inrimage.gforge.inria.fr/. The resulting mesh is shown in Figure 50.7.
File: examples/Mesh_3/mesh_3D_image.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Labeled_image_mesh_domain_3.h> #include <CGAL/make_mesh_3.h> #include <CGAL/Image_3.h> // Domain typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Labeled_image_mesh_domain_3<CGAL::Image_3,K> Mesh_domain; // Triangulation typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3; // Criteria typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria; // To avoid verbose function and named parameters call using namespace CGAL::parameters; int main() { // Loads image CGAL::Image_3 image; image.read("data/liver.inr.gz"); // Domain Mesh_domain domain(image); // Mesh criteria Mesh_criteria criteria(facet_angle=30, facet_size=6, facet_distance=4, cell_radius_edge=3, cell_size=8); // Meshing C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria); // Output std::ofstream medit_file("out.mesh"); c3t3.output_to_medit(medit_file); return 0; }
The following example shows how to use an analytical function as sizing field.
File: examples/Mesh_3/mesh_implicit_sphere_variable_size.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Implicit_mesh_domain_3.h> #include <CGAL/make_mesh_3.h> // Domain typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef K::FT FT; typedef K::Point_3 Point; typedef FT (Function)(const Point&); typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain; // Triangulation typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3; // Criteria typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria; // Sizing field struct Spherical_sizing_field { typedef ::FT FT; typedef Point Point_3; typedef Mesh_domain::Index Index; FT operator()(const Point_3& p, const int, const Index&) const { FT sq_d_to_origin = CGAL::squared_distance(p, Point(CGAL::ORIGIN)); return CGAL::abs( CGAL::sqrt(sq_d_to_origin)0.5 ) / 5. + 0.025; } }; // To avoid verbose function and named parameters call using namespace CGAL::parameters; // Function FT sphere_function (const Point& p) { return CGAL::squared_distance(p, Point(CGAL::ORIGIN))1; } int main() { // Domain (Warning: Sphere_3 constructor uses squared radius !) Mesh_domain domain(sphere_function, K::Sphere_3(CGAL::ORIGIN, 2.)); // Mesh criteria Spherical_sizing_field size; Mesh_criteria criteria(facet_angle=30, facet_size=0.1, facet_distance=0.025, cell_radius_edge_ratio=2, cell_size=size); // Mesh generation C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_exude(), no_perturb()); // Output std::ofstream medit_file("out.mesh"); c3t3.output_to_medit(medit_file); return 0; }
The following example shows how to use different size for different organs in a 3D medical image.
File: examples/Mesh_3/mesh_3D_image_variable_size.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Mesh_constant_domain_field_3.h> #include <CGAL/Labeled_image_mesh_domain_3.h> #include <CGAL/make_mesh_3.h> #include <CGAL/Image_3.h> // Domain typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Labeled_image_mesh_domain_3<CGAL::Image_3,K> Mesh_domain; // Triangulation typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3; // Criteria typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria; typedef CGAL::Mesh_constant_domain_field_3<Mesh_domain::R, Mesh_domain::Index> Sizing_field; // To avoid verbose function and named parameters call using namespace CGAL::parameters; int main() { // Loads image CGAL::Image_3 image; image.read("data/liver.inr.gz"); // Domain Mesh_domain domain(image); // Sizing field: set global size to 8 and kidney size (label 127) to 3 double kidney_size = 3.; int volume_dimension = 3; Sizing_field size(8); size.set_size(kidney_size, volume_dimension, domain.index_from_subdomain_index(127)); // Mesh criteria Mesh_criteria criteria(facet_angle=30, facet_size=6, facet_distance=2, cell_radius_edge_ratio=3, cell_size=size); // Meshing C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria); // Output std::ofstream medit_file("out.mesh"); c3t3.output_to_medit(medit_file); return 0; }
File: examples/Mesh_3/mesh_polyhedral_domain_with_features.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Delaunay_triangulation_3.h> #include <CGAL/Triangulation_vertex_base_with_info_3.h> #include <CGAL/Mesh_3/Robust_intersection_traits_3.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Polyhedral_mesh_domain_3.h> #include <CGAL/make_mesh_3.h> #include <CGAL/refine_mesh_3.h> #include <CGAL/IO/Polyhedron_iostream.h> #include <CGAL/Mesh_domain_with_polyline_features_3.h> typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Mesh_3::Robust_intersection_traits_3<K> Geom_traits; typedef Geom_traits::Point_3 Point_3; typedef CGAL::Polyhedron_3<Geom_traits> Polyhedron; typedef CGAL::Mesh_domain_with_polyline_features_3< CGAL::Polyhedral_mesh_domain_3<Polyhedron, Geom_traits> > MD_IntersectedLines; typedef CGAL::Mesh_triangulation_3<MD_IntersectedLines>::type Tr_MDI; typedef CGAL::Mesh_criteria_3<Tr_MDI> Mesh_criteria_MDI; typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr_MDI> C3t3_MDI; typedef C3t3_MDI::Facets_in_complex_iterator Facet_Iter_inComplex; typedef std::vector<Point_3> InterPolyline; typedef std::list<InterPolyline> InterPolylines; using namespace CGAL::parameters; int main() { //constructing the bounding box by triangles Polyhedron bounding_box; //top face bounding_box.make_triangle(Point_3(100,100,100), Point_3(0,100,100), Point_3(0,0,100)); bounding_box.make_triangle(Point_3(100,100,100), Point_3(0,0,100), Point_3(100,0,100)); //bottom face bounding_box.make_triangle(Point_3(100,100,0), Point_3(100,0,0), Point_3(0,0,0)); bounding_box.make_triangle(Point_3(100,100,0), Point_3(0,0,0), Point_3(0,100,0)); //front face bounding_box.make_triangle(Point_3(0,0,0), Point_3(100,0,100), Point_3(0,0,100)); bounding_box.make_triangle(Point_3(0,0,0), Point_3(100,0,0), Point_3(100,0,100)); //back face bounding_box.make_triangle(Point_3(0,100,0), Point_3(0,100,100), Point_3(100,100,100)); bounding_box.make_triangle(Point_3(100,100,100), Point_3(100,100,0), Point_3(0,100,0)); //left face bounding_box.make_triangle(Point_3(0,0,0), Point_3(0,0,100), Point_3(0,100,100)); bounding_box.make_triangle(Point_3(0,0,0), Point_3(0,100,100), Point_3(0,100,0)); //right face bounding_box.make_triangle(Point_3(100,100,100), Point_3(100,0,0), Point_3(100,100,0)); bounding_box.make_triangle(Point_3(100,100,100), Point_3(100,0,100), Point_3(100,0,0)); //constructing the input polyhedral face std::vector<Polyhedron*> AllSurfaces; Polyhedron NewPoly_1; NewPoly_1.make_triangle(Point_3(20,20,20), Point_3(25,20,20), Point_3(25,40,20)); NewPoly_1.make_triangle(Point_3(20,20,20), Point_3(25,40,20), Point_3(20,40,20)); AllSurfaces.push_back(&NewPoly_1); Polyhedron NewPoly_2; NewPoly_2.make_triangle(Point_3(25,40,20), Point_3(25,20,20), Point_3(60,20,20)); NewPoly_2.make_triangle(Point_3(60,40,20), Point_3(25,40,20), Point_3(60,20,20)); AllSurfaces.push_back(&NewPoly_2); std::ofstream file; file.open("input_1.off"); file << NewPoly_1; file.close(); file.open("input_2.off"); file << NewPoly_2; file.close(); file.open("bbox.off"); file << bounding_box; file.close(); //constructing the line feature InterPolylines interPolySet; InterPolyline interLine; interLine.push_back(Point_3(25,40,20)); interLine.push_back(Point_3(25,20,20)); interPolySet.push_back(interLine); //generating the mesh MD_IntersectedLines domain(AllSurfaces.begin(), AllSurfaces.end(), bounding_box); domain.add_features(interPolySet.begin(), interPolySet.end()); Mesh_criteria_MDI criteria(edge_size = 5, facet_angle=20, cell_size=5, facet_topology = CGAL::FACET_VERTICES_ON_SURFACE); C3t3_MDI c3t3 = CGAL::make_mesh_3<C3t3_MDI>(domain, criteria, no_exude(), no_perturb()); //output std::ofstream medit_file("out_1.mesh"); c3t3.output_to_medit(medit_file); medit_file.close(); } #if 0 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Polyhedral_mesh_domain_with_features_3.h> #include <CGAL/make_mesh_3.h> // Domain typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain; // Triangulation typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3< Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_segment_index> C3t3; // Criteria typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria; // To avoid verbose function and named parameters call using namespace CGAL::parameters; int main() { // Create domain Mesh_domain domain("data/fandisk.off"); // Get sharp features domain.detect_features(); // Mesh criteria Mesh_criteria criteria(edge_size = 0.025, facet_angle = 25, facet_size = 0.05, facet_distance = 0.005, cell_radius_edge_ratio = 3, cell_size = 0.05); // Mesh generation C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria); // Output std::ofstream medit_file("out.mesh"); c3t3.output_to_medit(medit_file); } #endif
File: examples/Mesh_3/mesh_two_implicit_spheres_with_balls.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Implicit_mesh_domain_3.h> #include <CGAL/Mesh_domain_with_polyline_features_3.h> #include <CGAL/make_mesh_3.h> // Kernel typedef CGAL::Exact_predicates_inexact_constructions_kernel K; // Domain typedef K::FT FT; typedef K::Point_3 Point; typedef FT (Function)(const Point&); typedef CGAL::Mesh_domain_with_polyline_features_3< CGAL::Implicit_mesh_domain_3<Function,K> > Mesh_domain; // Polyline typedef std::vector<Point> Polyline; typedef std::list<Polyline> Polylines; // Triangulation typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3< Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_segment_index> C3t3; // Criteria typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria; // To avoid verbose function and named parameters call using namespace CGAL::parameters; // Function FT sphere_function1 (const Point& p) { return CGAL::squared_distance(p, Point(CGAL::ORIGIN))2; } FT sphere_function2 (const Point& p) { return CGAL::squared_distance(p, Point(2, 0, 0))2; } FT sphere_function (const Point& p) { if(sphere_function1(p) < 0  sphere_function2(p) < 0) return 1; else return 1; } #include <cmath> int main() { // Domain (Warning: Sphere_3 constructor uses squared radius !) Mesh_domain domain(sphere_function, K::Sphere_3(Point(1, 0, 0), 6.)); // Mesh criteria Mesh_criteria criteria(edge_size = 0.15, facet_angle = 25, facet_size = 0.15, cell_radius_edge_ratio = 2, cell_size = 0.15); // Create edge that we want to preserve Polylines polylines (1); Polyline& polyline = polylines.front(); for(int i = 0; i < 360; ++i) { Point p (1, std::cos(i*CGAL_PI/180), std::sin(i*CGAL_PI/180)); polyline.push_back(p); } polyline.push_back(polyline.front()); // close the line // Insert edge in domain domain.add_features(polylines.begin(), polylines.end()); // Mesh generation without feature preservation C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, CGAL::parameters::no_features()); std::ofstream medit_file("outnoprotection.mesh"); c3t3.output_to_medit(medit_file); medit_file.close(); c3t3.clear(); // Mesh generation with feature preservation c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria); // Output medit_file.open("outwithprotection.mesh"); c3t3.output_to_medit(medit_file); medit_file.close(); return 0; }
In the previous examples, the mesh generation is launched through a call make_mesh_3(domain,criteria) with a minimal number of parameters. In such cases, the default optimization strategy is applied: after the Delaunay refinement process two optimization steps are performed, a perturbation and a sliver exudation. The following examples show how to disable default optimization steps and how to tune the parameters of optimization steps.
In this first example, we show how to disable the exudation step. The optimization phase after the refinement includes only a perturbation phase which is launched with no time bound and an objective of 10 degrees for the minimum dihedral angle of the mesh. The example shows two ways of achieving the same result. The first way issues a single call to make_mesh_3 with the required optimization process activated and tuned. In the second way, make_mesh_3 is first called without any optimization process and the resulting mesh is next optimized through a call to perturb_mesh_3 with tuned parameters.
File: examples/Mesh_3/mesh_optimization_example.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Labeled_image_mesh_domain_3.h> #include <CGAL/make_mesh_3.h> #include <CGAL/Image_3.h> // Domain typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Labeled_image_mesh_domain_3<CGAL::Image_3,K> Mesh_domain; // Triangulation typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3; // Mesh Criteria typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria; // To avoid verbose function and named parameters call using namespace CGAL::parameters; int main() { // Domain CGAL::Image_3 image; image.read("data/liver.inr.gz"); Mesh_domain domain(image); // Mesh criteria Mesh_criteria criteria(facet_angle=30, facet_size=5, facet_distance=1.5, cell_radius_edge=2, cell_size=7); // Mesh generation and optimization in one call (sliver_bound is the // targeted dihedral angle in degree) C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_exude(), perturb(sliver_bound=10, time_limit=15)); // Mesh generation and optimization in several call C3t3 c3t3_bis = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), no_exude()); CGAL::perturb_mesh_3(c3t3_bis, domain, time_limit=15); // Output std::ofstream medit_file("out.mesh"); c3t3.output_to_medit(medit_file); std::ofstream medit_file_bis("out_bis.mesh"); c3t3_bis.output_to_medit(medit_file_bis); return 0; }
In this second example, we show how to call the Lloyd optimization on the mesh, followed by a call to exudation. We set a time bound of 30s for the Lloyd optimization. We set a time bound of 10s and a sliver bound of 10 degrees for the exuder.
File: examples/Mesh_3/mesh_optimization_lloyd_example.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Mesh_triangulation_3.h> #include <CGAL/Mesh_complex_3_in_triangulation_3.h> #include <CGAL/Mesh_criteria_3.h> #include <CGAL/Labeled_image_mesh_domain_3.h> #include <CGAL/make_mesh_3.h> #include <CGAL/Image_3.h> // Domain typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Labeled_image_mesh_domain_3<CGAL::Image_3,K> Mesh_domain; // Triangulation typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3; // Mesh Criteria typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria; // To avoid verbose function and named parameters call using namespace CGAL::parameters; int main() { // Domain CGAL::Image_3 image; image.read("data/liver.inr.gz"); Mesh_domain domain(image); // Mesh criteria Mesh_criteria criteria(facet_angle=30, facet_distance=1.2, cell_radius_edge=2); // Mesh generation and optimization in one call C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, lloyd(time_limit=30), no_perturb(), exude(time_limit=10, sliver_bound=10)); // Mesh generation and optimization in several call C3t3 c3t3_bis = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), no_exude()); CGAL::lloyd_optimize_mesh_3(c3t3_bis, domain, time_limit=30); CGAL::exude_mesh_3(c3t3_bis, sliver_bound=10, time_limit=10); // Output std::ofstream medit_file("out.mesh"); c3t3.output_to_medit(medit_file); std::ofstream medit_file_bis("out_bis.mesh"); c3t3_bis.output_to_medit(medit_file_bis); return 0; }
We provide here some benchmarks of the performances of the mesh generation engine. The machine used is a PC running Linux64 with two Intel Xeon CPU X5450 clocked at 3.00 GHz with 32GB of RAM. The program has been compiled with g++ v4.3.2 with the O3 option. Note that our implementation does not take advantage of multicore architectures.
Those benchmarks have been done using Cgal v3.8.
We study the refinement part of the mesh generation engine in this section. We give the CPU time (measured by CGAL::Timer) using the 3 provided oracles. In all experiments, we produce well shaped elements: we set the facet angle bound and the radius edge bound to their theoretical limit (resp. 30 degrees and 2). We also use the same uniform sizing field for facets and cells.
We mesh an analytical sphere of radius 1.
Size bound  vertices nb  facets nb  tetrahedra nb  CPU Time (s)  vertices/second 
0.2  499  488  2,299  0.0240  20,800 
0.1  3,480  2,046  18,756  0.146  23,800 
0.05  25,556  8,274  149,703  1.50  17,000 
0.025  195,506  33,212  1,194,727  17.4  11,200 
0.0125  1,528,636  134,810  9,547,772  179  8,530 
We mesh a volume bounded by a closed triangulated surface made of about 50,000 vertices and 100,000 triangles. Picture 50.12 shows the mesh obtained when size is set to 0.005.
Size bound  vertices nb  facets nb  tetrahedra nb  CPU Time (s)  vertices/second 
0.04  423  717  1,332  0.488  866 
0.02  2,638  3,414  10,957  2.64  998 
0.01  18,168  15,576  90,338  13.9  1,310 
0.005  129,442  64,645  722,018  66.7  1,940 
0.0025  967,402  263,720  5,756,491  348  2,780 
We mesh image number 2 from the 3DIRCADb01^{1} public database. The size of this image is 512x512x172 voxels (about 45M voxels). The size of the voxels is 0.78mm x 0.78mm x 1.6mm. Picture 50.13 shows the mesh obtained for size set to 4.
Size bound (mm)  vertices nb  facets nb  tetrahedra nb  CPU Time (s)  vertices/second 
16  3,898  4,099  20,692  0.344  11,300 
8  34,117  27,792  199,864  3.09  11,000 
4  206,566  86,180  1,253,694  22.4  9,230 
2  1,546,196  329,758  9,617,278  199  7,780 
The Cgal mesh generation package implements a meshing engine based on the method of Delaunay refinement introduced by Chew [Che93] and Ruppert [Rup95] and pioneered in 3D by Shewchuk [She98b]. It uses the notion of restricted Delaunay triangulation to approximate 1dimensional curved features and curved surface patches and rely on the work of Boissonnat and Oudot [BO05] and Oudot et al. [ORY05] to achieve accurate representation of boundary and subdividing surfaces in the mesh. The mechanism of protecting balls, used to ensure a fair representation of 1dimensional features, if any, and the termination of the refinement process whatever may be the input geometry, in particular whatever small dihedral angles may form the boundary and subdivision surface patches, was pioneered by Cheng et al. [CDR07] and further experimented by Dey, Levine et al. [CDL07]. The optimization phase involves global optimization processes, a perturbation process and a sliver exudation process. The global optimizers are based on Lloyd smoothing [DFG99, DW02] and odt smoothing [Che04, ACSYD05], where odt means optimal Delaunay triangulation. The perturbation process is mainly based on the work of Tournois [Tou09] and Tournois et al. [TWAD09], while the exudation process is, the now famous, optimization by weighting described in Edelsbrunner et al. [CDE^{+}00].
From the beginning of 2009, most of the work has been performed by Stéphane Tayeb, in collaboration with Mariette Yvinec, Laurent Rineau, Pierre Alliez and Jane Tournois. First, Stéphane released the first public version of the package, implementing the specifications written by Laurent and Mariette.
The optimization processes are heavily based on the work of Jane Tournois and Pierre Alliez during the PhD of Jane advised by Pierre. The optimization phase was imported in the mesh generation package by Stéphane Tayeb and appeared first in release 3.6 of Cgal.
In collaboration with Laurent Rineau, Stéphane also added demos and examples. After some experiments on medical imaging data performed by Dobrina Boltcheva et al. [BYB09b, BYB09a], the handling of 1dimensional features was worked out by Laurent Rineau, Stéphane Tayeb and Mariette Yvinec. It appeared first in the release 3.8 of Cgal.
^{1}  available at http://www.ircad.fr/softwares/3Dircadb/3Dircadb1/index.php 