CGAL 6.1.1 - 3D Mesh Generation

Loading...
Searching...
No Matches
User Manual

Authors
  Pierre Alliez, Clément Jamin, Laurent Rineau, Stéphane Tayeb, Jane Tournois, Mariette Yvinec

Figure 63.1 Cut-view of a multi-domain 3D mesh generated from a segmented image.

# Introduction

This package is devoted to the generation of isotropic simplicial meshes discretizing 3D domains. The domain to be meshed is a subset of 3D space, required to be bounded. The domain 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 1-dimensional features (e.g. crease edges) and 0-dimensional 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 3-dimensional 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 1-dimensional 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 [9], [15], [16]. It uses the notion of restricted Delaunay triangulation to approximate 1-dimensional curves and surface patches [2]. Before the refinement, a mechanism of protecting balls is set up on 1-dimensional 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 [7], [8]. The Delaunay refinement is followed by a mesh optimization phase to remove slivers and provide a good quality mesh.

Optionally, the meshing and optimization algorithms support multi-core shared-memory architectures to take advantage of available parallelism.

## Input Domain

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 lower-dimensional 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, curves, 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. Curves may be for instance straight segments, parameterized curves or polylines. Each of those features will be accurately represented in the final mesh.

The 0 and 1-dimensional 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 1-dimensional 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 indices. If a subdomain is not connected, its different components receive the same index. Likewise different surface patches, 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.

The domain is passed 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 1-dimensional 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, 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 construct intersection points, if any. Lastly, if the input domain includes 1-dimensional 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, 1-dimensional features may be defined as segments and polyline segments.

## Output Mesh

The resulting mesh is output as a subcomplex of a weighted Delaunay 3D triangulation, in a class providing various iterators on mesh elements.

The 3D triangulation provides approximations of the subdomains, surface patches, curves, 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 1-dimensional exposed features are approximated by sequences of mesh edges and the 0-dimensional exposed features are represented by mesh vertices.

It is possible to extract the facets of the complex as a `FaceGraph`, using the function `facets_in_complex_3_to_triangle_mesh()`.

## Delaunay Refinement

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 1-dimensional 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 extent, the user may tune the Delaunay refinement to a prescribed trade-off 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 ratio between its circumradius and its shortest edge length.

## Protection of 0 and 1-dimensional Exposed Features

If the domain description includes 0 dimensional features, the corresponding points are inserted into the Delaunay triangulation from the start.

If the domain has 1-dimensional exposed features, the method of protecting balls [7], [8] 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 1-feature or the angles formed by two 1-features incident to a 0-feature.

According to this method, the 1-dimensional 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 1-features 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:

- that each segment joining two successive centers on a 1-dimensional feature will stay in the triangulation, thus ensuring an accurate approximation of the 1-dimensional features;
- 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.

## Optimization Phase

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, some sliver-shaped tetrahedra may still be present in the mesh at the end of the refinement process. The optimization phase aims to eliminate these slivers.

The optimization phase is a sequence of optimization processes, amongst the following available optimizers: an ODT-smoother, a Lloyd-smoother, a sliver perturber, and a sliver exuder.

The Lloyd and ODT-smoother are global optimizers moving the mesh vertices to minimize a mesh energy. Those optimizers are described respectively in [11], [10] and in [5], [1]. 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 ODT-smoother, 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.

Figure 63.2 The one dimensional illustration of the mesh energy minimized by ODT (left) and Lloyd (right) smoothers.

The Lloyd optimizer is known to be blind to the occurrence of slivers in the mesh while the ODT-smoother 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 the exuder applied afterwards, see Figure 63.3

The perturber and the exuder focus on improving the worst mesh elements. The perturber [18] improves the meshes by local changes in the vertices positions aiming to make slivers disappear. The exuder [6] chases the remaining slivers by re-weighting 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.

Figure 63.3 Compared effect of a global optimizer and the perturber. The left part shows the distribution of dihedral angles of mesh cells, right after Delaunay refinement (top), after some ODT smoothing (middle) and after perturbation (bottom). The numbers under the histograms give the measure in degrees of the smallest and biggest dihedral angles in the mesh.

## Post-processing

The tetrahedral mesh generation algorithm by Delaunay refinement provided by this package does not guarantee that all the vertices of the output mesh are actually present in the final mesh.

In most cases, all points are used, but if the geometry of the object has small features, compared to the size of the simplices (triangles and tetrahedra), it might be that the Delaunay facets that are selected in the restricted Delaunay triangulation miss some vertices of the triangulation. The concurrent version of the tetrahedral mesh generation algorithm also inserts a small set of auxiliary vertices that belong to the triangulation but are isolated from the complex at the end of the meshing process.

These so-called isolated vertices belong to the triangulation but not to any cell of the `C3T3`. They can be removed using the function `remove_isolated_vertices()` of `CGAL::Mesh_complex_3_in_triangulation_3`.

# Interface

As of CGAL 5.6, this package uses Named Parameters to set parameters. More details are provided in Upgrading Code using Boost Parameters to CGAL Named Function Parameters.

## The Global Functions

A 3D mesh generation process is launched through a call to one of the two following functions:

template <class C3T3, class MeshDomain, class MeshCriteria, class NamedParameters>
C3T3 make_mesh_3(const MeshDomain& domain,
const MeshCriteria& criteria,
const NamedParameters& np);
template <class C3T3, class MeshDomain, class MeshCriteria, class NamedParameters>
void refine_mesh_3(C3T3& c3t3,
const MeshDomain& domain,
const MeshCriteria& criteria,
const NamedParameters& np);
CGAL::refine_mesh_3
void refine_mesh_3(C3T3 &c3t3, const MeshDomain &domain, const MeshCriteria &criteria, const NamedParameters &np=parameters::default_values())
The function refine_mesh_3() is a 3D mesh generator.
Definition: refine_mesh_3.h:264
CGAL::make_mesh_3
C3T3 make_mesh_3(const MeshDomain &domain, const MeshCriteria &criteria, const NamedParameters &np=parameters::default_values())
The function make_mesh_3() is a 3D mesh generator.
Definition: make_mesh_3.h:607

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 1-dimensional features does not rely on Delaunay refinement, the function `refine_mesh_3()` has no parameter to preserve features.

The following sections describe the different template parameters (and their requirements) of these two global functions.

## The Data Structure

The template parameter `C3T3` is required to be a model of the concept `MeshComplex_3InTriangulation_3`, a data structure devised to represent a three dimensional complex embedded in a 3D triangulation. In both functions, an instance of type `C3T3` is used to maintain the current approximating simplicial mesh and to represent the final 3D mesh at the end of the procedure.

The embedding 3D triangulation is required to be the nested type `CGAL::Mesh_triangulation_3::type`, provided by the class template `CGAL::Mesh_triangulation_3`. The type for this triangulation is a wrapper around the class `CGAL::Regular_triangulation_3` whose vertex and cell base classes are respectively models of the concepts `MeshVertexBase_3` and `MeshCellBase_3`.

## The Domain Oracle and the Features Parameter

The template parameter `MeshDomain` is required to be a model of the concept `MeshDomain_3`. The argument `domain` of type `MeshDomain` is the sole link through which the domain to be discretized is known by the mesh generation algorithm.

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 1-dimensional features that have to be accurately represented in the final mesh, the template parameter `MeshDomain` is required to be of a model of the concept `MeshDomainWithFeatures_3`. The concept `MeshDomainWithFeatures_3` mainly provides the incidence graph of 0, 1 and 2-dimensional features, and a member function to construct sample points on curves.

Users whose domain is a model of `MeshDomainWithFeatures_3` can choose to have the corners and curves of the domain represented in the mesh or not, using the following parameters:

- `parameters::features(domain)` sets features according to the domain, i.e. 0 and 1-dimensional features are taken into account if `domain` is a `MeshDomainWithFeatures_3`.
- `parameters::no_features()` prevents the representation of 0 and 1-dimensional features in the mesh. This is useful to get a smooth and rough approximation of a domain with features.

## The Meshing Criteria

The template parameter `MeshCriteria` must be a model of the concept `MeshCriteria_3`, or a model of the refined concept `MeshCriteriaWithFeatures_3` if the domain has exposed features. The argument of type `MeshCriteria` passed to the mesh generator specifies the size and shape requirements for the tetrahedra in the mesh and for the triangles in the boundary surface mesh. These criteria condition the rules that drive the refinement process. At the end of the refinement process, mesh elements satisfy the criteria. Note that this may not be strictly true anymore after the optimization phase, but this last phase is devised to only improve the mesh quality.

The criteria for surface facets are governed by the four following parameters:

- `facet_angle`. This parameter controls the shape of surface facets. Specifically, it is a lower bound for the angle (in degrees) of surface facets. When boundary surfaces are smooth, the termination of the meshing process is guaranteed if this angular bound is at most 30 degrees [9].
- `facet_size`. This parameter controls the size of surface facets. Each surface facet has a surface Delaunay ball which is a ball circumscribing the surface facet and centered on the surface patch. The parameter `facet_size` is either a constant or a spatially variable scalar field, providing an upper bound for the radii of surface Delaunay balls.
- `facet_distance`. This parameter controls the approximation error of boundary and subdivision surfaces. Specifically, it is either a constant or a spatially variable scalar field. It provides an upper bound for the distance between the circumcenter of a surface facet and the center of a surface Delaunay ball of this facet.
- `facet_topology`. This parameter controls the set of topological constraints which have to be verified by each surface facet. By default, each vertex of a surface facet has to be located on a surface patch, on a curve, or on a corner. It can also be set to check whether the three vertices of a surface facet belongs to the same surface patch. This has to be done cautiously, as such a criterion needs that each intersection of input surface patches is an input 1-dimensional feature.

The criteria for mesh cells are governed by two parameters:

- `cell_radius_edge_ratio`. This parameter controls the shape of mesh cells (but can't filter slivers, as we discussed earlier). It is an upper bound for the ratio between the circumradius of a mesh tetrahedron and its shortest edge. There is a theoretical bound for this parameter: the Delaunay refinement process is guaranteed to terminate for values of `cell_radius_edge_ratio` bigger than 2.
- `cell_size`. This parameter controls the size of mesh tetrahedra. It is either a scalar or a spatially variable scalar field. It provides an upper bound on the circumradii of the mesh tetrahedra.

Figure 63.4 shows how the mesh generation process behaves with respect to these parameters.

Figure 63.4 Top: the mesh is obtained using the parameters \( (25,0.15,0.05)\) for the angular bound, radius bound and distance bound of surface facets and \( (4,0.2) \) for the radius-edge bound and radius bound of mesh cells. The result is a uniform mesh which contains tetrahedra of about the same size. Bottom left : the mesh is obtained by relaxing the size bound of tetrahedra and facets. The result is a small coarse mesh. Bottom middle : the mesh is obtained from the previous one by tightening the distance bound of surface facets to \( 0.01\). The result is then a graded 3D mesh with a dense surface mesh achieving a precise approximation. Bottom right : the mesh is obtained from the previous one by fixing radius bound of surface facets to \( 0.01\). The surface mesh is then denser to achieve the size bound.

If the domain has 1-dimensional exposed features, the `criteria` includes a sizing field to guide the sampling of 1-dimensional features with protecting balls centers.

- `edge_size`. This constant or variable scalar field is used as an upper bound for the distance between two protecting ball centers that are consecutive on a 1-feature. This parameter has to be set to a positive value when 1-dimensional features protection is used.

## The Optimization Parameters

The four additional parameters are optimization parameters. They control which optimization processes are performed and enable 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:

- `parameters::lloyd()` and `parameters::no_lloyd()` activate and deactivate the Lloyd smoother.
- `parameters::odt()` and `parameters::no_odt()` activate and deactivate the ODT-smoother.
- `parameters::perturb()` and `parameters::no_perturb()` activate and deactivate the perturber.
- `parameters::exude()` and `parameters::no_exude()` activate and deactivate the exuder.

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: ODT-smoother, Lloyd-smoother, perturber, and exuder.

template< class C3T3, class MeshDomain >
Mesh_optimization_return_code lloyd_optimize_mesh_3(C3T3& c3t3, const MeshDomain& domain);
template< class C3T3, class MeshDomain >
Mesh_optimization_return_code odt_optimize_mesh_3(C3T3& c3t3, const MeshDomain& domain);
template< class C3T3, class MeshDomain >
Mesh_optimization_return_code perturb_mesh_3(C3T3& c3t3, const MeshDomain& domain);
template< class C3T3 >
Mesh_optimization_return_code exude_mesh_3(C3T3& c3t3);
CGAL::Mesh_optimization_return_code
Mesh_optimization_return_code
CGAL::perturb_mesh_3
Mesh_optimization_return_code perturb_mesh_3(C3T3 &c3t3, const MeshDomain &domain, const NamedParameters &np=parameters::default_values())
The function perturb_mesh_3() is a mesh optimizer that improves the quality of a Delaunay mesh by cha...
Definition: perturb_mesh_3.h:105
CGAL::odt_optimize_mesh_3
Mesh_optimization_return_code odt_optimize_mesh_3(C3T3 &c3t3, const MeshDomain &domain, const NamedParameters &np=parameters::default_values())
The function odt_optimize_mesh_3() is a mesh optimization process based on the minimization of a glob...
Definition: odt_optimize_mesh_3.h:133
CGAL::lloyd_optimize_mesh_3
Mesh_optimization_return_code lloyd_optimize_mesh_3(C3T3 &c3t3, const MeshDomain &domain, const NamedParameters &np=parameters::default_values())
The function lloyd_optimize_mesh_3() is a mesh optimization process based on the minimization of a gl...
Definition: lloyd_optimize_mesh_3.h:138
CGAL::exude_mesh_3
Mesh_optimization_return_code exude_mesh_3(C3T3 &c3t3, const NamedParameters &np=parameters::default_values())
The function exude_mesh_3() performs a sliver exudation process on a Delaunay mesh.
Definition: exude_mesh_3.h:95

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.

## Parallel Algorithms

Enabling parallel meshing and optimization algorithms is achieved through setting the third template parameter of the `Mesh_triangulation_3` class to `Parallel_tag`, when defining the triangulation type. Note that when the user provides his/her own vertex and cell base classes, the `MeshVertexBase_3` and `MeshCellBase_3` concepts impose additional requirements.

Parallel algorithms require the executable to be linked against the Intel TBB library. To control the number of threads used, the user may use the tbb::task_scheduler_init class. See the TBB documentation for more details.

# Input/Output

Several file formats are supported for writing a mesh:

- Medit File Format, using `CGAL::IO::write_MEDIT()`
- VTK (VTU / VTP / legacy) File Formats, using `CGAL::IO::output_to_vtu()`
- Avizo File Format, using `CGAL::IO::output_to_avizo()`
- Tetgen File Format, using `CGAL::IO::output_to_tetgen()`

# Examples

## 3D Domains Bounded by Isosurfaces

### 3D Domains Bounded by Implicit Isosurfaces

The following code produces a 3D mesh for a domain whose boundary surface is an isosurface defined by an implicit function. Figure 63.5 shows a cut view of the resulting mesh.

Note the use of named parameters in the constructor of the `Mesh_criteria` instance.

File 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/Labeled_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::Labeled_mesh_domain_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
// Function
FT sphere_function (const Point& p)
{ return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-1; }
int main()
{
Mesh_domain domain =
Mesh_domain::create_implicit_mesh_domain( sphere_function,
K::Sphere_3(CGAL::ORIGIN, K::FT(2)));
// Mesh criteria
Mesh_criteria criteria(params::facet_angle(30).facet_size(0.1).facet_distance(0.025).
cell_radius_edge_ratio(2).cell_size(0.1));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return 0;
}
CGAL::Exact_predicates_inexact_constructions_kernel
CGAL::Labeled_mesh_domain_3
The class Labeled_mesh_domain_3 implements indexed domains.
Definition: Labeled_mesh_domain_3.h:354
CGAL::Mesh_complex_3_in_triangulation_3
CGAL::Mesh_criteria_3
The class gathers the refinement criteria for mesh tetrahedra and surface facets where surface facets...
Definition: Mesh_criteria_3.h:186
unspecified_type
CGAL::IO::write_MEDIT
void write_MEDIT(std::ostream &os, const T3 &t3, const NamedParameters &np=parameters::default_values())
CGAL::ORIGIN
const CGAL::Origin ORIGIN
CGAL::squared_distance
Kernel::FT squared_distance(Type1< Kernel > obj1, Type2< Kernel > obj2)
CGAL::parameters
CGAL::Parallel_tag
CGAL::Sequential_tag

Figure 63.5 Cut view of a 3D mesh produced from an implicit domain

## Meshing Multiple Domains

### Construction from a Vector of Implicit Functions

The following code produces a 3D mesh for a domain consisting of several subdomains. It uses `Implicit_multi_domain_to_labeling_function_wrapper` as model of `ImplicitFunction` which is required by `Labeled_mesh_domain_3`.

Figure 63.6 shows a view and a cut view of the resulting mesh.

File Mesh_3/mesh_implicit_domains.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_to_labeling_function_wrapper.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include "implicit_functions.h"
// IO
#include <CGAL/IO/File_medit.h>
namespace params = CGAL::parameters;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef FT_to_point_function_wrapper<double, K::Point_3> Function;
typedef CGAL::Implicit_multi_domain_to_labeling_function_wrapper<Function>
Function_wrapper;
typedef Function_wrapper::Function_vector Function_vector;
typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Mesh Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
typedef Mesh_criteria::Facet_criteria Facet_criteria;
typedef Mesh_criteria::Cell_criteria Cell_criteria;
int main()
{
// Define functions
Function f1(&torus_function);
Function f2(&sphere_function<3>);
Function_vector v;
v.push_back(f1);
v.push_back(f2);
// Domain (Warning: Sphere_3 constructor uses square radius !)
Mesh_domain domain(Function_wrapper(v), K::Sphere_3(CGAL::ORIGIN, CGAL::square(K::FT(5))),
params::relative_error_bound(1e-6));
// Set mesh criteria
Facet_criteria facet_criteria(30, 0.2, 0.02); // angle, size, approximation
Cell_criteria cell_criteria(2., 0.4); // radius-edge ratio, size
Mesh_criteria criteria(facet_criteria, cell_criteria);
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::no_exude().no_perturb());
// Perturbation (maximum cpu time: 10s, targeted dihedral angle: default)
CGAL::perturb_mesh_3(c3t3, domain, params::time_limit(10));
// Exudation
CGAL::exude_mesh_3(c3t3, params::time_limit(12));
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
return 0;
}
CGAL::Implicit_multi_domain_to_labeling_function_wrapper
The class Implicit_multi_domain_to_labeling_function_wrapper is a helping class to get a function wit...
Definition: Implicit_to_labeling_function_wrapper.h:169
CGAL::square
NT square(const NT &x)

Figure 63.6 View (left) and cut view (right) of a 3D mesh of a torus intersecting sphere.

### Construction from a Vector of Implicit Functions and a Vector of Strings

The following code produces a 3D mesh for a domain consisting of several subdomains too. Here, the set of subdomains is given by a vector of vector of signs, whereas the set was built automatically in the previous example. Each subdomain corresponds to a sign vector [s1, s2, ..., sn] where si is the sign of the function fi(p) at a point p of the subdomain.

Figure 63.7 shows a view of the resulting mesh.

File Mesh_3/mesh_implicit_domains_2.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_to_labeling_function_wrapper.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include "implicit_functions.h"
// IO
#include <CGAL/IO/File_medit.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef FT_to_point_function_wrapper<double, K::Point_3> Function;
typedef CGAL::Implicit_multi_domain_to_labeling_function_wrapper<Function>
Function_wrapper;
typedef Function_wrapper::Function_vector Function_vector;
typedef CGAL::Labeled_mesh_domain_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;
typedef Mesh_criteria::Facet_criteria Facet_criteria;
typedef Mesh_criteria::Cell_criteria Cell_criteria;
int main()
{
// Define functions
Function f1(&torus_function);
Function f2(&sphere_function<3>);
Function_vector v;
v.push_back(f1);
v.push_back(f2);
std::vector<std::string> vps;
vps.push_back("+-");
namespace params = CGAL::parameters;
Mesh_domain domain(Function_wrapper(v, vps),
K::Sphere_3(CGAL::ORIGIN, CGAL::square(K::FT(5))),
params::relative_error_bound(1e-6));
// Set mesh criteria
Facet_criteria facet_criteria(30, 0.2, 0.02); // angle, size, approximation
Cell_criteria cell_criteria(2., 0.4); // radius-edge ratio, size
Mesh_criteria criteria(facet_criteria, cell_criteria);
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::no_exude().no_perturb());
// Perturbation (maximum cpu time: 10s, targeted dihedral angle: default)
CGAL::perturb_mesh_3(c3t3, domain, params::time_limit(10));
// Exudation
CGAL::exude_mesh_3(c3t3, params::time_limit(12));
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return 0;
}

Figure 63.7 View of a 3D mesh of the relative complement of a torus in a sphere.

### Construction of a Hybrid Domain : From an Implicit and a Polyhedral Domain

The example code of Mesh_3/mesh_hybrid_mesh_domain.cpp produces a 3D mesh for a domain consisting of several subdomains too. Here, the set of subdomains is given by two subdomains of different types : a domain defined by an implicit sphere, and a domain defined by the triangulated surface of a cube.

Figure 63.8 shows the results with or without the protection of the 1D-features.

Figure 63.8 Left: the input hybrid domain.

Middle : the output mesh without any call to `add_features()`.

Right : the output mesh with added 1D-features at the intersection of subdomains.

## Polyhedral Domains

### 3D Polyhedral Domains

The following code produces a 3D mesh for a domain defined by a polyhedral surface. Figure 63.9 shows the resulting mesh.

File 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/Polyhedron_3.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/Polyhedral_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/refine_mesh_3.h>
#include <CGAL/IO/File_medit.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;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main(int argc, char*argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("meshes/elephant.off");
// Create input polyhedron
Polyhedron polyhedron;
std::ifstream input(fname);
input >> polyhedron;
if(input.fail()){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
input.close();
if (!CGAL::is_triangle_mesh(polyhedron)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(polyhedron);
// Mesh criteria (no cell_size set)
Mesh_criteria criteria(params::facet_angle(25).facet_size(0.15).facet_distance(0.008).
cell_radius_edge_ratio(3));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::no_perturb().no_exude());
// Output
std::ofstream medit_file("out_1.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
// Set tetrahedron size (keep cell_radius_edge_ratio), ignore facets
Mesh_criteria new_criteria(params::cell_radius_edge_ratio(3).cell_size(0.03));
// Mesh refinement (and make the output manifold)
CGAL::refine_mesh_3(c3t3, domain, new_criteria, params::manifold());
// Output
medit_file.open("out_2.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return EXIT_SUCCESS;
}
CGAL::Polyhedral_mesh_domain_3
The class Polyhedral_mesh_domain_3 implements a domain defined by a simplicial polyhedral surface.
Definition: Polyhedral_mesh_domain_3.h:153
CGAL::Polyhedron_3
CGAL::is_triangle_mesh
bool is_triangle_mesh(const FaceGraph &g)
CGAL::parameters::no_exude
unspecified_type no_exude()
The function parameters::no_exude() enables the user to tell the mesh generation functions make_mesh_...

Figure 63.9 View of 3D meshes produced from a polyhedral domain. (i) is a view of file out_1.mesh and (ii) is a view of file out_2.mesh. Code from subsection Mesh_3_subsection_examples_polyhedral generates these files.

### Remeshing a Polyhedral Domain with Surfaces

In the following example we have a "bounding polyhedron" which defines the meshing domain and two surfaces inside this domain. The surfaces inside the domain may be closed surfaces as well as surfaces with boundaries. In the case of a closed surface the volume delimited by this surface is also considered as inside the domain.

File Mesh_3/mesh_polyhedral_domain_with_surface_inside.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/Polyhedron_3.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Timer.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,CGAL::Sequential_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main(int argc, char*argv[])
{
std::cout.precision(17);
std::cerr.precision(17);
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("meshes/horizons.off");
std::ifstream input(fname);
const std::string fname2 = (argc>2)?argv[2]:CGAL::data_file_path("meshes/horizons-domain.off");
std::ifstream input2(fname2);
Polyhedron sm, smbounding;
input >> sm;
input2 >> smbounding;
if(input.fail()){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
CGAL::Timer t;
t.start();
// Create domain
Mesh_domain domain(sm, smbounding);
// Get sharp features
domain.detect_features();
// Mesh criteria
Mesh_criteria criteria(params::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,
params::no_perturb().no_exude());
std::cerr << t.time() << " sec." << std::endl;
// Output
dump_c3t3(c3t3, "out");
}
CGAL::Polyhedral_mesh_domain_with_features_3
The class Polyhedral_mesh_domain_with_features_3 implements a domain whose boundary is a simplicial p...
Definition: Polyhedral_mesh_domain_with_features_3.h:114
CGAL::Timer
CGAL::Timer::time
double time() const
CGAL::Timer::start
void start()

Figure 63.10 View of 3D meshes produced from a polyhedral domain with a nested surface.

### Remeshing a Polyhedral Surface

The following code creates a polyhedral domain, with only one polyhedron, and no "bounding polyhedron". The volumetric part of the domain may be empty and should not be meshed. In this case, it is recommended to use the `CGAL::parameters::surface_only()` parameter to speedup the meshing of the surface only. This enables to remesh a surface, and is equivalent to the function `make_surface_mesh()`.

File Mesh_3/remesh_polyhedral_surface.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_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;
// Polyhedron type
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
// 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_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main()
{
// Load a polyhedron
Polyhedron poly;
std::ifstream input(CGAL::data_file_path("meshes/lion-head.off"));
input >> poly;
if (!CGAL::is_triangle_mesh(poly)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// Create a polyhedral domain, with only one polyhedron,
// and no "bounding polyhedron".
// The volumetric part will be omitted by the use of `params::surface_only()`
Mesh_domain domain(poly);
// Get sharp features
domain.detect_features(); //includes detection of borders
// Mesh criteria
Mesh_criteria criteria(params::edge_size(0.025).
facet_angle(25).
facet_size(0.1).
facet_distance(0.001));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::surface_only());
// Output the facets of the c3t3 to an OFF file. The facets will not be
// oriented.
std::ofstream off_file("out.off");
c3t3.output_boundary_to_off(off_file);
return off_file.fail() ? EXIT_FAILURE : EXIT_SUCCESS;
}

Figure 63.11 View of a remeshed surface. (Left) input mesh (Right) output mesh. Code from subsection Remeshing a Polyhedral Surface generates the file.

## Domains From 3D Images

### 3D Domains Bounded by Isosurfaces in 3D Gray-Level Images

The following example produces a 3D mesh for a domain whose boundary surface is the isosurface associated to an isovalue inside the input gray-level 3D image. In the distribution you can also find the example Mesh_3/mesh_3D_gray_vtk_image.cpp which can deal with `*.nii` as well as `DICOM` files as input.

File Mesh_3/mesh_3D_gray_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_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
#include <functional>
typedef float Image_word_type;
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Labeled_mesh_domain_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;
namespace params = CGAL::parameters;
int main(int argc, char*argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("images/skull_2.9.inr");
// Load image
CGAL::Image_3 image;
if(!image.read(fname)){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
Mesh_domain domain =
Mesh_domain::create_gray_image_mesh_domain(image, params::iso_value(2.9f).value_outside(0.f));
// Mesh criteria
Mesh_criteria criteria(params::facet_angle(30).facet_size(6).facet_distance(2).
cell_radius_edge_ratio(3).cell_size(8));
// Meshing
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return 0;
}
CGAL::Image_3
The class Image_3 is a C++ wrapper around the InrImage library.
Definition: Image_3.h:11
CGAL::Image_3::read
bool read(const char *file)
Open a 3D image file.

### Domains From Segmented 3D Images

The following code produces a 3D mesh from a 3D image. The image is a segmented medical image in which each voxel is associated a label in accordance with the tissue the voxel belongs to. The domain is therefore a multi-domain where each subdomain corresponds to a specific tissue.

In the following example, the image is read from the file `liver.inr.gz` which is encoded in the format of the library Inrimage `https://www-pequan.lip6.fr/~bereziat/inrimage/`. The resulting mesh is shown in Figure 63.12.

File 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_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_mesh_domain_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main(int argc, char* argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("images/liver.inr.gz");
CGAL::Image_3 image;
if(!image.read(fname)){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image);
// Mesh criteria
Mesh_criteria criteria(params::facet_angle(30).facet_size(6).facet_distance(4).
cell_radius_edge_ratio(3).cell_size(8));
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return 0;
}

Figure 63.12 Cut view of a 3D mesh produced from a segmented liver image. Code from subsection Mesh_3_subsection_examples_3d_image generates this file.

### Domains From Segmented 3D Images, with Weights

When a segmented image is given as input, the generated mesh surface sometimes sticks too closely to the voxels surface, causing an aliasing effect. A solution to generate a smooth and accurate output surface was described by Stalling et al in [17]. It consists in generating a second input image, made of integer coefficients called weights, and use those weights to define smoother domain boundaries. The 3D image of weights can be generated using `CGAL::Mesh_3::generate_label_weights()` as shown in the following example.

File Mesh_3/mesh_3D_weighted_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_mesh_domain_3.h>
#include <CGAL/Mesh_3/generate_label_weights.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
#include <CGAL/IO/File_binary_mesh_3.h>
#include <CGAL/tags.h>
// Domain
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
// Triangulation
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain,
CGAL::Default,
CGAL::Parallel_if_available_tag>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
// Criteria
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
namespace params = CGAL::parameters;
int main(int argc, char* argv[])
{
const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("images/420.inr");// liver.inr.gz");
CGAL::Image_3 image;
if(!image.read(fname)){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
const float sigma = (std::max)(image.vx(), (std::max)(image.vy(), image.vz()));
CGAL::Image_3 img_weights =
CGAL::Mesh_3::generate_label_weights(image, sigma);
Mesh_domain domain
= Mesh_domain::create_labeled_image_mesh_domain(image,
params::weights(img_weights).
relative_error_bound(1e-6));
Mesh_criteria criteria(params::facet_angle(30).facet_size(6).facet_distance(0.5).
cell_radius_edge_ratio(3).cell_size(8));
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::no_exude(), params::no_perturb());
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
std::ofstream bin_file("out.binary.cgal", std::ios_base::binary);
CGAL::IO::save_binary_file(bin_file, c3t3);
return 0;
}
CGAL::Mesh_3::generate_label_weights
CGAL::Image_3 generate_label_weights(const CGAL::Image_3 &image, const float &sigma)
Free function that generates a CGAL::Image_3 of weights associated to each voxel of image,...
Definition: generate_label_weights.h:315
CGAL::IO::save_binary_file
bool save_binary_file(std::ostream &os, const C3T3 &c3t3)
CGAL::Default
CGAL::Mesh_triangulation_3
The class Mesh_triangulation_3 is a class template which provides the triangulation type to be used f...
Definition: Mesh_triangulation_3.h:169
CGAL::Parallel_if_available_tag

Figure 63.13 Surface of the output mesh generated with a very small `facet_distance` without the weights (left, 25563 vertices) and with the weights (right, 19936 vertices).

### Domains From 3D Images, with a Custom Initialization

The example Mesh_3/mesh_3D_image_with_image_initialization.cpp is a modification of Mesh_3/mesh_3D_image.cpp. The goal of that example is to show how the default initialization of the triangulation, using random rays, can be replaced by a new implementation. In this case, the initialization detects all connected components in the 3D segmented image, and inserts points in the triangulation for each connected component.

// Mesh generation with a custom initialization that places points
// on the surface of each connected component of the image.
CGAL::Construct_initial_points_labeled_image<C3t3, Mesh_domain>
img_pts_generator(image, domain, CGAL::parameters::verbose(true));
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::initial_points_generator(img_pts_generator));
const auto& tr = c3t3.triangulation();
CGAL::Construct_initial_points_labeled_image
Functor for generating initial points in labeled images.
Definition: Construct_initial_points_labeled_image.h:100

The parameter `CGAL::parameters::initial_points_generator` is used. It expects a functor that returns a set of points for the mesh initialization (following the concept `InitialPointsGenerator_3`). The functor `Construct_initial_points_labeled_image` is used in this example. It constructs points using the API of the mesh domain, as follows. First, the functor `construct_intersection` is created:

typename MeshDomain::Construct_intersection construct_intersection =
domain_.construct_intersection_object();

then the `Mesh_domain::Intersection` object (a `tuple` with three elements) is constructed using a call to the functor `construct_intersection`

const typename MeshDomain::Intersection intersect =
construct_intersection(test_segment);

and eventually `intersect_index` is the element #1 of `intersect`.

const Point_3& intersect_point = std::get<0>(intersect);
const Index& intersect_index = std::get<1>(intersect);
Index

The dimension of the underlying simplex is stored as the element #2 of `intersect`.

The result of the custom initialization can be seen in Figure 63.14. The generated 3D image contains a big sphere at the center, and 50 smaller spheres, generated randomly. Without the custom initialization, only the biggest component (the sphere at the center) was initialized and meshed. With the custom initialization, the initial `c3t3` object contains points on all connected components, and all spheres are meshed.

Figure 63.14 Left: the mesh without the custom initialization, only the big sphere at the center is meshed

Right: the mesh generated after the initialization of all connected components

Note that the example Mesh_3/mesh_3D_image_with_image_initialization.cpp also shows how to create a 3D image using the undocumented API of CGAL_ImageIO.

CGAL::Image_3 image = random_labeled_image(number_of_spheres);

The code of the function `random_labeled_image()` is in the header file Mesh_3/random_labeled_image.h.

The example Mesh_3/mesh_3D_image_with_custom_initialization.cpp features a custom functor that initializes the triangulation.

A struct `Custom_Initial_points_generator` that places 1D-feature points alongside a line is created.

Finally, the example Mesh_3/mesh_3D_image_with_initial_points.cpp features a point container that initializes the triangulation using the meshing parameter `parameters::initial_points()`.

## Using Variable Sizing Field

### Sizing Field as an Analytical Function

The following example shows how to use an analytical function as sizing field.

File 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/Labeled_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::Labeled_mesh_domain_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::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;
}
};
namespace params = CGAL::parameters;
// Function
FT sphere_function (const Point& p)
{ return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-1; }
int main()
{
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain
(sphere_function, K::Sphere_3(CGAL::ORIGIN, K::FT(2))
);
// Mesh criteria
Spherical_sizing_field size;
Mesh_criteria criteria(params::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, params::no_exude().no_perturb());
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return 0;
}
CGAL::abs
NT abs(const NT &x)
CGAL::sqrt
NT sqrt(const NT &x)
CGAL::parameters::no_perturb
unspecified_type no_perturb()
The function parameters::no_perturb() enables the user to tell mesh generation global functions make_...

Figure 63.15 Cut view of a 3D mesh produced from an implicit sphere with non-constant sizing field. Code from subsection Mesh_3_subsubsection_examples_sphere_variable generates this file.

### Different Sizing Field for Different Subdomains

The following example shows how to use different size for different organs in a 3D medical image.

File 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_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_mesh_domain_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::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;
namespace params = CGAL::parameters;
int main(int argc, char* argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("images/liver.inr.gz");
// Loads image
CGAL::Image_3 image;
if(!image.read(fname)){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
// Domain
Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_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(params::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");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return 0;
}
CGAL::Mesh_constant_domain_field_3
The class Mesh_constant_domain_field_3 is a model of concept MeshDomainField_3.
Definition: Mesh_constant_domain_field_3.h:51

Figure 63.16 View of a 3D mesh produced from a 3D image with different size for different organs. Code from subsection Mesh_3_subsubsection_examples_3d_image_variable generates this file.

### Lipschitz Sizing Field

The following example shows how to use another custom sizing function, that is `k`-Lipschitz. For each subdomain, the user provides the parameter `k`, a minimal size and maximal size for cells.

File Mesh_3/mesh_polyhedral_domain_with_lipschitz_sizing.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_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Mesh_3/experimental/Lipschitz_sizing_polyhedron.h>
//Kernel
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
// Domain
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
// Sizing field
typedef CGAL::Mesh_3::Lipschitz_sizing<K, Mesh_domain, Mesh_domain::AABB_tree> Lip_sizing;
namespace params = CGAL::parameters;
int main(int argc, char*argv[])
{
const std::string fname = (argc>1) ? argv[1] : CGAL::data_file_path("meshes/fandisk.off");
std::ifstream input(fname);
Polyhedron polyhedron;
input >> polyhedron;
if (input.fail()){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
if (!CGAL::is_triangle_mesh(polyhedron)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(polyhedron);
// Get sharp features
domain.detect_features();
// Create Lipschitz sizing field
Lip_sizing lip_sizing(domain, &domain.aabb_tree());
FT min_size = 0.02;
lip_sizing.add_parameters_for_subdomain(1, //subdomain id
0.3, //k
min_size,//min_size
0.5); //max_size
// Mesh criteria
Mesh_criteria criteria(params::edge_size(min_size).
facet_angle(25).
facet_size(min_size).
facet_distance(0.005).
cell_radius_edge_ratio(3).
cell_size(lip_sizing));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return EXIT_SUCCESS;
}

Figure 63.17 View of a 3D mesh produced with a Lipschitz sizing field, for a single subdomain. Code from subsection Mesh_3_subsubsection_examples_3d_image_lip_sizing generates this mesh.

## Meshing Domains with Sharp Features

### 3D Polyhedral Domain with Edges

The following example shows how to generate a mesh from a polyhedral surface. The output mesh conforms to the sharp features of the input surface.

File Mesh_3/mesh_polyhedral_domain_with_features.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_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/IO/output_to_vtu.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main(int argc, char*argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("meshes/fandisk.off");
std::ifstream input(fname);
Polyhedron polyhedron;
input >> polyhedron;
if(input.fail()){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
if (!CGAL::is_triangle_mesh(polyhedron)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(polyhedron);
// Get sharp features
domain.detect_features();
// Mesh criteria
Mesh_criteria criteria(params::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 file("out.vtu");
CGAL::IO::output_to_vtu(file, c3t3);
// Could be replaced by:
// CGAL::IO::write_MEDIT(file, c3t3);
return EXIT_SUCCESS;
}
CGAL::IO::output_to_vtu
void output_to_vtu(std::ostream &os, const C3T3 &c3t3, Mode mode=BINARY)

Figure 63.18 View of a 3D mesh with sharp features. Code from subsection Mesh_3_subsection_example_polyhedral_with_edges generates this mesh.

### Implicit Domain With 1D Features

The following example shows how to generate a mesh from an implicit domain. We add by hand the intersection of the spheres as a sharp feature.

File 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/Labeled_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::Labeled_mesh_domain_3<K> > Mesh_domain;
// Polyline
typedef std::vector<Point> Polyline_3;
typedef std::list<Polyline_3> Polylines;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = 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 =
Mesh_domain::create_implicit_mesh_domain(sphere_function,
K::Sphere_3(Point(1, 0, 0), 6.));
// Mesh criteria
Mesh_criteria criteria(params::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_3& 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,
params::no_features());
std::ofstream medit_file("out-no-protection.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
c3t3.clear();
// Mesh generation with feature preservation
c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Output
medit_file.open("out-with-protection.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return 0;
}
CGAL::Mesh_domain_with_polyline_features_3
The class Mesh_domain_with_polyline_features_3 enables the user to add some 0- and 1-dimensional feat...
Definition: Mesh_domain_with_polyline_features_3.h:535

Figure 63.19 View of a 3D mesh with sharp features generated from two intersected implicit spheres. On the left, one can see the mesh without feature preservation, and on the right the mesh with feature preservation.

### Domains from Segmented 3D Images, with 1D Features

The example Mesh_3/mesh_3D_image_with_features.cpp is a modification of Mesh_3/mesh_3D_image.cpp. That example shows how to generate a mesh from a 3D labeled image (also known as "a segmented image"), that has 2D surfaces that intersect the box corresponding to the boundary of the image. The intersection of the 2D surface with the bounding box of the image is composed of 1D curves, and must be defined as 1D-features in the domain.

The first modification is the type of the mesh domain. Instead of being `Labeled_mesh_domain_3`, it is a `Mesh_domain_with_polyline_features_3` templated by a `Labeled_mesh_domain_3`.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Labeled_mesh_domain_3<K> Image_domain;
typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> Mesh_domain;

In the main() function, the domain is created with an additional argument - a dedicated functor that computes the one-dimensional features, which are then added to the domain.

Mesh_domain domain
= Mesh_domain::create_labeled_image_mesh_domain(image,
params::features_detector = CGAL::Mesh_3::Detect_features_in_image());
CGAL::Mesh_3::Detect_features_in_image
Functor for feature detection in labeled images.
Definition: Detect_features_in_image.h:263

The `CGAL::Mesh_3::Detect_features_in_image` functor is defined in its own header file. It computes the one-dimensional features that correspond to the intersections of the bounding box of the image with the surfaces defined by the image, as well as polylines that lie at the intersection of three or more subdomains (including the outside). It then constructs a graph of these polyline features. The named constructor adds this feature graph to the domain for later feature protection. The original feature detection algorithm was described in [12], which provides a list of possible voxel configurations. The feature detection implemented in CGAL generalizes this description.

The example Mesh_3/mesh_3D_image_with_features.cpp shows how user-specified input polylines can further be added as 1D features to the mesh domain.

Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image,
params::features_detector = CGAL::Mesh_3::Detect_features_on_image_bbox(),
params::input_features = std::cref(features_inside));//use std::cref to avoid a copy
CGAL::Mesh_3::Detect_features_on_image_bbox
Functor for feature detection in labeled images.
Definition: Detect_features_on_image_bbox.h:76

In the meshing criteria, if 1D features are added to the domain, the user can define the parameter `edge_size` of the criteria class `Mesh_criteria_3`, as follows, to set up an upper bound on the length of the edges of the 1D-mesh corresponding to the 1D-features.

Mesh_criteria criteria(params::edge_size = 6.,
params::facet_angle = 30,
params::facet_size = 6,
params::facet_distance = 4,
params::cell_radius_edge_ratio = 3,
params::cell_size = 8);

The rest of the example is similar to Mesh_3/mesh_3D_image.cpp.

Figure 63.20 shows the results without or with the protection of the 1D-features.

Figure 63.20 Left: the mesh without any call to `add_features()`.

Middle: the mesh with only the 1D-features computed by `CGAL::Mesh_3::Detect_features_on_image_bbox()`.

Right: the mesh with added 1D-features in the interior of the bounding box of the image.

### Polyhedral Complex

The following example shows how to generate a mesh from a polyhedral complex that forms a bounded domain. The domain is defined by a group of polyhedral surfaces that are conformal. For each facet of the input polyhedral surfaces, the ids of incident subdomains are known. The input surfaces intersect along polylines that are considered as sharp features, and protected as for all domain types in this Section.

File Mesh_3/mesh_polyhedral_complex.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_complex_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <cstdlib>
#include <cassert>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
typedef CGAL::Polyhedral_complex_mesh_domain_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
const char* const filenames[] = {
"meshes/patch-01.off",
"meshes/patch-13.off",
"meshes/patch-20.off",
"meshes/patch-21.off",
"meshes/patch-23.off",
"meshes/patch-30.off"
};
const std::pair<int, int> incident_subdomains[] = {
std::make_pair(0, 1),
std::make_pair(1, 3),
std::make_pair(2, 0),
std::make_pair(2, 1),
std::make_pair(2, 3),
std::make_pair(3, 0),
};
int main()
{
const std::size_t nb_patches = sizeof(filenames) / sizeof(const char*);
assert(sizeof(incident_subdomains) == nb_patches * sizeof(std::pair<int, int>));
std::vector<Polyhedron> patches(nb_patches);
for(std::size_t i = 0; i < nb_patches; ++i) {
std::ifstream input(CGAL::data_file_path(filenames[i]));
if(!(input >> patches[i])) {
std::cerr << "Error reading " << CGAL::data_file_path(filenames[i]) << " as a polyhedron!\n";
return EXIT_FAILURE;
}
}
// Create domain
Mesh_domain domain(patches.begin(), patches.end(),
incident_subdomains, incident_subdomains+nb_patches);
domain.detect_features(); //includes detection of borders
// Mesh criteria
Mesh_criteria criteria(params::edge_size(8).
facet_angle(25).facet_size(8).facet_distance(0.2).
cell_radius_edge_ratio(3).cell_size(10));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
return EXIT_SUCCESS;
}
CGAL::Polyhedral_complex_mesh_domain_3
The class Polyhedral_complex_mesh_domain_3 implements a domain defined by a collection of polyhedral ...
Definition: Polyhedral_complex_mesh_domain_3.h:120

Figure 63.21 shows the results without or with the protection of the 1D-features.

Figure 63.21 Left: the mesh without any call to `detect_features()`. The mesh is overrefined next to sharp features, because the mesher tries to reach the target `facet_distance`, by inserting centers of surface Delaunay balls, that are not likely to be on sharp edges.

Right: the mesh with added 1D-features along sharp edges and edges shared by at least two of the input polyhedral surfaces forming the complex.

### Sizing Field for Feature Edges

The following example shows how to generate a mesh from a polyhedral complex or a polyhedral surface with polyline features. Polyline features are covered by a set of so-called "protecting balls" which sizes are highly related to the edge length, and are driven by the size component of `Mesh_edge_criteria_3` (see Section Protection of 0 and 1-dimensional Exposed Features). The ideal size can be computed using `Sizing_field_with_aabb_tree` that helps start the feature protection and one-dimensional meshing process with a good initial guess. To fit the protecting balls requirements, no protecting ball can have its radius larger than half of the distance from the corresponding vertex (its center), to surface patches the current polyline feature does not belong to.

File Mesh_3/mesh_polyhedral_domain_with_features_sizing.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/Surface_mesh.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Sizing_field_with_aabb_tree.h>
// Domain
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Polyhedron = CGAL::Surface_mesh<K::Point_3>;
using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K, Polyhedron>;
using Features_sizing_field = CGAL::Sizing_field_with_aabb_tree<K, Mesh_domain>;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main(int argc, char*argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("meshes/fandisk.off");
std::ifstream input(fname);
Polyhedron polyhedron;
input >> polyhedron;
if(input.fail()){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(polyhedron);
// Get sharp features
domain.detect_features();
// Mesh criteria
Features_sizing_field edges_sizing_field(0.07, domain);
Mesh_criteria criteria(params::edge_size(edges_sizing_field).
facet_distance(0.0072).
cell_radius_edge_ratio(3));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::no_exude().no_perturb());
// Output
CGAL::dump_c3t3(c3t3, "out_sizing_field_with_aabb_tree");
return EXIT_SUCCESS;
}
CGAL::Sizing_field_with_aabb_tree
The class Sizing_field_with_aabb_tree is a model of concept MeshDomainField_3.
Definition: Sizing_field_with_aabb_tree.h:73

### Approximation Criterion for Feature Edges

The following example shows how to generate a mesh with an approximation error criterion. Polyline features are covered by a set of so-called "protecting balls" which sizes are highly related to the `edge_size` and `edge_distance` criteria. This edge distance is driven by the distance component of `Mesh_edge_criteria_3` (see Section Protection of 0 and 1-dimensional Exposed Features). The parameter `edge_distance` enables the subdivision of edges that are too far away from the corresponding input sharp features. It may result in an output mesh with shorter edges along the input polyline features.
File Mesh_3/mesh_polyhedral_domain_with_edge_distance.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_with_features_3.h>
#include <CGAL/make_mesh_3.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main(int argc, char*argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("meshes/u_arch.off");
std::ifstream input(fname);
Polyhedron polyhedron;
input >> polyhedron;
if(input.fail()){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
if (!CGAL::is_triangle_mesh(polyhedron)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(polyhedron);
// Get sharp features
domain.detect_features();
// Mesh criteria
Mesh_criteria criteria(params::edge_size(4).edge_distance(0.01).
facet_angle(25).facet_size(4).facet_distance(0.01).
cell_radius_edge_ratio(3).cell_size(4));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Output
//CGAL::dump_c3t3(c3t3, "out");
return EXIT_SUCCESS;
}

Figure 63.22 Left: the mesh without the `edge_distance` criterion.

Right: the mesh with the `edge_distance` criterion set to the same value as `facet_distance`.

## Tuning Mesh Optimization

In the previous examples, the mesh generation is launched through a call `make_mesh_3()` 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.

### Disabling Exudation and Tuning Perturbation

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 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_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_mesh_domain_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Mesh Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main(int argc, char* argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("images/liver.inr.gz");
// Domain
CGAL::Image_3 image;
if(!image.read(fname)){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image);
// Mesh criteria
Mesh_criteria criteria(params::facet_angle(30).facet_size(5).facet_distance(1.5).
cell_radius_edge_ratio(2).cell_size(7));
// Mesh generation and optimization in one call (sliver_bound is the
// targeted dihedral angle in degrees)
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::no_exude().
perturb(params::sliver_bound(10).time_limit(15)));
// Mesh generation and optimization in several call
C3t3 c3t3_bis = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::no_perturb().no_exude());
CGAL::perturb_mesh_3(c3t3_bis, domain, params::time_limit(15));
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
std::ofstream medit_file_bis("out_bis.mesh");
CGAL::IO::write_MEDIT(medit_file_bis, c3t3_bis);
medit_file_bis.close();
return 0;
}

### Using Lloyd Global Optimization

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 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_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_mesh_domain_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Mesh Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main(int argc, char*argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("images/liver.inr.gz");
// Domain
CGAL::Image_3 image;
if(!image.read(fname)){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image);
// Mesh criteria
Mesh_criteria criteria(params::facet_angle(30).facet_distance(1.2).
cell_radius_edge_ratio(2));
// Mesh generation and optimization in one call
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::lloyd(params::time_limit(30)).
no_perturb().
exude(params::time_limit(10).sliver_bound(10)));
// Mesh generation and optimization in several call
C3t3 c3t3_bis = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::no_perturb().no_exude());
CGAL::lloyd_optimize_mesh_3(c3t3_bis, domain, params::time_limit(30));
CGAL::exude_mesh_3(c3t3_bis, params::sliver_bound(10), params::time_limit(10));
// Output
std::ofstream medit_file("out.mesh");
CGAL::IO::write_MEDIT(medit_file, c3t3);
medit_file.close();
std::ofstream medit_file_bis("out_bis.mesh");
CGAL::IO::write_MEDIT(medit_file_bis, c3t3_bis);
medit_file_bis.close();
return 0;
}

# Performance

We provide here some benchmarks of the performance of the mesh generation algorithms.

## Delaunay Refinement

The computer used for benchmarking 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. These benchmarks have been done using CGAL v3.8. Note that these benchmarks were obtained with the sequential version of the algorithm, which does not take advantage of multi-core architectures. See the next section for performance of parallel algorithms.

We study the refinement part of the mesh generation engine in this section. We give the CPU time (measured by `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.

### Implicit Function

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

### Polyhedral Domain

Figure 63.23 View of polyhedral mesh generation result (size = 0.005).

We mesh a volume bounded by a closed triangulated surface made of about 50,000 vertices and 100,000 triangles. Figure 63.23 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

### 3D Image

Figure 63.24 View of 3d image mesh generation result (size = 4).

We mesh image number 2 from the 3D-IRCADb-01 (available at https://www.ircad.fr/research/data-sets/liver-segmentation-3d-ircadb-01/) 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. Figure 63.24 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

## Parallel Performance

We provide below speed-up charts generated using the parallel version of the meshing algorithms of CGAL 4.5. The machine used is a PC running Windows 7 64-bits with two 6-core Intel Xeon CPU X5660 clocked at 2.80 GHz with 32GB of RAM. The program has been compiled with Microsoft Visual C++ 2012 in Release mode.

Figure 63.25 shows mesh refinement speed-up, and Figure 63.26 shows Lloyd optimization speed-up. ODT optimization exhibits similar speed-up.

Figure 63.25 Facet refinement speed-up (left) and cell refinement speed-up (right), compared to the sequential version of the algorithm.

Figure 63.26 Lloyd optimization speed-up, compared to the sequential version of the algorithm. Parallel ODT exhibits similar performance.

# Design and Implementation History

## Theoretical Foundations

The CGAL mesh generation package implements a meshing engine based on the method of Delaunay refinement introduced by Chew [9] and Ruppert [15] and pioneered in 3D by Shewchuk [16]. It uses the notion of restricted Delaunay triangulation to approximate 1-dimensional curved features and curved surface patches and rely on the work of Boissonnat and Oudot [2] and Oudot et al. [13] to achieve accurate representation of boundary and subdividing surfaces in the mesh. The mechanism of protecting balls, used to ensure a fair representation of 1-dimensional 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. [8] and further experimented by Dey, Levine et al. [7]. The optimization phase involves global optimization processes, a perturbation process and a sliver exudation process. The global optimizers are based on Lloyd smoothing [11], [10] and odt smoothing [5], [1], where odt means optimal Delaunay triangulation. The perturbation process is mainly based on the work of Tournois [20] and Tournois et al. [19], while the exudation process is, the now famous, optimization by weighting described in Edelsbrunner et al. [6].

## Implementation History

Work on the package `Mesh_3` started during the PhD thesis of Laurent Rineau advised by Mariette Yvinec. A code prototype, together with a first version of design and specifications [14] came out of their collaboration.

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. [4], [3], the handling of 1-dimensional features was worked out by Laurent Rineau, Stéphane Tayeb and Mariette Yvinec. It appeared first in the release 3.8 of CGAL.

In 2013, Clément Jamin made the meshing and optimization algorithms parallel on multi-core shared-memory architectures.