CGAL 6.1 - Polygon Mesh Processing

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

Authors
  David Coeurjolly, Jaques-Olivier Lachaud, Konstantinos Katrioplas, Sébastien Loriot, Ivan Pađen, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, Sébastien Valette, Léo Valque, and Ilker O. Yaz

# Introduction

This package implements a collection of methods and classes for polygon mesh processing, ranging from basic operations on simplices, to complex geometry processing algorithms. The implementation of this package mainly follows algorithms and references given in Botsch et al.'s book on polygon mesh processing [6].

## Polygon Mesh

A polygon mesh is a consistent and orientable surface mesh, that can have one or more boundaries. The faces are simple polygons. The edges are segments. Each edge connects two vertices, and is shared by two faces (including the null face for boundary edges). A polygon mesh can have any number of connected components, and also some self-intersections. In this package, a polygon mesh is considered to have the topology of a 2-manifold.

## API

This package follows the BGL API described in CGAL and the Boost Graph Library. It can thus be used either with `Polyhedron_3`, `Surface_mesh`, or any class model of the concept `FaceGraph`. Each function or class of this package details the requirements on the input polygon mesh.

Named Parameters are used to deal with optional parameters. The page Named Parameters describes their usage.

## Outline

The algorithms described in this manual are organized in sections:

- Meshing : meshing algorithms, including triangulation of non-triangulated meshes, refinement, optimization by fairing, remeshing of triangulated surface meshes and smoothing algorithms.
- Corefinement and Boolean Operations : methods to corefine triangle meshes and to compute boolean operations out of corefined closed triangle meshes.
- Hole Filling : available hole filling algorithms, which can possibly be combined with refinement and fairing.
- Predicates : predicates that can be evaluated on the processed polygon. mesh, which includes point location and self intersection tests.
- Orientation : checking or fixing the orientation of a polygon soup.
- Combinatorial Repair : repair of polygon meshes and polygon soups.
- Geometric Repair : repair of the geometry of polygon meshes.
- Computing Normals : normal computation at vertices and on faces of a polygon mesh.
- Computing Curvatures : computing curvatures (mean, gaussian, principal) on a polygon mesh.
- Slicer : functor able to compute the intersections of a polygon mesh with arbitrary planes (slicer).
- Connected Components : methods to deal with connected components of a polygon mesh (extraction, marks, removal, ...).

## Reading and Writing Polygon Meshes

In all functions of this package, the polygon meshes are required to be models of the graph concepts defined in the package CGAL and the Boost Graph Library. Using common graph concepts enables having common input/output functions for all the models of these concepts. The page I/O Functions provides an exhaustive description of the available I/O functions. In addition, this package offers the function `CGAL::Polygon_mesh_processing::IO::read_polygon_mesh()`, which can perform some repairing if the input data do not represent a manifold surface.

# Meshing

A surface patch can be refined by inserting new vertices and flipping edges to get a triangulation. Using a criterion presented in [14], the density of triangles near the boundary of the patch is approximated by the refinement function. The validity of the mesh is enforced by flipping edges. An edge is flipped only if the opposite edge does not exist in the original mesh and if no degenerate triangles are generated.

A region of the surface mesh (e.g. the refined region), can be faired to obtain a tangentially continuous and smooth surface patch. The region to be faired is defined as a range of vertices that are relocated. The fairing step minimizes a linear bi-Laplacian system with boundary constraints, described in [5]. The visual results of aforementioned steps are depicted by Figure 75.14 (c and d).

## API

### Meshing

Refinement and fairing functions can be applied to an arbitrary region on a triangle mesh, using :

- `CGAL::Polygon_mesh_processing::refine()` : given a set of facets on a mesh, refines the region.
- `CGAL::Polygon_mesh_processing::fair()` : given a set of vertices on a mesh, fairs the region.

Fairing needs a sparse linear solver and we recommend the use of Eigen 3.2 or later. Note that fairing might fail if fixed vertices, which are used as boundary conditions, do not suffice to solve the constructed linear system.

Many algorithms require as input meshes in which all the faces have the same degree, or even are triangles. Hence, one may want to triangulate all polygon faces of a mesh.

This package provides the function `CGAL::Polygon_mesh_processing::triangulate_faces()` that triangulates all faces of the input polygon mesh. An approximated support plane is chosen for each face, orthogonal to the normal vector computed by `CGAL::Polygon_mesh_processing::compute_face_normal()`. Then, the triangulation of each face is the one obtained by building a `CGAL::Constrained_Delaunay_triangulation_2` in this plane. This choice is made because the constrained Delaunay triangulation is the triangulation that, given the edges of the face to be triangulated, maximizes the minimum angle of all the angles of the triangles in the triangulation.

### Remeshing

#### Local Isotropic Incremental Remeshing

The incremental triangle-based isotropic remeshing algorithm introduced by Botsch et al [4], [6] is implemented in this package. This algorithm incrementally performs simple operations such as edge splits, edge collapses, edge flips, and Laplacian smoothing. All the vertices of the remeshed patch are reprojected to the original surface to keep a good approximation of the input.

A triangulated region of a polygon mesh can be remeshed using the function `CGAL::Polygon_mesh_processing::isotropic_remeshing()`, as illustrated by Figure 75.1. The algorithm has two parameters: the sizing field object for the remeshed surface patch, and the number of iterations of the abovementioned sequence of operations.

The sizing field establishes the local target edge length for the remeshed surface. Two sizing fields are provided: a uniform and a curvature-adaptive sizing field. With `CGAL::Polygon_mesh_processing::Uniform_sizing_field`, all triangle edges are targeted to have equal lengths. With `CGAL::Polygon_mesh_processing::Adaptive_sizing_field`, triangle edge lengths depend on the local curvature – shorter edges appear in regions with a higher curvature and vice versa. The outline of the adaptive sizing field algorithm is available in [9]. The distinction between uniform and adaptive sizing fields is depicted in Figure 75.2.

As the number of iterations increases, the mesh tends to be smoother and closer to the target edge length.

An additional option has been added to protect (i.e. not modify) some given polylines. In some cases, those polylines are too long, and reaching the desired target edge length while protecting them is not possible and leads to an infinite loop of edge splits in the incident faces. To avoid that pitfall, the function `CGAL::Polygon_mesh_processing::split_long_edges()` should be called on the list of constrained edges before remeshing.

Figure 75.1 Isotropic remeshing. (a) Triangulated input surface mesh. (b) Surface uniformly and entirely remeshed. (c) Selection of a range of faces to be remeshed. (d) Surface mesh with the selection uniformly remeshed.

Figure 75.2 Sizing fields in isotropic remeshing. (a) Uniform sizing field. (b) Curvature-based adaptive sizing field.

#### Delaunay-Based Surface Remeshing

The mesh generation algorithm implemented in the 3D Mesh Generation package can be used to remesh a given triangulated surface mesh. The algorithm, based on Delaunay refinement of a restricted Delaunay triangulation, generates a triangle surface mesh that provably has the required properties on simplices size, surface approximation, shape of facets, and topology of the surface. A set of edges from the input mesh can be given as constraints to build and protect a set of polylines features, that are resampled while preserving the topology of the input feature graph.

This solution for triangle surface remeshing is available in this package with the function `CGAL::Polygon_mesh_processing::surface_Delaunay_remeshing()`. All the meshing criteria defined for `CGAL::make_mesh_3()` can be used directly for size of simplices, shape of facets, surface approximation, topology, and one-dimensional features sampling.

An example of how to remesh a given triangulated surface mesh with the Delaunay refinement algorithm while preserving the detected sharp edges is given in Polygon_mesh_processing/delaunay_remeshing_example.cpp

#### Approximated Discrete Centroidal Voronoi Diagram (ACVD) Remeshing

This remeshing algorithm uses clustering on polygonal meshes as to approximate a Centroidal Voronoi Diagram construction. It is inspired by the method presented in [18] and further developed in [3] and [19]. The algorithm is similar to Lloyd's algorithm (or k-means) where random input vertices are picked to initialize clusters of vertices, which are then grown to minimize a particular energy, until convergence. Upon convergence, output vertices are computed from the vertices of each cluster. The function `CGAL::Polygon_mesh_processing::approximated_centroidal_Voronoi_diagram_remeshing()` takes as input a triangle mesh and an expected number of output vertices, and updates the mesh with the remeshed version. Note that the final vertex count may exceed the input parameter if the input mesh is not closed, or if the specified vertex budget was insufficient to generate a manifold mesh. In this case, the uniformity of the output triangles may be affected. Note that there is no guarantee that the output mesh will have the same topology as the input. For example, if a small handle is entirely contained in a cluster then this handle will not be part of the output. If the input mesh contains sharp features or corners, it is possible to use quadric error metrics to either move output vertices (fast but can produce bad looking triangles), of to use quadric error metrics directly into the energy formulation of each cluster (slower but produces better quality triangles). Additionally, adaptive remeshing based on surface curvature is possible.

Figure 75.3 ACVD remeshing. From left to right: Input; with default parameters; using adaptative option; using QEM based energy.

#### Remeshing of (Almost) Planar Patches

When many triangles are used to describe a planar region of a model, one might wish to simplify the mesh in this region to use few elements, or even a single large polygonal face when the region makes up a simply connected patch. This can be achieved using the function `CGAL::Polygon_mesh_processing::remesh_planar_patches()`. This function performs the detection of the planar regions, using geometric predicates for coplanarity and collinearity checks. If these tests are performed exactly, the planar regions can be unexpectedly small due to the input in fact not being perfectly planar. To palliate this, it is possible to specify a threshold on the angle between adjacent faces (resp. segments) such that they are considered coplanar (resp. collinear). However, this tolerance threshold is only local and there is no global control, which can have undesired effects, such as in the classic example of a densely sampled circle arc where all points are eventually found to be almost collinear). To circumvent this situation, we provide the function `CGAL::Polygon_mesh_processing::remesh_almost_planar_patches()` , which expects the segmentation into planar patches and corners to be provided by the user. Such segmentation can be obtained using the function `CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces()`, which uses the region growing algorithm to detect planar regions in a mesh with global and local criteria. Similarly, the function `CGAL::Polygon_mesh_processing::detect_corners_of_regions()` can be used to detect corner vertices on the border of the planar regions detected by running the region growing algorithm on border segments of the patch.

(a) | (b)

Figure 75.4 Remeshing of planar patches in two models: (a) planar patches in the cheese model are retriangulated (b) the remeshed version of this model contains 14 vertices and has not been retriangulated, the patch ids have been assigned to the input and the output meshes, allowing an identical color scheme.

Figure 75.5 Remeshing of planar patches using region growing for the detection of planar patches. From left to right: input mesh, remeshed version, remeshed version using the same angle threshold but larger approximation error.

### Smoothing

Smoothing of a triangulated mesh region can be achieved with algorithms that aim at either mesh smoothing or shape smoothing. While mesh smoothing is achieved by improving the quality of triangles based on criteria such as angle and area, shape smoothing is designed to be intrinsic, depending as little as possible on the discretization and smoothing the shape alone without optimizing the shape of the triangles.

- Mesh smoothing by angle and area optimization: `CGAL::Polygon_mesh_processing::angle_and_area_smoothing()` moves vertices to optimize geometry around each vertex: it can try to equalize the angles between incident edges, or (and) move vertices so that areas of adjacent triangles tend to equalize. Border vertices are considered constrained and do not move at any step of the procedure. No vertices are inserted at any time. Angle and area smoothing algorithms are based on Surazhsky and Gotsman [16]. Since area smoothing considers only areas as a smoothing criterion, it may result in long and skinny triangles. To palliate this phenomenon, area smoothing is followed by an (optional) step of Delaunay-based edge flips. In any case, area smoothing is guaranteed to improve the spatial distribution of the vertices over the area that is being smoothed. A simple example can be found in Polygon_mesh_processing/mesh_smoothing_example.cpp.

Figure 75.6 Mesh smoothing of the closed surface blobby, containing self-intersections (circled in red). For each smoothing combination, 10 iterations were applied. From left to right: (a) Input mesh; (b) Smoothing based on areas without using Delaunay flips; (c) Smoothing based on areas with Delaunay flips; (d) Smoothing based on angles; (e) Smoothing based on angles and areas, with Delaunay flips.

Figure 75.7 Statistics for the various combinations of mesh smoothing.

- Mesh smoothing by tangential relaxation: `CGAL::Polygon_mesh_processing::tangential_relaxation()` moves vertices following an area-based Laplacian smoothing scheme, performed at each vertex in an estimated tangent plane to the surface. The full algorithm is described in [6]. The example Polygon_mesh_processing/tangential_relaxation_example.cpp shows how this mesh relaxation function can be used.
- Shape smoothing: `CGAL::Polygon_mesh_processing::smooth_shape()` moves vertices towards a weighted barycenter of their neighbors along the mean curvature flow. The curvature flow algorithm for shape smoothing is based on Desbrun et al. [8] and on Kazhdan et al. [10]. The algorithm uses the mean curvature flow to calculate the translation of vertices along the surface normal with a speed equal to the mean curvature of the area that is being smoothed. This means that vertices on sharp corners slide faster. If the region around a vertex is flat, this vertex does not move (zero curvature). To avoid the formation of undesirable neck pinches (cylindrical surface areas that form singularities) the algorithm slows down the evolution in cylindrical regions. The smoothed shape converges to a sphere while staying conformally equivalent to its original shape. A simple example can be found in Polygon_mesh_processing/shape_smoothing_example.cpp.

Figure 75.8 Shape smoothing of the devil model, using the mean curvature flow with a time step equal to 0.05 and constraining border vertices (located at the neck, where the mesh is open).

### Extrusion

This package provides two functions to extrude a triangle mesh with boundaries. A first one extrudes by offsetting each vertex by a user provided vector. A second one is more general: It takes two users provided functors which extrude "up" and "down". In both cases the boundaries are connected by a triangle strip. Note that the extrusion may generate self intersecting surfaces.

## Meshing Examples

### Refine and Fair a Region on a Triangle Mesh

The following example calls the functions `CGAL::Polygon_mesh_processing::refine()` and `CGAL::Polygon_mesh_processing::fair()` for some selected regions on the input triangle mesh.

File Polygon_mesh_processing/refine_fair_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/refine.h>
#include <CGAL/Polygon_mesh_processing/fair.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <iostream>
#include <iterator>
#include <map>
#include <utility>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Polyhedron_3<Kernel> Mesh;
typedef Mesh::Vertex_handle Vertex_handle;
namespace PMP = CGAL::Polygon_mesh_processing;
// extract vertices which are at most k (inclusive)
// far from vertex v in the graph of edges
void extract_k_ring(Vertex_handle v,
int k,
std::vector<Vertex_handle>& qv)
{
std::map<Vertex_handle, int> D;
qv.push_back(v);
D[v] = 0;
std::size_t current_index = 0;
int dist_v;
while (current_index < qv.size() && (dist_v = D[qv[current_index]]) < k)
{
v = qv[current_index++];
Mesh::Halfedge_around_vertex_circulator e(v->vertex_begin()), e_end(e);
do {
Vertex_handle new_v = e->opposite()->vertex();
if (D.insert(std::make_pair(new_v, dist_v + 1)).second)
qv.push_back(new_v);
} while (++e != e_end);
}
}
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off");
Mesh poly;
if(!PMP::IO::read_polygon_mesh(filename, poly) || !CGAL::is_triangle_mesh(poly))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
std::vector<Mesh::Facet_handle> new_facets;
std::vector<Vertex_handle> new_vertices;
PMP::refine(poly, faces(poly),
std::back_inserter(new_facets),
std::back_inserter(new_vertices),
CGAL::parameters::density_control_factor(2.));
std::ofstream refined_off("refined.off");
refined_off.precision(17);
refined_off << poly;
refined_off.close();
std::cout << "Refinement added " << new_vertices.size() << " vertices." << std::endl;
Mesh::Vertex_iterator v = poly.vertices_begin();
std::advance(v, 82/*e.g.*/);
std::vector<Vertex_handle> region;
extract_k_ring(v, 12/*e.g.*/, region);
bool success = PMP::fair(poly, region);
std::cout << "Fairing : " << (success ? "succeeded" : "failed") << std::endl;
std::ofstream faired_off("faired.off");
faired_off.precision(17);
faired_off << poly;
faired_off.close();
return 0;
}
CGAL::Exact_predicates_inexact_constructions_kernel
CGAL::Polyhedron_3
CGAL::is_triangle_mesh
bool is_triangle_mesh(const FaceGraph &g)
Kernel

### Triangulate a Polygon Mesh

Triangulating a polygon mesh can be achieved through the function `CGAL::Polygon_mesh_processing::triangulate_faces()` as shown in the following example.

File Polygon_mesh_processing/triangulate_faces_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/boost/graph/helpers.h>
#include <iostream>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/P.off");
const char* outfilename = (argc > 2) ? argv[2] : "P_tri.off";
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Error: Invalid input." << std::endl;
return EXIT_FAILURE;
}
if(is_empty(mesh))
{
std::cerr << "Warning: empty file?" << std::endl;
return EXIT_SUCCESS;
}
if(!CGAL::is_triangle_mesh(mesh))
std::cout << "Input mesh is not triangulated." << std::endl;
else
std::cout << "Input mesh is triangulated." << std::endl;
PMP::triangulate_faces(mesh);
// Confirm that all faces are triangles.
for(boost::graph_traits<Mesh>::face_descriptor f : faces(mesh))
{
if(!CGAL::is_triangle(halfedge(f, mesh), mesh))
std::cerr << "Error: non-triangular face left in mesh." << std::endl;
}
CGAL::IO::write_polygon_mesh(outfilename, mesh, CGAL::parameters::stream_precision(17));
return EXIT_SUCCESS;
}
CGAL::Surface_mesh
Kernel::Point_3
CGAL::is_triangle
bool is_triangle(typename boost::graph_traits< FaceGraph >::halfedge_descriptor hd, const FaceGraph &g)
CGAL::is_empty
bool is_empty(const FaceGraph &g)
CGAL::IO::write_polygon_mesh
bool write_polygon_mesh(const std::string &fname, Graph &g, const NamedParameters &np=parameters::default_values())

An additional parameter, named `visitor` can be used to track the how faces are triangulated into subfaces. The following examples shows how to use it.

File Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/boost/graph/copy_face_graph.h>
#include <iostream>
#include <fstream>
#include <string>
#include <unordered_map>
#include <utility>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Mesh;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
class Insert_iterator
{
typedef std::unordered_map<face_descriptor,face_descriptor> Container;
Container& container;
public:
Insert_iterator(Container &c)
: container(c) {}
Insert_iterator&
operator=(const std::pair<face_descriptor, face_descriptor>& p)
{
container[p.second] = p.first;
return *this;
}
Insert_iterator&
operator*() { return *this; }
Insert_iterator
operator++(int) { return *this; }
};
struct Visitor : public CGAL::Polygon_mesh_processing::Triangulate_faces::Default_visitor<Mesh>
{
typedef std::unordered_map<face_descriptor,face_descriptor> Container;
Container& container;
face_descriptor qfd;
Visitor(Container& container)
: container(container)
{}
void before_subface_creations(face_descriptor fd)
{
std::cout << "split : " << fd << " into:" << std::endl;
Container::iterator it = container.find(fd);
qfd = it->second;
container.erase(it);
}
void after_subface_created(face_descriptor fd)
{
std::cout << " " << fd;
container[fd]=qfd;
}
void after_subface_creations()
{
std::cout << std::endl;
}
};
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/P.off");
std::ifstream input(filename);
Mesh mesh;
if (!input || !(input >> mesh) || mesh.is_empty())
{
std::cerr << "Not a valid off file." << std::endl;
return 1;
}
std::unordered_map<face_descriptor,face_descriptor> t2q;
Mesh copy;
CGAL::copy_face_graph(mesh, copy, CGAL::parameters::face_to_face_output_iterator(Insert_iterator(t2q)));
Visitor v(t2q);
CGAL::Polygon_mesh_processing::triangulate_faces(copy,
CGAL::parameters::visitor(v));
for(std::unordered_map<face_descriptor,face_descriptor>::iterator it = t2q.begin(); it != t2q.end(); ++it){
std::cout << it->first << " " << it->second << std::endl;
}
return 0;
}
CGAL::Insert_iterator::Insert_iterator
Insert_iterator(Container &c)
Container
operator*
Vector_2< Kernel > operator*(const Vector_2< Kernel > &v, const Kernel::RT &s)
CGAL::Polygon_mesh_processing::triangulate_faces
bool triangulate_faces(FaceRange face_range, PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
triangulates given faces of a polygon mesh.
Definition: triangulate_faces.h:373
CGAL::copy_face_graph
void copy_face_graph(const SourceMesh &sm, TargetMesh &tm, const NamedParameters1 &np1=parameters::default_values(), const NamedParameters2 &np2=parameters::default_values())
CGAL::Polygon_mesh_processing::Triangulate_faces::Default_visitor
Default new face visitor model of PMPTriangulateFaceVisitor.
Definition: triangulate_faces.h:48

### Isotropic Remeshing of a Region on a Polygon Mesh

The following example shows a complete example of how the isotropic remeshing function can be used. First, the border of the polygon mesh is collected. Since the boundary edges will be considered as constrained and protected in this example, the function `split_long_edges()` is called first on these edges.

Once this is done, remeshing is run on all the surface, with protection of constraints activated, for 3 iterations.

File Polygon_mesh_processing/isotropic_remeshing_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <boost/iterator/function_output_iterator.hpp>
#include <iostream>
#include <string>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
typedef boost::graph_traits<Mesh>::edge_descriptor edge_descriptor;
namespace PMP = CGAL::Polygon_mesh_processing;
struct halfedge2edge
{
halfedge2edge(const Mesh& m, std::vector<edge_descriptor>& edges)
: m_mesh(m), m_edges(edges)
{}
void operator()(const halfedge_descriptor& h) const
{
m_edges.push_back(edge(h, m_mesh));
}
const Mesh& m_mesh;
std::vector<edge_descriptor>& m_edges;
};
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/pig.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
double target_edge_length = (argc > 2) ? std::stod(std::string(argv[2])) : 0.04;
unsigned int nb_iter = (argc > 3) ? std::stoi(std::string(argv[3])) : 10;
std::cout << "Split border...";
std::vector<edge_descriptor> border;
PMP::border_halfedges(faces(mesh), mesh, boost::make_function_output_iterator(halfedge2edge(mesh, border)));
PMP::split_long_edges(border, target_edge_length, mesh);
std::cout << "done." << std::endl;
std::cout << "Start remeshing of " << filename
<< " (" << num_faces(mesh) << " faces)..." << std::endl;
PMP::isotropic_remeshing(faces(mesh), target_edge_length, mesh,
CGAL::parameters::number_of_iterations(nb_iter)
.protect_constraints(true)); //i.e. protect border, here
CGAL::IO::write_polygon_mesh("out.off", mesh, CGAL::parameters::stream_precision(17));
std::cout << "Remeshing done." << std::endl;
return 0;
}

### Extrusion of a Triangle Mesh

The following example extrudes in both directions by projecting points along the vertex normal. In case the function `Bottom` was the identity this would have been a one-sided extrusion.

File Polygon_mesh_processing/extrude.cpp

#include <CGAL/Surface_mesh.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/extrude.h>
#include <CGAL/IO/polygon_mesh_io.h>
#include <iostream>
#include <fstream>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef Kernel::Vector_3 Vector;
typedef boost::property_map<Mesh, CGAL::vertex_point_t>::type VPMap;
typedef Mesh::template Property_map<vertex_descriptor, Vector> VNMap;
struct Bottom
{
Bottom(VPMap pmap, VNMap nmap, double vlen)
: pmap(pmap), nmap(nmap), vlen(vlen)
{}
void operator()(const vertex_descriptor& vin, const vertex_descriptor vout) const
{
put(pmap, vout, get(pmap, vout) - vlen* get(nmap, vin));
}
VPMap pmap;
VNMap nmap;
double vlen;
};
struct Top
{
Top(VPMap pmap, VNMap nmap, double vlen)
: pmap(pmap), nmap(nmap), vlen(vlen)
{}
void operator()(const vertex_descriptor& vin, const vertex_descriptor vout) const
{
put(pmap, vout, get(pmap, vout) + vlen* get(nmap, vin));
}
VPMap pmap;
VNMap nmap;
double vlen;
};
int main(int argc, char* argv[])
{
Mesh in, out;
std::string filename = (argc > 1) ? std::string(argv[1]) : CGAL::data_file_path("meshes/cube-ouvert.off");
double vlen = (argc > 2) ? std::stod(argv[2]) : 0.1;
CGAL::IO::read_polygon_mesh(filename, in);
VNMap vnormals = in.template add_property_map<vertex_descriptor, Vector>("v:normal", CGAL::NULL_VECTOR).first;
CGAL::Polygon_mesh_processing::compute_vertex_normals(in, vnormals);
Bottom bottom(get(CGAL::vertex_point,out), vnormals, vlen);
Top top(get(CGAL::vertex_point,out), vnormals, vlen);
CGAL::Polygon_mesh_processing::extrude_mesh(in, out, bottom, top);
filename = filename.substr(filename.find_last_of("/") + 1, filename.length() - 1);
filename = filename.substr(0, filename.find_last_of("."));
filename = filename + "_" + std::to_string(vlen) + ".off";
CGAL::IO::write_polygon_mesh(filename, out, CGAL::parameters::stream_precision(17));
return 0;
}
Kernel::Vector_3
CGAL::Polygon_mesh_processing::extrude_mesh
void extrude_mesh(const InputMesh &input, OutputMesh &output, const BottomFunctor &bot, const TopFunctor &top, const NamedParameters1 &np_in=parameters::default_values(), const NamedParameters2 &np_out=parameters::default_values())
performs a generalized extrusion of input and puts it in output.
Definition: extrude.h:164
CGAL::Polygon_mesh_processing::compute_vertex_normals
void compute_vertex_normals(const PolygonMesh &pmesh, VertexNormalMap vertex_normals, const NamedParameters &np=parameters::default_values())
computes the outward unit vector normal for all vertices of the polygon mesh.
Definition: compute_normal.h:787
CGAL::IO::read_polygon_mesh
bool read_polygon_mesh(const std::string &fname, Graph &g, const NamedParameters &np=parameters::default_values())
CGAL::NULL_VECTOR
const CGAL::Null_vector NULL_VECTOR

# Corefinement and Boolean Operations

## Definitions

Corefinement Given two triangulated surface meshes, the corefinement operation consists in refining both meshes so that their intersection polylines are a subset of edges in both refined meshes.

Figure 75.9 Corefinement of two triangulated surface meshes. (Left) Input meshes; (Right) The two input meshes corefined. The common edges of the two meshes are drawn in green.

Volume bounded by a triangulated surface mesh Given a closed triangulated surface mesh, each connected component splits the 3D space into two subspaces. The vertex sequence of each face of a component is seen either clockwise or counterclockwise from these two subspaces. The subspace that sees the sequence clockwise (resp. counterclockwise) is on the negative (resp. positive) side of the component. Given a closed triangulated surface mesh `tm` with no self-intersections, the connected components of `tm` divide the 3D space into subspaces. We say that `tm` bounds a volume if each subspace lies exclusively on the positive (or negative) side of all the incident connected components of `tm`. The volume bounded by `tm` is the union of all subspaces that are on negative sides of their incident connected components of `tm`.

Figure 75.10 Volumes bounded by a triangulated surface mesh: The figure shows meshes representing three nested spheres (three connected components). The left side of the picture shows a clipped triangulated surface mesh, with the two possible orientations of the faces for which a volume is bounded by the mesh. The positive and negative sides of each connected component is displayed in light and dark blue, respectively. The right part of the picture shows clipped tetrahedral meshes of the corresponding bounded volumes.

## Corefinement

The corefinement of two triangulated surface meshes can be done using the function `CGAL::Polygon_mesh_processing::corefine()`. It takes as input the two triangulated surface meshes to corefine. If constrained edge maps are provided, edges belonging to the intersection of the input meshes will be marked as constrained. In addition, if an edge that was marked as constrained is split during the corefinement, sub-edges will be marked as constrained as well.

## Boolean Operations

Figure 75.11 Let `C` and `S` be the volumes bounded by the triangulated surface meshes of a cube and a sphere, respectively. From left to right, the picture shows the triangulated surface meshes bounding the union of `C` and `S`, `C` minus `S`, the intersection of `C` and `S` and `S` minus `C`.

The corefinement of two triangulated surface meshes can naturally be used for computing Boolean operations on volumes. Considering two triangulated surface meshes, each bounding a volume, the functions `CGAL::Polygon_mesh_processing::corefine_and_compute_union()`, `CGAL::Polygon_mesh_processing::corefine_and_compute_intersection()` and `CGAL::Polygon_mesh_processing::corefine_and_compute_difference()` respectively compute the union, the intersection and the difference of the two volumes. If several Boolean operations must be computed at the same time, the function `corefine_and_compute_boolean_operations()` should be used.

There is no restriction on the topology of the input volumes. However, there are some requirements on the input to guarantee that the operation is possible. First, the input meshes must not self-intersect. Second, the operation is possible only if the output can be bounded by a manifold triangulated surface mesh. In particular this means that the output volume has no part with zero thickness. Mathematically speaking, the intersection with an infinitesimally small ball centered in the output volume is a topological ball. At the surface level this means that no non-manifold vertex or edge is allowed in the output. For example, it is not possible to compute the union of two cubes that are disjoint but sharing an edge. In case you have to deal with such scenarios, you should consider using the package 3D Boolean Operations on Nef Polyhedra.

It is possible to update the input so that it contains the result (in-place operation). In that case the whole mesh will not be copied and only the region around the intersection polyline will be modified. In case the Boolean operation is not possible, the input mesh will nevertheless be corefined.

## Kernel and Validity of the Output

The corefinement operation (which is also internally used in the three Boolean operations) will correctly change the topology of the input surface mesh if the point type used in the point property maps of the input meshes is from a CGAL Kernel with exact predicates. If that kernel does not have exact constructions, the embedding of the output surface mesh might have self-intersections. In case of consecutive operations, it is thus recommended to use a point property map with points from a kernel with exact predicates and exact constructions (such as `CGAL::Exact_predicates_exact_constructions_kernel`).

In practice, this means that with exact predicates and inexact constructions, edges will be split at each intersection with a triangle but the position of the intersection point might create self-intersections due to the limited precision of floating point numbers.

## Clipping

As a natural extension, some clipping functionalities with a volume bounded by a closed mesh and a halfspace (defined by the negative side of a plane to be consistent with the outward normal convention) are offered. The functions `CGAL::Polygon_mesh_processing::clip()` and `CGAL::Polygon_mesh_processing::split()` have some options to select whether the clipping should be done at the volume or surface level, and also if the clipper should be considered as compact or not. This is illustrated on Figure 75.12 and Figure 75.13.

Figure 75.12 Clipping a cube with a halfspace. From left to right: (i) initial cube and the plane defining the clipping halfspace; (ii) `clip_volume=false`: clipping of the surface mesh (boundary edges depicted in red); (iii) `clip_volume=true`: clipping of the volume bounded by the surface mesh.

Figure 75.13 Clipping a cube with a halfspace: compactivity of the clipper (`clip_volume=false` in both cases). From left to right: (i) initial cube and the plane defining the clipping halfspace, note that a whole face of the cube (2 triangles) is exactly contained in the plane; (ii) `use_compact_clipper=true`: clipping of the surface mesh with a compact halfspace: coplanar faces are part of the output; (iii) `use_compact_clipper=false`: clipping of the surface mesh with a non-compact halfspace: coplanar faces are not part of the output.

## Examples

### Computing the Union of Two Volumes

File Polygon_mesh_processing/corefinement_mesh_union.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <iostream>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off");
const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off");
Mesh mesh1, mesh2;
if(!PMP::IO::read_polygon_mesh(filename1, mesh1) || !PMP::IO::read_polygon_mesh(filename2, mesh2))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
Mesh out;
bool valid_union = PMP::corefine_and_compute_union(mesh1,mesh2, out);
if(valid_union)
{
std::cout << "Union was successfully computed\n";
CGAL::IO::write_polygon_mesh("union.off", out, CGAL::parameters::stream_precision(17));
return 0;
}
std::cout << "Union could not be computed\n";
return 1;
}

### Boolean Operation and Local Remeshing

This example is similar to the previous one, but here we subtract a volume and update the first input triangulated surface mesh (in-place operation). The edges that are on the intersection of the input meshes are marked and the region around them is remeshed isotropically while preserving the intersection polyline.
File Polygon_mesh_processing/corefinement_difference_remeshed.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/boost/graph/selection.h>
#include <fstream>
#include <iostream>
#include <iterator>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
typedef boost::graph_traits<Mesh>::edge_descriptor edge_descriptor;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
namespace PMP = CGAL::Polygon_mesh_processing;
namespace params = CGAL::parameters;
struct Vector_pmap_wrapper
{
std::vector<bool>& vect;
Vector_pmap_wrapper(std::vector<bool>& v) : vect(v) {}
friend bool get(const Vector_pmap_wrapper& m, face_descriptor f)
{
return m.vect[f];
}
friend void put(const Vector_pmap_wrapper& m, face_descriptor f, bool b)
{
m.vect[f]=b;
}
};
int main(int argc, char* argv[])
{
const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off");
const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off");
Mesh mesh1, mesh2;
if(!PMP::IO::read_polygon_mesh(filename1, mesh1) || !PMP::IO::read_polygon_mesh(filename2, mesh2))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
//create a property on edges to indicate whether they are constrained
Mesh::Property_map<edge_descriptor,bool> is_constrained_map =
mesh1.add_property_map<edge_descriptor,bool>("e:is_constrained", false).first;
// update mesh1 to contain the mesh bounding the difference
// of the two input volumes.
bool valid_difference =
PMP::corefine_and_compute_difference(mesh1,
mesh2,
mesh1,
params::default_values(), // default parameters for mesh1
params::default_values(), // default parameters for mesh2
params::edge_is_constrained_map(is_constrained_map));
if (valid_difference)
{
std::cout << "Difference was successfully computed\n";
CGAL::IO::write_polygon_mesh("difference.off", mesh1, CGAL::parameters::stream_precision(17));
}
else
{
std::cout << "Difference could not be computed\n";
return 1;
}
// collect faces incident to a constrained edge
std::vector<face_descriptor> selected_faces;
std::vector<bool> is_selected(num_faces(mesh1), false);
for(edge_descriptor e : edges(mesh1))
if (is_constrained_map[e])
{
// insert all faces incident to the target vertex
for(halfedge_descriptor h : halfedges_around_target(halfedge(e,mesh1),mesh1))
{
if (!is_border(h, mesh1) )
{
face_descriptor f=face(h, mesh1);
if ( !is_selected[f] )
{
selected_faces.push_back(f);
is_selected[f]=true;
}
}
}
}
// increase the face selection
CGAL::expand_face_selection(selected_faces, mesh1, 2,
Vector_pmap_wrapper(is_selected), std::back_inserter(selected_faces));
std::cout << selected_faces.size()
<< " faces were selected for the remeshing step\n";
// remesh the region around the intersection polylines
PMP::isotropic_remeshing(selected_faces, 0.02, mesh1,
params::edge_is_constrained_map(is_constrained_map));
CGAL::IO::write_polygon_mesh("difference_remeshed.off", mesh1, CGAL::parameters::stream_precision(17));
return 0;
}
CGAL::is_border
bool is_border(typename boost::graph_traits< FaceGraph >::halfedge_descriptor hd, const FaceGraph &g)
CGAL::halfedges_around_target
Iterator_range< Halfedge_around_target_iterator< Graph > > halfedges_around_target(typename boost::graph_traits< Graph >::halfedge_descriptor h, const Graph &g)
CGAL::expand_face_selection
OutputIterator expand_face_selection(const FaceRange &selection, FaceGraph &fg, unsigned int k, IsFaceSelectedPMap is_selected, OutputIterator out)
CGAL::parameters

### Robustness of Consecutive Operations

This example computes the intersection of two volumes and then does the union of the result with one of the input volumes. This operation is in general not possible when using inexact constructions. Instead of using a mesh with a point from a kernel with exact constructions, the exact points are a property of the mesh vertices that we can reuse in a later operations. With that property, we can manipulate a mesh with points having floating point coordinates but benefit from the robustness provided by the exact constructions.
File Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <iostream>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_predicates_exact_constructions_kernel EK;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef Mesh::Property_map<vertex_descriptor,EK::Point_3> Exact_point_map;
namespace PMP = CGAL::Polygon_mesh_processing;
namespace params = CGAL::parameters;
struct Exact_vertex_point_map
{
// typedef for the property map
typedef boost::property_traits<Exact_point_map>::value_type value_type;
typedef boost::property_traits<Exact_point_map>::reference reference;
typedef boost::property_traits<Exact_point_map>::key_type key_type;
typedef boost::read_write_property_map_tag category;
// exterior references
Exact_point_map exact_point_map;
Mesh* tm_ptr;
// Converters
CGAL::Cartesian_converter<K, EK> to_exact;
CGAL::Cartesian_converter<EK, K> to_input;
Exact_vertex_point_map()
: tm_ptr(nullptr)
{}
Exact_vertex_point_map(const Exact_point_map& ep, Mesh& tm)
: exact_point_map(ep)
, tm_ptr(&tm)
{
for (Mesh::Vertex_index v : vertices(tm))
exact_point_map[v]=to_exact(tm.point(v));
}
friend
reference get(const Exact_vertex_point_map& map, key_type k)
{
CGAL_precondition(map.tm_ptr!=nullptr);
return map.exact_point_map[k];
}
friend
void put(const Exact_vertex_point_map& map, key_type k, const EK::Point_3& p)
{
CGAL_precondition(map.tm_ptr!=nullptr);
map.exact_point_map[k]=p;
// create the input point from the exact one
map.tm_ptr->point(k)=map.to_input(p);
}
};
int main(int argc, char* argv[])
{
const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off");
const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off");
Mesh mesh1, mesh2;
if(!PMP::IO::read_polygon_mesh(filename1, mesh1) || !PMP::IO::read_polygon_mesh(filename2, mesh2))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
Exact_point_map mesh1_exact_points =
mesh1.add_property_map<vertex_descriptor,EK::Point_3>("v:exact_point").first;
Exact_point_map mesh2_exact_points =
mesh2.add_property_map<vertex_descriptor,EK::Point_3>("v:exact_point").first;
Exact_vertex_point_map mesh1_vpm(mesh1_exact_points, mesh1);
Exact_vertex_point_map mesh2_vpm(mesh2_exact_points, mesh2);
if ( PMP::corefine_and_compute_intersection(mesh1,
mesh2,
mesh1,
params::vertex_point_map(mesh1_vpm),
params::vertex_point_map(mesh2_vpm),
params::vertex_point_map(mesh1_vpm) ) )
{
if ( PMP::corefine_and_compute_union(mesh1,
mesh2,
mesh2,
params::vertex_point_map(mesh1_vpm),
params::vertex_point_map(mesh2_vpm),
params::vertex_point_map(mesh2_vpm) ) )
{
std::cout << "Intersection and union were successfully computed\n";
CGAL::IO::write_polygon_mesh("inter_union.off", mesh2, CGAL::parameters::stream_precision(17));
return 0;
}
std::cout << "Union could not be computed\n";
return 1;
}
std::cout << "Intersection could not be computed\n";
return 1;
}
CGAL::Cartesian_converter
CGAL::Exact_predicates_exact_constructions_kernel

# Hole Filling

This package provides an algorithm for filling one closed hole that is either in a triangulated surface mesh or defined by a sequence of points that describe a polyline. The main steps of the algorithm are described in [14] and can be summarized as follows.

First, the largest patch triangulating the boundary of the hole is generated without introducing any new vertex. The patch is selected so as to minimize a quality function evaluated for all possible triangular patches. The quality function first minimizes the worst dihedral angle between patch triangles, then the total surface area of the patch as a tiebreaker. Following the suggestions in [21], the performance of the algorithm is significantly improved by narrowing the search space to faces of a 3D Delaunay triangulation of the hole boundary vertices, from all possible patches, while searching for the best patch with respect to the aforementioned quality criteria.

For some complicated input hole boundary, the generated patch may have self-intersections. After hole filling, the generated patch can be refined and faired using the meshing functions `CGAL::Polygon_mesh_processing::refine()` and `CGAL::Polygon_mesh_processing::fair()` described in Section Meshing.

Figure 75.14 Results of the main steps of the algorithm. From left to right: (a) the hole, (b) the hole after its triangulation, (c) after triangulation and refinement, (d) after triangulation, refinement and fairing.

## API

This package provides four functions for hole filling:

- `triangulate_hole_polyline()` : given a sequence of points defining the hole, triangulates the hole.
- `triangulate_hole()` : given a border halfedge on the boundary of the hole on a mesh, triangulates the hole.
- `triangulate_and_refine_hole()` : in addition to `triangulate_hole()` the generated patch is refined.
- `triangulate_refine_and_fair_hole()` : in addition to `triangulate_and_refine_hole()` the generated patch is also faired.

## Examples

### Triangulate a Polyline

The following example triangulates a hole described by an input polyline.

File Polygon_mesh_processing/triangulate_polyline_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
#include <CGAL/utility.h>
#include <vector>
#include <iterator>
#include <cassert>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
int main()
{
std::vector<Point> polyline;
polyline.push_back(Point( 1.,0.,0.));
polyline.push_back(Point( 0.,1.,0.));
polyline.push_back(Point(-1.,0.,0.));
polyline.push_back(Point( 1.,1.,0.));
// repeating first point (i.e. polyline.push_back(Point(1.,0.,0.)) ) is optional
// any type, having Type(int, int, int) constructor available, can be used to hold output triangles
typedef CGAL::Triple<int, int, int> Triangle_int;
std::vector<Triangle_int> patch;
patch.reserve(polyline.size() -2); // there will be exactly n-2 triangles in the patch
CGAL::Polygon_mesh_processing::triangulate_hole_polyline(
polyline,
std::back_inserter(patch));
for(std::size_t i = 0; i < patch.size(); ++i)
{
std::cout << "Triangle " << i << ": "
<< patch[i].first << " " << patch[i].second << " " << patch[i].third
<< std::endl;
}
// note that no degenerate triangles are generated in the patch
std::vector<Point> polyline_collinear;
polyline_collinear.push_back(Point(1.,0.,0.));
polyline_collinear.push_back(Point(2.,0.,0.));
polyline_collinear.push_back(Point(3.,0.,0.));
polyline_collinear.push_back(Point(4.,0.,0.));
std::vector<Triangle_int> patch_will_be_empty;
CGAL::Polygon_mesh_processing::triangulate_hole_polyline(polyline_collinear,
back_inserter(patch_will_be_empty));
assert(patch_will_be_empty.empty());
return 0;
}
CGAL::Triple
CGAL::Polygon_mesh_processing::triangulate_hole_polyline
OutputIterator triangulate_hole_polyline(const PointRange1 &points, const PointRange2 &third_points, OutputIterator out, const NamedParameters &np=parameters::default_values())
creates triangles to fill the hole defined by points in the range points.
Definition: triangulate_hole.h:721

### Hole Filling From the Border of the Hole

If the input polygon mesh has a hole or more than one hole, it is possible to iteratively fill them by detecting border edges (i.e. with only one incident non-null face) after each hole filling step.

Holes are filled one after the other, and the process stops when there is no border edge left.

This process is illustrated by the example below, where holes are iteratively filled, refined and faired. Optionally, only holes not exceeding a certain diameter or number of edges can be filled. This example assumes that the mesh is stored in a `Surface_mesh` datastructure. Analogous examples when using the `Polyhedron_3` class and a few others are part of the code base.

File Polygon_mesh_processing/hole_filling_example_SM.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <boost/lexical_cast.hpp>
#include <iostream>
#include <iterator>
#include <string>
#include <tuple>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
namespace PMP = CGAL::Polygon_mesh_processing;
bool is_small_hole(halfedge_descriptor h, Mesh & mesh,
double max_hole_diam, int max_num_hole_edges)
{
int num_hole_edges = 0;
CGAL::Bbox_3 hole_bbox;
for (halfedge_descriptor hc : CGAL::halfedges_around_face(h, mesh))
{
const Point& p = mesh.point(target(hc, mesh));
hole_bbox += p.bbox();
++num_hole_edges;
// Exit early, to avoid unnecessary traversal of large holes
if (num_hole_edges > max_num_hole_edges) return false;
if (hole_bbox.xmax() - hole_bbox.xmin() > max_hole_diam) return false;
if (hole_bbox.ymax() - hole_bbox.ymin() > max_hole_diam) return false;
if (hole_bbox.zmax() - hole_bbox.zmin() > max_hole_diam) return false;
}
return true;
}
// Incrementally fill the holes that are no larger than given diameter
// and with no more than a given number of edges (if specified).
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/mech-holes-shark.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
// Both of these must be positive in order to be considered
double max_hole_diam = (argc > 2) ? boost::lexical_cast<double>(argv[2]): -1.0;
int max_num_hole_edges = (argc > 3) ? boost::lexical_cast<int>(argv[3]) : -1;
unsigned int nb_holes = 0;
std::vector<halfedge_descriptor> border_cycles;
// collect one halfedge per boundary cycle
PMP::extract_boundary_cycles(mesh, std::back_inserter(border_cycles));
for(halfedge_descriptor h : border_cycles)
{
if(max_hole_diam > 0 && max_num_hole_edges > 0 &&
!is_small_hole(h, mesh, max_hole_diam, max_num_hole_edges))
continue;
std::vector<face_descriptor> patch_facets;
std::vector<vertex_descriptor> patch_vertices;
bool success = std::get<0>(PMP::triangulate_refine_and_fair_hole(mesh,
h,
CGAL::parameters::face_output_iterator(std::back_inserter(patch_facets))
.vertex_output_iterator(std::back_inserter(patch_vertices))));
std::cout << "* Number of facets in constructed patch: " << patch_facets.size() << std::endl;
std::cout << " Number of vertices in constructed patch: " << patch_vertices.size() << std::endl;
std::cout << " Is fairing successful: " << success << std::endl;
++nb_holes;
}
std::cout << std::endl;
std::cout << nb_holes << " holes have been filled" << std::endl;
CGAL::IO::write_polygon_mesh("filled_SM.off", mesh, CGAL::parameters::stream_precision(17));
std::cout << "Mesh written to: filled_SM.off" << std::endl;
return 0;
}
CGAL::Bbox_3
CGAL::Bbox_3::ymin
double ymin() const
CGAL::Bbox_3::xmax
double xmax() const
CGAL::Bbox_3::zmin
double zmin() const
CGAL::Bbox_3::zmax
double zmax() const
CGAL::Bbox_3::ymax
double ymax() const
CGAL::Bbox_3::xmin
double xmin() const
CGAL::halfedges_around_face
Iterator_range< Halfedge_around_face_iterator< Graph > > halfedges_around_face(typename boost::graph_traits< Graph >::halfedge_descriptor h, const Graph &g)

Figure 75.15 Holes in the fork model are filled with triangle patches.

An additional parameter, named `visitor` can be used to track which phase of the algorithm is being run. This can allow a user to stop the algorithm if a timeout needs to be implemented. For an example, see Polygon_mesh_processing/hole_filling_visitor_example.cpp .

## Performance

The hole filling algorithm has a complexity which depends on the number of vertices. While [14] has a running time of \(O(n^3)\) , [21] in most cases has running time of \(O(n \log n)\). We benchmarked the function `triangulate_refine_and_fair_hole()` for the two meshes below (as well as two more meshes with smaller holes). The machine used was a PC running Windows 10 with an Intel Core i7 CPU clocked at 2.70 GHz. The program was compiled with the Visual C++ 2013 compiler with the O2 option, which maximizes speed.

Figure 75.16 The elephant on the left/right has a hole with 963/7657 vertices.

The following running times were observed:

# vertices | without Delaunay (sec.) | with Delaunay (sec.)
--- | --- | ---
565 | 8.5 | 0.03
774 | 21 | 0.035
967 | 43 | 0.06
7657 | na | 0.4

# Predicates

This package provides several predicates to be evaluated with respect to a triangle mesh.

## Intersections Detection

Intersection tests between triangle meshes and/or polylines can be done using `CGAL::Polygon_mesh_processing::do_intersect()` . Additionally, the function `CGAL::Polygon_mesh_processing::intersecting_meshes()` records all pairs of intersecting meshes in a range.

### Self Intersections

Self intersections within a triangle mesh can be detected by calling the function `CGAL::Polygon_mesh_processing::does_self_intersect()`. Additionally, the function `CGAL::Polygon_mesh_processing::self_intersections()` reports all pairs of intersecting triangles.

### Self Intersections Example

The following example illustrates the detection of self intersection in the `pig.off` mesh. The detected self-intersection is illustrated on Figure 75.17.

File Polygon_mesh_processing/self_intersections_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Real_timer.h>
#include <CGAL/tags.h>
#include <iostream>
#include <iterator>
#include <string>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/pig.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
std::cout << "Using parallel mode? " << std::is_same<CGAL::Parallel_if_available_tag, CGAL::Parallel_tag>::value << std::endl;
CGAL::Real_timer timer;
timer.start();
bool intersecting = PMP::does_self_intersect<CGAL::Parallel_if_available_tag>(mesh, CGAL::parameters::vertex_point_map(get(CGAL::vertex_point, mesh)));
std::cout << (intersecting ? "There are self-intersections." : "There is no self-intersection.") << std::endl;
std::cout << "Elapsed time (does self intersect): " << timer.time() << std::endl;
timer.reset();
std::vector<std::pair<face_descriptor, face_descriptor> > intersected_tris;
PMP::self_intersections<CGAL::Parallel_if_available_tag>(faces(mesh), mesh, std::back_inserter(intersected_tris));
std::cout << intersected_tris.size() << " pairs of triangles intersect." << std::endl;
std::cout << "Elapsed time (self intersections): " << timer.time() << std::endl;
return EXIT_SUCCESS;
}

Figure 75.17 Detecting self-intersections on a triangle mesh. The intersecting triangles are displayed in dark grey and red on the right image.

## Side of Triangle Mesh

The class `CGAL::Side_of_triangle_mesh` provides a functor that tests whether a query point is inside, outside, or on the boundary of the domain bounded by a given closed triangle mesh.

A point is said to be on the bounded side of the domain bounded by the input triangle mesh if an odd number of surfaces is crossed when walking from the point to infinity. The input triangle mesh is expected to contain no self-intersections and to be free from self-inclusions.

The algorithm can handle the case of a triangle mesh with several connected components, and returns correct results. In case of self-inclusions, the ray intersections parity test is performed, and the execution will not fail. However, the user should be aware that the predicate alternately considers sub-volumes to be on the bounded and unbounded sides of the input triangle mesh.

### Inside Test Example

File Polygon_mesh_processing/point_inside_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Side_of_triangle_mesh.h>
#include <iostream>
#include <limits>
#include <string>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point;
typedef CGAL::Polyhedron_3<K> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
double max_coordinate(const Mesh& poly)
{
double max_coord = -std::numeric_limits<double>::infinity();
for(Mesh::Vertex_handle v : vertices(poly))
{
Point p = v->point();
max_coord = (std::max)(max_coord, CGAL::to_double(p.x()));
max_coord = (std::max)(max_coord, CGAL::to_double(p.y()));
max_coord = (std::max)(max_coord, CGAL::to_double(p.z()));
}
return max_coord;
}
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off");
Mesh poly;
if(!PMP::IO::read_polygon_mesh(filename, poly) || CGAL::is_empty(poly) || !CGAL::is_triangle_mesh(poly))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
CGAL::Side_of_triangle_mesh<Mesh, K> inside(poly);
double size = max_coordinate(poly);
unsigned int nb_points = 100;
std::vector<Point> points;
points.reserve(nb_points);
CGAL::Random_points_in_cube_3<Point> gen(size);
for (unsigned int i = 0; i < nb_points; ++i)
points.push_back(*gen++);
std::cout << "Test " << nb_points << " random points in cube "
<< "[-" << size << "; " << size <<"]" << std::endl;
int nb_inside = 0;
int nb_boundary = 0;
for (std::size_t i = 0; i < nb_points; ++i)
{
CGAL::Bounded_side res = inside(points[i]);
if (res == CGAL::ON_BOUNDED_SIDE) { ++nb_inside; }
if (res == CGAL::ON_BOUNDARY) { ++nb_boundary; }
}
std::cerr << "Total query size: " << points.size() << std::endl;
std::cerr << " " << nb_inside << " points inside " << std::endl;
std::cerr << " " << nb_boundary << " points on boundary " << std::endl;
std::cerr << " " << points.size() - nb_inside - nb_boundary << " points outside " << std::endl;
return 0;
}
CGAL::Random_points_in_cube_3
CGAL::Side_of_triangle_mesh
This class provides an efficient point location functionality with respect to a domain bounded by one...
Definition: Side_of_triangle_mesh.h:76
CGAL::to_double
double to_double(const NT &x)
CGAL::Bounded_side
Bounded_side
CGAL::ON_BOUNDARY
ON_BOUNDARY
CGAL::ON_BOUNDED_SIDE
ON_BOUNDED_SIDE

## Polyhedral Envelope Containment Check

The class `CGAL::Polyhedral_envelope` provides functors to check if a query point, segment, or triangle is fully contained in a polyhedral envelope of a triangle mesh or of a triangle soup. In the following, input triangles will refer either to the triangles of a mesh or of a soup.

The polyhedral envelope is a conservative approximation of the Minkowski sum envelope of a set of triangles with a sphere of radius \( \epsilon \). The latter has cylindrical and spherical patches at convex edges and vertices of the input triangles.

Given a distance \( \delta = \epsilon / \sqrt(3)\) we can associate a prism to each triangle by intersecting two halfspaces parallel to the triangle, three halfspaces orthogonal to the triangle and parallel to the edges, and additionally halfspaces for clipping obtuse angles, with the face normal corresponding to the bisector of the angle. These halfspaces are at distance \( \delta \) and such that they contain the triangle.

The polyhedral envelope of a set of triangles with a tolerance \( \epsilon \) then is the union of the prisms of all faces with \( \delta = \epsilon / \sqrt(3) \).

Figure 75.18 The prism for a single triangle, the polyhedral envelope as well as the Minkowski sum envelope for a triangle mesh.

The polyhedral envelope is guaranteed to be inside the Minkowski sum envelope. This containment test is exact for the polyhedral envelope, and conservative for the Minkowswi sum envelope: If a query is inside the polyhedral envelope, we can be sure that it is also in the Minkowski sum envelope, but if it is outside the polyhedral envelope we do not know where it is with respect to the Minkowski sum envelope.

The algorithm of Wang et al. [20] for the polyhedral envelope containment check can be summarized as follows. The prisms of the faces of the input triangles are stored in an AABB tree, which is used to quickly identify the prisms whose bounding box overlaps with the query.

For a query point the algorithm checks if it is inside one of these prisms. For a query segment or triangle the algorithms checks if the query is completely covered. The details of how to check this covering can be found in the paper.

The polyhedral envelope containment check is used by the class `Surface_mesh_simplification::Polyhedral_envelope_filter` of the package Triangulated Surface Mesh Simplification, in order to simplify a triangle mesh within a given tolerance.

### Polyhedral Envelope Example

The following example shows how to construct the polyhedral envelope for a `Surface_mesh` and perform queries.

File Polygon_mesh_processing/polyhedral_envelope.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedral_envelope.h>
#include <CGAL/Surface_mesh.h>
#include <iostream>
#include <fstream>
int main(int argc, char* argv[])
{
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;
typedef CGAL::Surface_mesh<Point_3> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef CGAL::Polyhedral_envelope<Kernel> Envelope;
std::ifstream in((argc>1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"));
Mesh tmesh;
in >> tmesh;
double eps = (argc>2) ? std::stod(std::string(argv[2])) : 0.2;
Envelope envelope(tmesh, eps);
int i = (argc>3) ? std::stoi(std::string(argv[3])) : 0;
int j = (argc>4) ? std::stoi(std::string(argv[4])) : 100;
int k = (argc>5) ? std::stoi(std::string(argv[5])) : 200;
if(envelope(tmesh.point(vertex_descriptor(i)),
tmesh.point(vertex_descriptor(j)),
tmesh.point(vertex_descriptor(k)))){
std::cout << "inside polyhedral envelope" << std::endl;
}
return 0;
}
CGAL::Polyhedral_envelope
This class can be used to check if a query point, segment, or triangle is inside or outside a polyhed...
Definition: Polyhedral_envelope.h:105

As the polyhedral envelope does not need any connectivity information, the same check can be performed just with a triangle soup.

File Polygon_mesh_processing/polyhedral_envelope_of_triangle_soup.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedral_envelope.h>
#include <CGAL/IO/OFF.h>
#include <iostream>
#include <fstream>
#include <vector>
int main(int argc, char* argv[])
{
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;
typedef CGAL::Polyhedral_envelope<Kernel> Envelope;
std::ifstream in((argc>1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"));
double eps = (argc>2) ? std::stod(std::string(argv[2])) : 0.2;
std::vector<Point_3> points;
std::vector<std::vector<std::size_t> > polygons;
CGAL::IO::read_OFF(in, points, polygons);
Envelope envelope(points, polygons, eps);
int i = (argc>3) ? std::stoi(std::string(argv[3])) : 0;
int j = (argc>4) ? std::stoi(std::string(argv[4])) : 100;
int k = (argc>5) ? std::stoi(std::string(argv[5])) : 200;
if (envelope(points[i], points[j],points[k]))
{
std::cout << "inside polyhedral envelope" << std::endl;
}
return 0;
}
CGAL::IO::read_OFF
bool read_OFF(std::istream &is, Graph &g, const NamedParameters &np=parameters::default_values())

A triangle mesh can also be used as query, to check if a remeshed version is within the polyhedral envelope of an input mesh.

File Polygon_mesh_processing/polyhedral_envelope_mesh_containment.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedral_envelope.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Surface_mesh.h>
#include <algorithm>
#include <iostream>
#include <fstream>
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;
typedef CGAL::Surface_mesh<Point_3> Mesh;
typedef CGAL::Polyhedral_envelope<Kernel> Envelope;
std::ifstream in((argc>1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"));
Mesh tmesh;
in >> tmesh;
// remesh the input using the longest edge size as target edge length
Mesh query = tmesh;
Mesh::Edge_iterator longest_edge_it =
std::max_element(edges(query).begin(), edges(query).end(),
[&query](Mesh::Edge_index e1, Mesh::Edge_index e2)
{
return PMP::edge_length(halfedge(e1, query), query) <
PMP::edge_length(halfedge(e2, query), query);
});
PMP::isotropic_remeshing(faces(tmesh), PMP::edge_length(halfedge(*longest_edge_it, query), query), query);
// construct the polyhedral envelope
const double eps = (argc>2) ? std::stod(std::string(argv[2])) : 0.01;
Envelope envelope(tmesh, eps);
// check is the remeshed version is inside the polyhedral envelope of the input mesh
if ( envelope(query) )
std::cout << "Remeshing is inside the polyhedral envelope\n";
else
std::cout << "Remeshing is not inside the polyhedral envelope\n";
std::ofstream("remeshed.off") << query;
return 0;
}

## Shape Predicates

Badly shaped or, even worse, completely degenerate elements of a polygon mesh are problematic in many algorithms which one might want to use on the mesh. This package offers a toolkit of functions to detect such undesirable elements.

- `CGAL::Polygon_mesh_processing::is_degenerate_edge()`, to detect if an edge is degenerate (that is, if its two vertices share the same geometric location).
- `CGAL::Polygon_mesh_processing::is_degenerate_triangle_face()`, to detect if a face is degenerate (that is, if its three vertices are collinear).
- `CGAL::Polygon_mesh_processing::degenerate_edges()`, to collect degenerate edges within a range of edges.
- `CGAL::Polygon_mesh_processing::degenerate_faces()`, to collect degenerate faces within a range of faces.
- `CGAL::Polygon_mesh_processing::is_cap_triangle_face()`
- `CGAL::Polygon_mesh_processing::is_needle_triangle_face()`

## Surface Location Functions

To ease the manipulation of points on a surface, CGAL offers a multitude of functions based upon a different representation of a point on a polygon mesh: the point is represented as a pair of a face of the polygon mesh and a triplet of barycentric coordinates. This definition enables a robust handling of polylines between points living in the same face: for example, two 3D segments created by four points within the same face that should intersect might not actually intersect due to inexact computations. However, manipulating these same points through their barycentric coordinates can instead be done, and intersections computed in the barycentric space will not suffer from the same issues. Furthermore, this definition is only dependent on the intrinsic dimension of the surface (i.e. 2) and not on the ambient dimension within which the surface is embedded.

The functions of the group Location Functions offer the following functionalities: location computations (`CGAL::Polygon_mesh_processing::locate()`, and similar) given a point, finding the nearest point on a mesh given a point or a ray (`CGAL::Polygon_mesh_processing::locate_with_AABB_tree()`, and similar), and location-based predicates (for example, `CGAL::Polygon_mesh_processing::is_on_face_border()`).

The example Polygon_mesh_processing/locate_example.cpp presents a few of these functions.

# Orientation

This package offers multiple functions to compute consistent face orientations for set of faces (Section Polygon Soups) and polygon meshes (Section Polygon Meshes). Section Orientation Examples offers an example of combination of these functions.

## Polygon Soups

When the faces of a polygon mesh are given but the connectivity is unknown, this set of faces is called a polygon soup.

Before running any of the algorithms on a polygon soup, one should ensure that the polygons are consistently oriented. To do so, this package provides the function `CGAL::Polygon_mesh_processing::orient_polygon_soup()`, described in [1].

To deal with polygon soups that cannot be converted to a combinatorially manifold surface, some points must be duplicated. Because a polygon soup does not have any connectivity (each point has as many occurrences as the number of polygons it belongs to), duplicating one point (or a pair of points) amounts to duplicating the polygon to which it belongs. The duplicated points are either an endpoint of an edge incident to more than two polygons, an endpoint of an edge between two polygons with incompatible orientations (during the re-orientation process), or more generally a point p at which the intersection of an infinitesimally small ball centered at p with the polygons incident to it is not a topological disk.

Once the polygon soup is consistently oriented, with possibly duplicated (or more) points, the connectivity can be recovered and made consistent to build a valid polygon mesh. The function `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()` performs this mesh construction step.

Inversely, a polygon soup can be constructed from a polygon mesh, using the function `CGAL::Polygon_mesh_processing::polygon_mesh_to_polygon_soup()`.

## Polygon Meshes

This package provides functions dealing with the orientation of faces in a closed polygon mesh.

- The function `CGAL::Polygon_mesh_processing::orient()` makes each connected component of a closed polygon mesh outward- or inward-oriented.
- The function `CGAL::Polygon_mesh_processing::orient_to_bound_a_volume()` orients the connected components of a closed polygon mesh so that it bounds a volume (see Definitions for the precise definition).
- The function `CGAL::Polygon_mesh_processing::is_outward_oriented()` checks whether an oriented polygon mesh is oriented such that the normals to all faces are oriented towards the outside of the domain bounded by the input polygon mesh.
- The function `CGAL::Polygon_mesh_processing::reverse_face_orientations()` reverses the orientation of halfedges around faces. As a consequence, the normal computed for each face (see Section Computing Normals) is also reversed.
- The function `CGAL::Polygon_mesh_processing::volume_connected_components()` provides information about the 3D arrangement of the surface connected components in a given triangle mesh. It comes with many named parameter options making it also a more general version of `is_outward_oriented()`.
- The function `CGAL::Polygon_mesh_processing::duplicate_non_manifold_edges_in_polygon_soup()` duplicates points and edges to make a soup orientable, without changing the orientation of the faces.
- The function `CGAL::Polygon_mesh_processing::orient_triangle_soup_with_reference_triangle_mesh()` takes an input mesh as a reference and orients the triangles of a soup according to it.
- The function `CGAL::Polygon_mesh_processing::merge_reversible_connected_components()` merges the connected components of a polygon mesh if possible.

## Orientation Examples

This example shows how to generate a mesh from a polygon soup. The first step is to get a soup of consistently oriented faces, before rebuilding the connectivity. In this example, some orientation tests are performed on the output polygon mesh to illustrate Section Orientation.

File Polygon_mesh_processing/orient_polygon_soup_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedron_items_with_id_3.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Polyhedron_3<K, CGAL::Polyhedron_items_with_id_3> Mesh;
// Optional visitor for orientating a polygon soup to demonstrate usage for some functions.
// inherits from the default class as some functions are not overloaded
struct Visitor : public CGAL::Polygon_mesh_processing::Default_orientation_visitor
{
void non_manifold_edge(std::size_t id1, std::size_t id2, std::size_t nb_poly)
{
std::cout << "The edge " << id1 << ", " << id2 << " is not manifold: " << nb_poly << " incident polygons." << std::endl;
}
void non_manifold_vertex(std::size_t id, std::size_t nb_cycles)
{
std::cout << "The vertex " << id << " is not manifold: " << nb_cycles << " connected components of vertices in the link." << std::endl;
}
void duplicated_vertex(std::size_t v1, std::size_t v2)
{
std::cout << "The vertex " << v1 << " has been duplicated, its new id is " << v2 << "." << std::endl;
}
void vertex_id_in_polygon_replaced(std::size_t p_id, std::size_t i1, std::size_t i2)
{
std::cout << "In the polygon " << p_id << ", the index " << i1 << " has been replaced by " << i2 << "." << std::endl;
}
void polygon_orientation_reversed(std::size_t p_id)
{
std::cout << "The polygon " << p_id << " has been reversed." << std::endl;
}
};
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/tet-shuffled.off");
std::vector<K::Point_3> points;
std::vector<std::vector<std::size_t> > polygons;
if(!CGAL::IO::read_polygon_soup(filename, points, polygons) || points.empty())
{
std::cerr << "Cannot open file " << std::endl;
return EXIT_FAILURE;
}
Visitor visitor;
CGAL::Polygon_mesh_processing::orient_polygon_soup(points, polygons, CGAL::parameters::visitor(visitor));
Mesh mesh;
CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh);
// Number the faces because 'orient_to_bound_a_volume' needs a face <--> index map
int index = 0;
for(Mesh::Face_iterator fb=mesh.facets_begin(), fe=mesh.facets_end(); fb!=fe; ++fb)
fb->id() = index++;
if(CGAL::is_closed(mesh))
CGAL::Polygon_mesh_processing::orient_to_bound_a_volume(mesh);
std::ofstream out("tet-oriented1.off");
out.precision(17);
out << mesh;
out.close();
CGAL::Polygon_mesh_processing::reverse_face_orientations(mesh);
std::ofstream out2("tet-oriented2.off");
out2.precision(17);
out2 << mesh;
out2.close();
return EXIT_SUCCESS;
}
CGAL::IO::read_polygon_soup
bool read_polygon_soup(const std::string &fname, PointRange &points, PolygonRange &polygons, const NamedParameters &np=parameters::default_values())
CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh
void polygon_soup_to_polygon_mesh(const PointRange &points, const PolygonRange &polygons, PolygonMesh &out, const NamedParameters_PS &np_ps=parameters::default_values(), const NamedParameters_PM &np_pm=parameters::default_values())
builds a polygon mesh from a soup of polygons.
Definition: polygon_soup_to_polygon_mesh.h:330
CGAL::Polygon_mesh_processing::orient_polygon_soup
bool orient_polygon_soup(PointRange &points, PolygonRange &polygons, const NamedParameters &np=parameters::default_values())
tries to consistently orient a soup of polygons in 3D space.
Definition: orient_polygon_soup.h:540
CGAL::Polygon_mesh_processing::reverse_face_orientations
void reverse_face_orientations(PolygonMesh &pmesh)
reverses for each face the order of the vertices along the face boundary.
Definition: orientation.h:274
CGAL::Polygon_mesh_processing::orient_to_bound_a_volume
void orient_to_bound_a_volume(TriangleMesh &tm, const NamedParameters &np=parameters::default_values())
orients the connected components of tm to make it bound a volume.
Definition: orientation.h:1347
CGAL::is_closed
bool is_closed(const FaceGraph &g)
CGAL::Polygon_mesh_processing::Default_orientation_visitor
Default visitor model of PMPPolygonSoupOrientationVisitor.
Definition: orient_polygon_soup.h:47

This example shows how to correctly repair and orient a soup to get a mesh from a reference :

File Polygon_mesh_processing/orientation_pipeline_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup_extension.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/algorithm.h>
#include <CGAL/Timer.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <algorithm>
#include <cstdlib>
#include <iostream>
#include <string>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef CGAL::Surface_mesh<Point_3> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char** argv)
{
const std::string input_filename = (argc < 2) ? CGAL::data_file_path("meshes/blobby-shuffled.off") : argv[1];
const std::string reference_filename = (argc < 2) ? CGAL::data_file_path("meshes/blobby.off") : argv[2];
std::vector<Point_3> points;
std::vector<std::vector<std::size_t> > polygons;
if(!CGAL::IO::read_polygon_soup(input_filename, points, polygons) ||
points.size() == 0 || polygons.size() == 0)
{
std::cerr << "Error: can not read input file.\n";
return 1;
}
Mesh ref1;
if(!PMP::IO::read_polygon_mesh(reference_filename, ref1))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
std::cout << "Is the soup a polygon mesh ? : " << PMP::is_polygon_soup_a_polygon_mesh(polygons) << std::endl;
PMP::orient_triangle_soup_with_reference_triangle_mesh<CGAL::Sequential_tag>(ref1, points, polygons);
std::cout << "And now ? : " << PMP::is_polygon_soup_a_polygon_mesh(polygons) << std::endl;
PMP::duplicate_non_manifold_edges_in_polygon_soup(points, polygons);
std::cout << "And now ? : " << PMP::is_polygon_soup_a_polygon_mesh(polygons) << std::endl;
Mesh poly;
PMP::polygon_soup_to_polygon_mesh(points, polygons, poly);
typedef boost::property_map<Mesh, CGAL::dynamic_face_property_t<std::size_t> >::type Fccmap;
Fccmap fccmap = get(CGAL::dynamic_face_property_t<std::size_t>(), poly);
std::cout << PMP::connected_components(poly, fccmap) << " CCs before merge." << std::endl;
PMP::merge_reversible_connected_components(poly);
std::cout<<PMP::connected_components(poly, fccmap) << " remaining CCs." << std::endl;
return 0;
}
CGAL::dynamic_face_property_t

# Combinatorial Repair

## Polygon Soup Repairing

To ensure that a polygon soup can be oriented (see Section Polygon Soups) and transformed into a workable polygon mesh, it might be necessary to preprocess the data to remove combinatorial and geometrical errors. This package offers the following functions:

- `CGAL::Polygon_mesh_processing::merge_duplicate_points_in_polygon_soup()`,
- `CGAL::Polygon_mesh_processing::merge_duplicate_polygons_in_polygon_soup()`,
- `CGAL::Polygon_mesh_processing::remove_isolated_points_in_polygon_soup()`,

as well as the function `CGAL::Polygon_mesh_processing::repair_polygon_soup()`, which bundles the previous functions and an additional handful of repairing techniques to obtain an as-clean-as-possible polygon soup.

## Stitching

When handling polygon meshes, it might happen that a mesh has several edges and vertices that are duplicated. For those edges and vertices, the connectivity of the mesh is incomplete, if not considered incorrect.

Stitching the borders of a polygon mesh can be done to fix some of the duplication. It consists in two main steps. First, border edges that are geometrically identical but duplicated are detected and paired. Then, they are "stitched" together so that edges and vertices duplicates are removed from the mesh, and each of these remaining edges is incident to exactly two faces.

The functions `CGAL::Polygon_mesh_processing::stitch_boundary_cycle()`, `CGAL::Polygon_mesh_processing::stitch_boundary_cycles()`, and `CGAL::Polygon_mesh_processing::stitch_borders()` can perform such repairing operations: the first two functions can be used to stitch halfedges that are part of the same boundary(ies), whereas the third function is more generic and can also stitch halfedges that live on different borders.

The input mesh should be manifold; otherwise, stitching is not guaranteed to succeed.

### Stitching Example

The following example applies the stitching operation to a simple quad mesh with duplicated border edges.

File Polygon_mesh_processing/stitch_borders_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/stitch_borders.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <iostream>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Polyhedron_3<K> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/quads_to_stitch.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
std::cout << "Before stitching : " << std::endl;
std::cout << "\t Number of vertices :\t" << mesh.size_of_vertices() << std::endl;
std::cout << "\t Number of halfedges :\t" << mesh.size_of_halfedges() << std::endl;
std::cout << "\t Number of facets :\t" << mesh.size_of_facets() << std::endl;
PMP::stitch_borders(mesh);
std::cout << "Stitching done : " << std::endl;
std::cout << "\t Number of vertices :\t" << mesh.size_of_vertices() << std::endl;
std::cout << "\t Number of halfedges :\t" << mesh.size_of_halfedges() << std::endl;
std::cout << "\t Number of facets :\t" << mesh.size_of_facets() << std::endl;
CGAL::IO::write_polygon_mesh("mesh_stitched.off", mesh, CGAL::parameters::stream_precision(17));
return 0;
}
CGAL::Polyhedron_3::size_of_vertices
size_type size_of_vertices() const

## Polygon Mesh Manifoldness

Non-manifold vertices can be detected using the function `CGAL::Polygon_mesh_processing::is_non_manifold_vertex()`. The function `CGAL::Polygon_mesh_processing::duplicate_non_manifold_vertices()` can be used to attempt to create a combinatorially manifold surface mesh by splitting any non-manifold vertex into as many vertices as there are manifold sheets at this geometric position. Note however that the mesh will still not be manifold from a geometric point of view, as the positions of the new vertices introduced at a non-manifold vertex are identical to the input non-manifold vertex.

### Manifoldness Repair Example

In the following example, a non-manifold configuration is artificially created and fixed with the help of the functions described above.

File Polygon_mesh_processing/manifoldness_repair_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/boost/graph/iterator.h>
#include <iostream>
#include <iterator>
#include <string>
#include <vector>
namespace PMP = CGAL::Polygon_mesh_processing;
namespace NP = CGAL::parameters;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
void merge_vertices(vertex_descriptor v_keep, vertex_descriptor v_rm, Mesh& mesh)
{
std::cout << "merging vertices " << v_keep << " and " << v_rm << std::endl;
for(halfedge_descriptor h : CGAL::halfedges_around_target(v_rm, mesh))
set_target(h, v_keep, mesh); // to ensure that no halfedge points at the deleted vertex
remove_vertex(v_rm, mesh);
}
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh) || CGAL::is_empty(mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
// Artificially create non-manifoldness for the sake of the example by merging some vertices
vertex_descriptor v0 = *(vertices(mesh).begin());
vertex_descriptor v1 = *(--(vertices(mesh).end()));
merge_vertices(v0, v1, mesh);
// Count non manifold vertices
int counter = 0;
for(vertex_descriptor v : vertices(mesh))
{
if(PMP::is_non_manifold_vertex(v, mesh))
{
std::cout << "vertex " << v << " is non-manifold" << std::endl;
++counter;
}
}
std::cout << counter << " non-manifold occurrence(s)" << std::endl;
// Fix manifoldness by splitting non-manifold vertices
std::vector<std::vector<vertex_descriptor> > duplicated_vertices;
std::size_t new_vertices_nb = PMP::duplicate_non_manifold_vertices(mesh,
NP::output_iterator(
std::back_inserter(duplicated_vertices)));
std::cout << new_vertices_nb << " vertices have been added to fix mesh manifoldness" << std::endl;
for(std::size_t i=0; i<duplicated_vertices.size(); ++i)
{
std::cout << "Non-manifold vertex " << duplicated_vertices[i].front() << " was fixed by creating";
for(std::size_t j=1; j<duplicated_vertices[i].size(); ++j)
std::cout << " " << duplicated_vertices[i][j];
std::cout << std::endl;
}
return EXIT_SUCCESS;
}

## Duplicated Vertices in Boundary Cycles

Similarly to the problematic configuration described in the previous section, another issue that can be present in a polygon mesh is the occurrence of a "pinched" hole, that is the configuration where, when starting from a border halfedge and walking the halfedges of this border, a geometric position appears more than once (although, with different vertices) before reaching the initial border halfedge again. The functions `CGAL::Polygon_mesh_processing::merge_duplicated_vertices_in_boundary_cycle()` and `CGAL::Polygon_mesh_processing::merge_duplicated_vertices_in_boundary_cycle()`, which merge vertices at identical positions, can be used to repair this configuration.

# Geometric Repair

## Self-intersection Resolution (Autorefinement) in Triangle Soups

Given a soup of triangles, a self-intersection is defined as the intersection of two triangles from the soup such that the intersection is not defined by the convex hull of one, two or three shared vertices. In other words, it is an intersection that happens in the interior of one of the two triangles, or in the interior of one their edges, except if identical points are associated to different vertices of the triangle soup which would then also includes overlaps of duplicated points.

The function `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` provides a way to refine a triangle soup using the intersections of the triangles from the soup. In particular, if some points are duplicated they will be merged. Note that if a kernel with exact predicates but inexact constructions is used, some new self-intersections might be introduced due to the rounding of the coordinates of intersection points. The `apply_iterative_snap_rounding` option can be used to resolve this issue. When set to `true`, it ensures that the coordinates are rounded to fit in `double` with potential additional subdivisions, preventing any self-intersections from occurring.

## Removal of Almost Degenerate Triangle Faces

Triangle faces of a mesh made up of almost collinear points are badly shaped elements that might not be desirable to have in a mesh. The function `CGAL::Polygon_mesh_processing::remove_almost_degenerate_faces()` enables removing such elements, with user-defined parameters to qualify what almost means (`cap_threshold` and `needle_threshold`). As some badly shaped elements are inevitable (the triangulation of a long cylinder with only vertices on the top and bottom circles for example), extra parameters can be passed to prevent the removal of such elements (`collapse_length_threshold` and `flip_triangle_height_threshold`).

# Computing Normals

This package provides methods to compute normals on the polygon mesh. The normal can either be computed for each single face, or estimated for each vertex, as the average of its incident face normals. These computations are performed with :

- `CGAL::Polygon_mesh_processing::compute_face_normal()`
- `CGAL::Polygon_mesh_processing::compute_vertex_normal()`

Furthermore, we provide functions to compute all the normals to faces, or to vertices, or to both :

- `CGAL::Polygon_mesh_processing::compute_face_normals()`
- `CGAL::Polygon_mesh_processing::compute_vertex_normals()`
- `CGAL::Polygon_mesh_processing::compute_normals()`.

Property maps are used to record the computed normals.

## Normals Computation Examples

Property maps are an API introduced in the boost library that allows to associate values to keys. In the following examples we associate a normal vector to each vertex and to each face.

### Normals Computation for a Surface Mesh

The following example illustrates how to compute the normals to faces and vertices and store them in property maps provided by the class `Surface_mesh`.

File Polygon_mesh_processing/compute_normals_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <iostream>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point;
typedef K::Vector_3 Vector;
typedef CGAL::Surface_mesh<Point> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
auto vnormals = mesh.add_property_map<vertex_descriptor, Vector>("v:normal", CGAL::NULL_VECTOR).first;
auto fnormals = mesh.add_property_map<face_descriptor, Vector>("f:normal", CGAL::NULL_VECTOR).first;
PMP::compute_normals(mesh, vnormals, fnormals);
std::cout << "Vertex normals :" << std::endl;
for(vertex_descriptor vd: vertices(mesh))
std::cout << vnormals[vd] << std::endl;
std::cout << "Face normals :" << std::endl;
for(face_descriptor fd: faces(mesh))
std::cout << fnormals[fd] << std::endl;
return 0;
}
CGAL::Surface_mesh::add_property_map
std::pair< Property_map< I, T >, bool > add_property_map(std::string name=std::string(), const T t=T())

### Normals Computation for a Polyhedron_3

The following example illustrates how to compute the normals to faces and vertices and store them in ordered or unordered maps as the class `Polyhedron_3` does not provide storage for the normals.

File Polygon_mesh_processing/compute_normals_example_Polyhedron.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
// #include <CGAL/Unique_hash_map.h>
// #include <boost/unordered_map.hpp>
#include <boost/property_map/property_map.hpp>
#include <iostream>
#include <map>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point;
typedef K::Vector_3 Vector;
typedef CGAL::Polyhedron_3<K> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
std::map<face_descriptor,Vector> fnormals;
std::map<vertex_descriptor,Vector> vnormals;
// Instead of std::map you may use std::unordered_map, boost::unordered_map
// or CGAL::Unique_hash_map
// CGAL::Unique_hash_map<face_descriptor,Vector> fnormals;
// boost::unordered_map<vertex_descriptor,Vector> vnormals;
PMP::compute_normals(mesh, boost::make_assoc_property_map(vnormals),
boost::make_assoc_property_map(fnormals));
std::cout << "Face normals :" << std::endl;
for(face_descriptor fd: faces(mesh))
std::cout << fnormals[fd] << std::endl;
std::cout << "Vertex normals :" << std::endl;
for(vertex_descriptor vd: vertices(mesh))
std::cout << vnormals[vd] << std::endl;
return 0;
}

# Computing Curvatures

This package provides methods to compute curvatures on polygonal meshes based on Interpolated Corrected Curvatures on Polyhedral Surfaces [11]. This includes mean curvature, Gaussian curvature, principal curvatures and directions. These can be computed on triangle meshes, quad meshes, and meshes with n-gon faces (for n-gons, the centroid must be inside the n-gon face). The algorithms used prove to work well in general. Also, on meshes with noise on vertex positions, they give accurate results, under the condition that the correct vertex normals are provided.

It is worth noting that the Principal Curvatures and Directions can also be estimated using the Estimation of Local Differential Properties of Point-Sampled Surfaces package, which estimates the local differential quantities of a surface at a point using a local polynomial fitting (fitting a d-jet). Unlike the Interpolated Corrected Curvatures, the Jet Fitting method discards topological information and thus can be used on point clouds as well.

## Brief Background

Surface curvatures are quantities that describe the local geometry of a surface. They are important in many geometry processing applications. As surfaces are 2-dimensional objects (embedded in 3D), they can bend in 2 independent directions. These directions are called principal directions, and the amount of bending in each direction is called the principal curvature: \( k_1 \) and \( k_2 \) (denoting max and min curvatures). Curvature is usually expressed as scalar quantities like the mean curvature \( H \) and the Gaussian curvature \( K \) which are defined in terms of the principal curvatures.

The algorithms are based on the two papers [12] and [11]. They introduce a new way to compute curvatures on polygonal meshes. The main idea in [12] is based on decoupling the normal information from the position information, which is useful for dealing with digital surfaces, or meshes with noise on vertex positions. [11] introduces some extensions to this framework, as it uses linear interpolation on the corrected normal vector field to derive new closed-form equations for the corrected curvature measures. These interpolated curvature measures are the first step for computing the curvatures. For a triangle \( \tau_{ijk} \), with vertices i, j, k:

\begin{align*} \mu^{(0)}(\tau_{ijk}) = &\frac{1}{2} \langle \bar{\mathbf{u}} \mid (\mathbf{x}_j - \mathbf{x}_i) \times (\mathbf{x}_k - \mathbf{x}_i) \rangle, \\ \mu^{(1)}(\tau_{ijk}) = &\frac{1}{2} \langle \bar{\mathbf{u}} \mid (\mathbf{u}_k - \mathbf{u}_j) \times \mathbf{x}_i + (\mathbf{u}_i - \mathbf{u}_k) \times \mathbf{x}_j + (\mathbf{u}_j - \mathbf{u}_i) \times \mathbf{x}_k \rangle, \\ \mu^{(2)}(\tau_{ijk}) = &\frac{1}{2} \langle \mathbf{u}_i \mid \mathbf{u}_j \times \mathbf{u}_k \rangle, \\ \mu^{\mathbf{X},\mathbf{Y}}(\tau_{ijk}) = & \frac{1}{2} \big\langle \bar{\mathbf{u}} \big| \langle \mathbf{Y} | \mathbf{u}_k -\mathbf{u}_i \rangle \mathbf{X} \times (\mathbf{x}_j - \mathbf{x}_i) \big\rangle -\frac{1}{2} \big\langle \bar{\mathbf{u}} \big| \langle \mathbf{Y} | \mathbf{u}_j -\mathbf{u}_i \rangle \mathbf{X} \times (\mathbf{x}_k - \mathbf{x}_i) \big\rangle, \end{align*}

where \( \langle \cdot \mid \cdot \rangle \) denotes the usual scalar product, \( \bar{\mathbf{u}}=\frac{1}{3}( \mathbf{u}_i + \mathbf{u}_j + \mathbf{u}_k )\).

The first measure \( \mu^{(0)} \) is the area measure of the triangle, and the measures \( \mu^{(1)} \) and \( \mu^{(2)} \) are the mean and Gaussian corrected curvature measures of the triangle. The last measure \( \mu^{\mathbf{X},\mathbf{Y}} \) is the anisotropic corrected curvature measure of the triangle. The anisotropic measure is later used to compute the principal curvatures and directions through an eigenvalue solver.

The interpolated curvature measures are then computed for each vertex \( v \) as the sum of the curvature measures of the faces in a ball around \( v \) weighted by the inclusion ratio of the triangle in the ball. This ball radius is an optional (named) parameter of the function. There are 3 cases for the ball radius passed value:

- A positive value is passed: it is naturally used as the radius of the ball.
- 0 is passed, a small epsilon (`average_edge_length * 1e-6`) is used (to account for the convergence of curvatures at infinitely small balls).
- It is not specified (or negative), the sum is instead computed over the incident faces of the vertex \( v \).

To get the final curvature value for a vertex \( v \), the respective interpolated curvature measure is divided by the interpolated area measure.

\[ \mu^{(k)}( B ) = \sum_{\tau : \text{triangle} } \mu^{(k)}( \tau ) \frac{\mathrm{Area}( \tau \cap B )}{\mathrm{Area}(\tau)}. \]

## API

The implementation is generic in terms of mesh data structure. It can be used on `Surface_mesh`, `Polyhedron_3` and other polygonal mesh structures based on the concept `FaceGraph`.

These computations are performed using (on all vertices of the mesh) `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures()` where function named parameters are used to select the curvatures (and possibly directions) to be computed. An overload function with the same name but taking a given vertex is also available in case the computation should be done only for that vertex.

## Results

First, Figure 75.19 illustrates various curvature measures on a triangular mesh.

(a) | (b)

Figure 75.19 Mean curvature, Gaussian curvature, minimal principal curvature direction and maximal principal curvature direction on a mesh (ball radius set to 0.04).

(a) | (b)

Figure 75.20 When changing the integration ball radius, we obtain a scale space of curvature measure that can be used to tackle possible noise in the input as illustrated in the second row (mean curvature only with fixed colormap ranges and ball radii in {0.02,0.03,0.04,0.05}).

## Performance

The implemented algorithms exhibit a linear complexity in the number of faces of the mesh. It is worth noting that we pre-computed the vertex normals and passed them as a named parameter to the function to better estimate the performance of the curvature computation. For the data reported in the following table, we used a machine with an Intel Core i7-8750H CPU @ 2.20GHz, 16GB of RAM, on Windows 11, 64 bits and compiled with Visual Studio 2019.

Ball Radius | Computation | Spot (6k faces) | Bunny (144K faces) | Stanford Dragon (871K faces) | Old Age or Winter (6M faces)
--- | --- | --- | --- | --- | ---
vertex 1-ring faces (default) | Mean Curvature | < 0.001 s | 0.019 s | 0.11 s | 2.68 s
Gaussian Curvature | < 0.001 s | 0.017 s | 0.10 s | 2.77 s
Principal Curvatures & Directions | 0.002 s | 0.044 s | 0.25 s | 3.98 s
All (optimized for shared computations) | 0.003 s | 0.049 s | 0.28 s | 4.52 s
r = 0.1 * avg_edge_length | Mean Curvature | 0.017 s | 0.401 s | 2.66 s | 22.29 s
Gaussian Curvature | 0.018 s | 0.406 s | 2.63 s | 21.61 s
Principal Curvatures & Directions | 0.019 s | 0.430 s | 2.85 s | 23.55 s
All (optimized for shared computations) | 0.017 s | 0.428 s | 2.89 s | 24.16 s
r = 0.5 * avg_edge_length | Mean Curvature | 0.024 s | 0.388 s | 3.18 s | 22.79 s
Gaussian Curvature | 0.024 s | 0.392 s | 3.21 s | 23.58 s
Principal Curvatures & Directions | 0.027 s | 0.428 s | 3.41 s | 24.44 s
All (optimized for shared computations) | 0.025 s | 0.417 s | 3.44 s | 23.93 s

Figure 75.21 Performance of the curvature computation on various meshes (in seconds). The first 4 rows show the performance of the default value for the ball radius, which is using the 1-ring of neighboring faces around each vertex, instead of actually approximating the inclusion ratio of the faces in a ball of certain radius. The other rows show a ball radius of 0.1 (and 0.5) scaled by the average edge length of the mesh. It is clear that using the 1-ring of faces is much faster, but it might not be as effective when used on a noisy input mesh.

Property Maps are used to record the computed curvatures as shown in examples. In the following examples, for each property map, we associate a curvature value to each vertex.

## Interpolated Corrected Curvatures on a Surface Mesh Example

The following example illustrates how to compute the curvatures on vertices and store them in the property maps provided by the class `Surface_mesh`.

File Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/property_map.h>
#include <boost/graph/graph_traits.hpp>
#include <iostream>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel;
typedef CGAL::Surface_mesh<Epic_kernel::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
int main(int argc, char* argv[])
{
Mesh smesh;
const std::string filename = (argc > 1) ?
argv[1] :
CGAL::data_file_path("meshes/sphere.off");
if (!CGAL::IO::read_polygon_mesh(filename, smesh))
{
std::cerr << "Invalid input file." << std::endl;
return EXIT_FAILURE;
}
// creating and tying surface mesh property maps for curvatures (with defaults = 0)
bool created = false;
Mesh::Property_map<vertex_descriptor, Epic_kernel::FT>
mean_curvature_map, Gaussian_curvature_map;
std::tie(mean_curvature_map, created) =
smesh.add_property_map<vertex_descriptor, Epic_kernel::FT>("v:mean_curvature_map", 0);
assert(created);
std::tie(Gaussian_curvature_map, created) =
smesh.add_property_map<vertex_descriptor, Epic_kernel::FT>("v:Gaussian_curvature_map", 0);
assert(created);
// we use a tuple of 2 scalar values and 2 vectors for principal curvatures and directions
Mesh::Property_map<vertex_descriptor, PMP::Principal_curvatures_and_directions<Epic_kernel>>
principal_curvatures_and_directions_map;
std::tie(principal_curvatures_and_directions_map, created) =
smesh.add_property_map<vertex_descriptor, PMP::Principal_curvatures_and_directions<Epic_kernel>>
("v:principal_curvatures_and_directions_map", { 0, 0,
Epic_kernel::Vector_3(0,0,0),
Epic_kernel::Vector_3(0,0,0) });
assert(created);
PMP::interpolated_corrected_curvatures(smesh,
CGAL::parameters::vertex_mean_curvature_map(mean_curvature_map)
.vertex_Gaussian_curvature_map(Gaussian_curvature_map)
.vertex_principal_curvatures_and_directions_map(principal_curvatures_and_directions_map)
// uncomment to use an expansion ball radius of 0.5 to estimate the curvatures
// .ball_radius(0.5)
);
for (vertex_descriptor v : vertices(smesh))
{
auto PC = principal_curvatures_and_directions_map[v];
std::cout << v.idx() << ": HC = " << mean_curvature_map[v]
<< ", GC = " << Gaussian_curvature_map[v] << "\n"
<< ", PC = [ " << PC.min_curvature << " , " << PC.max_curvature << " ]\n";
}
return 0;
}

## Interpolated Corrected Curvatures on a Polyhedron Example

The following example illustrates how to compute the curvatures on vertices and store them in dynamic property maps as the class `Polyhedron_3` does not provide storage for the curvatures.

File Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/property_map.h>
#include <boost/graph/graph_traits.hpp>
#include <iostream>
#include <string>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel;
typedef CGAL::Polyhedron_3<Epic_kernel> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
int main(int argc, char* argv[])
{
Mesh polyhedron;
const std::string filename = (argc > 1) ?
argv[1] :
CGAL::data_file_path("meshes/sphere.off");
if (!CGAL::IO::read_polygon_mesh(filename, polyhedron))
{
std::cerr << "Invalid input file." << std::endl;
return EXIT_FAILURE;
}
// define property map to store curvature value and directions
boost::property_map<Mesh, CGAL::dynamic_vertex_property_t<Epic_kernel::FT>>::type
mean_curvature_map = get(CGAL::dynamic_vertex_property_t<Epic_kernel::FT>(), polyhedron),
Gaussian_curvature_map = get(CGAL::dynamic_vertex_property_t<Epic_kernel::FT>(), polyhedron);
boost::property_map<Mesh, CGAL::dynamic_vertex_property_t<PMP::Principal_curvatures_and_directions<Epic_kernel>>>::type
principal_curvatures_and_directions_map =
get(CGAL::dynamic_vertex_property_t<PMP::Principal_curvatures_and_directions<Epic_kernel>>(), polyhedron);
PMP::interpolated_corrected_curvatures(polyhedron,
CGAL::parameters::vertex_mean_curvature_map(mean_curvature_map)
.vertex_Gaussian_curvature_map(Gaussian_curvature_map)
.vertex_principal_curvatures_and_directions_map(principal_curvatures_and_directions_map)
// uncomment to use an expansion ball radius of 0.5 to estimate the curvatures
// .ball_radius(0.5)
);
int i = 0;
for (vertex_descriptor v : vertices(polyhedron))
{
auto PC = get(principal_curvatures_and_directions_map, v);
std::cout << i << ": HC = " << get(mean_curvature_map, v)
<< ", GC = " << get(Gaussian_curvature_map, v) << "\n"
<< ", PC = [ " << PC.min_curvature << " , " << PC.max_curvature << " ]\n";
i++;
}
}
CGAL::dynamic_vertex_property_t

## Interpolated Corrected Curvatures on a Vertex Example

The following example illustrates how to compute the curvatures on a specific vertex.

File Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Surface_mesh.h>
#include <boost/graph/graph_traits.hpp>
#include <iostream>
#include <string>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel;
typedef CGAL::Surface_mesh<Epic_kernel::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
int main(int argc, char* argv[])
{
// instantiating and reading mesh
Mesh smesh;
const std::string filename = (argc > 1) ?
argv[1] :
CGAL::data_file_path("meshes/sphere.off");
if (!CGAL::IO::read_polygon_mesh(filename, smesh))
{
std::cerr << "Invalid input file." << std::endl;
return EXIT_FAILURE;
}
// loop over vertices and use vertex_descriptor to compute a curvature on one vertex
for (vertex_descriptor v : vertices(smesh))
{
double h, g;
PMP::Principal_curvatures_and_directions<Epic_kernel> p;
PMP::interpolated_corrected_curvatures(v,
smesh,
CGAL::parameters::vertex_mean_curvature(std::ref(h))
.vertex_Gaussian_curvature(std::ref(g))
.vertex_principal_curvatures_and_directions(std::ref(p)));
// we can also specify a ball radius for expansion and a user defined vertex normals map using
// named parameters. Refer to interpolated_corrected_curvatures_SM.cpp to see example usage.
std::cout << v.idx() << ": HC = " << h
<< ", GC = " << g << "\n"
<< ", PC = [ " << p.min_curvature << " , " << p.max_curvature << " ]\n";
}
return 0;
}

## Discrete Curvatures

The package also provides methods to compute the standard, non-interpolated discrete mean and Gaussian curvatures on triangle meshes, based on the work of Meyer et al. [15]. These curvatures are computed at each vertex of the mesh, and are based on the angles of the incident triangles. The functions are:

- `CGAL::Polygon_mesh_processing::discrete_mean_curvature()`
- `CGAL::Polygon_mesh_processing::discrete_mean_curvatures()`
- `CGAL::Polygon_mesh_processing::discrete_Gaussian_curvature()`
- `CGAL::Polygon_mesh_processing::discrete_Gaussian_curvatures()`

# Slicer

The `CGAL::Polygon_mesh_slicer` is an operator that intersects a triangle surface mesh with a plane. It records the intersection as a set of polylines since the intersection can be made of more than one connected component. The degenerate case where the intersection is a single point is handled.

Figure 75.22 shows the polylines returned by the slicing operation for a triangle mesh and a set of parallel planes.

Figure 75.22 Slicing a mesh. A triangle mesh (left) and the polylines computed by the mesh slicer by intersecting a set of parallel planes (right).

## Slicer Example

The example below illustrates how to use the mesh slicer for a given triangle mesh and a plane. Two constructors are used in the example for pedagogical purposes.

File Polygon_mesh_processing/mesh_slicer_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_slicer.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/AABB_halfedge_graph_segment_primitive.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits_3.h>
#include <iostream>
#include <list>
#include <string>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef std::vector<K::Point_3> Polyline_type;
typedef std::list<Polyline_type> Polylines;
typedef CGAL::AABB_halfedge_graph_segment_primitive<Mesh> HGSP;
typedef CGAL::AABB_traits_3<K, HGSP> AABB_traits;
typedef CGAL::AABB_tree<AABB_traits> AABB_tree;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh) || CGAL::is_empty(mesh) || !CGAL::is_triangle_mesh(mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
// Slicer constructor from the mesh
CGAL::Polygon_mesh_slicer<Mesh, K> slicer(mesh);
Polylines polylines;
slicer(K::Plane_3(0, 0, 1, -0.4), std::back_inserter(polylines));
std::cout << "At z = 0.4, the slicer intersects "
<< polylines.size() << " polylines" << std::endl;
polylines.clear();
slicer(K::Plane_3(0, 0, 1, 0.2), std::back_inserter(polylines));
std::cout << "At z = -0.2, the slicer intersects "
<< polylines.size() << " polylines" << std::endl;
polylines.clear();
// Use the Slicer constructor from a pre-built AABB_tree
AABB_tree tree(edges(mesh).first, edges(mesh).second, mesh);
CGAL::Polygon_mesh_slicer<Mesh, K> slicer_aabb(mesh, tree);
slicer_aabb(K::Plane_3(0, 0, 1, -0.4), std::back_inserter(polylines));
std::cout << "At z = 0.4, the slicer intersects "
<< polylines.size() << " polylines" << std::endl;
polylines.clear();
return 0;
}
CGAL::AABB_halfedge_graph_segment_primitive
CGAL::AABB_traits_3
CGAL::AABB_tree
CGAL::Polygon_mesh_slicer
Function object that computes the intersection of a plane with a triangulated surface mesh.
Definition: Polygon_mesh_slicer.h:89
CGAL::AABB_traits
AABB_traits_3< GeomTraits, AABBPrimitive, BboxMap > AABB_traits

# Connected Components

This package provides functions to enumerate and store the connected components of a polygon mesh. The connected components can be either closed and geometrically separated, or separated by border or user-specified constraint edges.

First, the function `CGAL::Polygon_mesh_processing::connected_component()` collects all the faces that belong to the same connected component as the face that is given as a parameter.

Then, `CGAL::Polygon_mesh_processing::connected_components()` collects all the connected components, and fills a property map with the indices of the different connected components.

The functions `CGAL::Polygon_mesh_processing::keep_connected_components()` and `CGAL::Polygon_mesh_processing::remove_connected_components()` enable the user to keep and remove only a selection of connected components, provided either as a range of faces that belong to the desired connected components or as a range of connected component ids (one or more per connected component).

When a triangle mesh has no boundary, it partitions the 3D space in different volumes. The function `CGAL::Polygon_mesh_processing::volume_connected_components()` can be used to assign to each face an id per volume defined by the surface connected components.

Finally, it can be useful to quickly remove some connected components, for example for noisy data where small connected components should be discarded in favor of major connected components. The function `CGAL::Polygon_mesh_processing::keep_largest_connected_components()` enables the user to keep only a given number from the largest connected components. The size of a connected component is given by the sum of the sizes of the faces it contains; by default, the size of a face is `1` and thus the size of a connected component is equal to the number of faces it contains. However, it is also possible to pass a face size map, such that the size of the face is its area, for example. Similarly to the previous function, the function `CGAL::Polygon_mesh_processing::keep_large_connected_components()` can be used to discard all connected components whose size is below a user-defined threshold.

Also, the function `CGAL::Polygon_mesh_processing::split_connected_components()` enables the user to split the connected components of a polygon mesh in as many polygon meshes.

## Connected Components Example

The first example shows how to record the connected components of a polygon mesh. In particular, we provide an example for the optional parameter `EdgeConstraintMap`, a property map that returns information about an edge being a constraint or not. A constraint provides a mean to demarcate the border of a connected component, and prevents the propagation of a connected component index to cross it.

File Polygon_mesh_processing/connected_components_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/connected_components.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <boost/iterator/function_output_iterator.hpp>
#include <boost/property_map/property_map.hpp>
#include <iostream>
#include <iterator>
#include <map>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Compare_dihedral_angle_3 Compare_dihedral_angle_3;
typedef CGAL::Surface_mesh<Point> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
template <typename G>
struct Constraint
{
typedef typename boost::graph_traits<G>::edge_descriptor edge_descriptor;
typedef boost::readable_property_map_tag category;
typedef bool value_type;
typedef bool reference;
typedef edge_descriptor key_type;
Constraint()
:g_(NULL)
{}
Constraint(G& g, double bound)
: g_(&g), bound_(bound)
{}
value_type operator[](edge_descriptor e) const
{
const G& g = *g_;
return compare_(g.point(source(e, g)),
g.point(target(e, g)),
g.point(target(next(halfedge(e, g), g), g)),
g.point(target(next(opposite(halfedge(e, g), g), g), g)),
bound_) == CGAL::SMALLER;
}
friend inline
value_type get(const Constraint& m, const key_type k)
{
return m[k];
}
const G* g_;
Compare_dihedral_angle_3 compare_;
double bound_;
};
template <typename PM>
struct Put_true
{
Put_true(const PM pm)
:pm(pm)
{}
template <typename T>
void operator()(const T& t)
{
put(pm, t, true);
}
PM pm;
};
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby_3cc.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
const double bound = std::cos(0.75 * CGAL_PI);
std::vector<face_descriptor> cc;
face_descriptor fd = *faces(mesh).first;
PMP::connected_component(fd,
mesh,
std::back_inserter(cc));
std::cerr << "Connected components without edge constraints" << std::endl;
std::cerr << cc.size() << " faces in the CC of " << fd << std::endl;
// Instead of writing the faces into a container, you can set a face property to true
typedef Mesh::Property_map<face_descriptor, bool> F_select_map;
F_select_map fselect_map =
mesh.add_property_map<face_descriptor, bool>("f:select", false).first;
PMP::connected_component(fd,
mesh,
boost::make_function_output_iterator(Put_true<F_select_map>(fselect_map)));
std::cerr << "\nConnected components with edge constraints (dihedral angle < 3/4 pi)" << std::endl;
Mesh::Property_map<face_descriptor, std::size_t> fccmap =
mesh.add_property_map<face_descriptor, std::size_t>("f:CC").first;
std::size_t num = PMP::connected_components(mesh,
fccmap,
CGAL::parameters::edge_is_constrained_map(Constraint<Mesh>(mesh, bound)));
std::cerr << "- The graph has " << num << " connected components (face connectivity)" << std::endl;
typedef std::map<std::size_t/*index of CC*/, unsigned int/*nb*/> Components_size;
Components_size nb_per_cc;
for(face_descriptor f : faces(mesh)){
nb_per_cc[ fccmap[f] ]++;
}
for(const Components_size::value_type& cc : nb_per_cc){
std::cout << "\t CC #" << cc.first
<< " is made of " << cc.second << " faces" << std::endl;
}
std::cerr << "- We keep only components which have at least 4 faces" << std::endl;
PMP::keep_large_connected_components(mesh,
4,
CGAL::parameters::edge_is_constrained_map(Constraint<Mesh>(mesh, bound)));
std::cerr << "- We keep the two largest components" << std::endl;
PMP::keep_largest_connected_components(mesh,
2,
CGAL::parameters::edge_is_constrained_map(Constraint<Mesh>(mesh, bound)));
return 0;
}
CGAL::opposite
Oriented_side opposite(const Oriented_side &o)
CGAL::SMALLER
SMALLER

The second example shows how to use the class template `Face_filtered_graph` which enables to treat one or several connected components as a face graph.

File Polygon_mesh_processing/face_filtered_graph_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/connected_components.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/boost/graph/Face_filtered_graph.h>
#include <boost/property_map/property_map.hpp>
#include <iostream>
#include <map>
#include <string>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Mesh;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
typedef boost::graph_traits<Mesh>::faces_size_type faces_size_type;
typedef Mesh::Property_map<face_descriptor, faces_size_type> FCCmap;
typedef CGAL::Face_filtered_graph<Mesh> Filtered_graph;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby_3cc.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
FCCmap fccmap = mesh.add_property_map<face_descriptor, faces_size_type>("f:CC").first;
faces_size_type num = PMP::connected_components(mesh,fccmap);
std::cerr << "- The graph has " << num << " connected components (face connectivity)" << std::endl;
std::cout << "The faces in component 0 are:" << std::endl;
Filtered_graph ffg(mesh, 0, fccmap);
for(boost::graph_traits<Filtered_graph>::face_descriptor f : faces(ffg))
std::cout << f << std::endl;
if(num > 1)
{
std::vector<faces_size_type> components;
components.push_back(0);
components.push_back(1);
std::cout << "The faces in components 0 and 1 are:" << std::endl;
ffg.set_selected_faces(components, fccmap);
for(Filtered_graph::face_descriptor f : faces(ffg))
std::cout << f << std::endl;
}
return 0;
}
CGAL::Face_filtered_graph

# Hausdorff Distance

This package provides methods to compute (approximate) distances between meshes and point sets.

## Approximate Hausdorff Distance

The function `approximate_Hausdorff_distance()` computes an approximation of the Hausdorff distance from a mesh `tm1` to a mesh `tm2`. Given a a sampling of `tm1`, it computes the distance to `tm2` of the farthest sample point to `tm2` [7]. The symmetric version (`approximate_symmetric_Hausdorff_distance()`) is the maximum of the two non-symmetric distances. Internally, points are sampled using `sample_triangle_mesh()` and the distance to each sample point is computed using `max_distance_to_triangle_mesh()`. The quality of the approximation depends on the quality of the sampling and the runtime depends on the number of sample points. Three sampling methods with different parameters are provided (see Figure 75.23).

Figure 75.23 Sampling of a triangle mesh using different sampling methods. From left to right: (a) Grid sampling, (b) Monte-Carlo sampling with fixed number of points per face and per edge, (c) Monte-Carlo sampling with a number of points proportional to the area/length, and (d) Uniform random sampling. The four pictures represent the sampling on the same portion of a mesh, parameters were adjusted so that the total number of points sampled in faces (blue points) and on edges (red points) are roughly the same. Note that when using the random uniform sampling some faces/edges may not contain any point, but this method is the only one that allows to exactly match a given number of points.

The function `approximate_max_distance_to_point_set()` computes an approximation of the Hausdorff distance from a mesh to a point set. For each triangle, a lower and upper bound of the Hausdorff distance to the point set are computed. Triangles are refined until the difference between the bounds is lower than a user-defined precision threshold.

### Approximate Hausdorff Distance Example

In the following example, a mesh is isotropically remeshed and the approximate distance between the input and the output is computed.

File Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/distance.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#define TAG CGAL::Parallel_if_available_tag
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int, char**)
{
Mesh tm1, tm2;
CGAL::make_tetrahedron(Point(.0,.0,.0),
Point(2,.0,.0),
Point(1,1,1),
Point(1,.0,2),
tm1);
tm2 = tm1;
CGAL::Polygon_mesh_processing::isotropic_remeshing(tm2.faces(),.05, tm2);
std::cout << "Approximated Hausdorff distance: "
<< CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance
<TAG>(tm1, tm2, CGAL::parameters::number_of_points_per_area_unit(4000))
<< std::endl;
return 0;
}
CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance
double approximate_Hausdorff_distance(const TriangleMesh &tm1, const TriangleMesh &tm2, const NamedParameters1 &np1=parameters::default_values(), const NamedParameters2 &np2=parameters::default_values())
computes the approximate Hausdorff distance from tm1 to tm2 by returning the distance of the farthest...
Definition: distance.h:1173
CGAL::Polygon_mesh_processing::isotropic_remeshing
void isotropic_remeshing(const FaceRange &faces, SizingFunction &sizing, PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
remeshes a triangulated region of a polygon mesh.
Definition: remesh.h:214
CGAL::make_tetrahedron
boost::graph_traits< Graph >::halfedge_descriptor make_tetrahedron(const P &p0, const P &p1, const P &p2, const P &p3, Graph &g)

### Max Distance Between Point Set and Surface Example

In Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp, a triangulated surface mesh is constructed from a point set using the Poisson reconstruction algorithm , and the distance between the point set and the reconstructed surface is computed with the following code:

// computes the approximation error of the reconstruction
double max_dist =
CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set
(output_mesh,
CGAL::make_range (boost::make_transform_iterator
(points.begin(), CGAL::Property_map_to_unary_function<Point_map>()),
boost::make_transform_iterator
(points.end(), CGAL::Property_map_to_unary_function<Point_map>())),
4000);
std::cout << "Max distance to point_set: " << max_dist << std::endl;
CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set
double approximate_max_distance_to_point_set(const TriangleMesh &tm, const PointRange &points, const double precision, const NamedParameters &np=parameters::default_values())
returns an approximation of the distance between points and the point lying on tm that is the farthes...
Definition: distance.h:1252
CGAL::Property_map_to_unary_function

## Bounded Hausdorff Distance

The function `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` computes an estimate of the Hausdorff distance of two triangle meshes which is bounded by a user-given error bound. Given two meshes `tm1` and `tm2`, it follows the procedure given by [17]. Namely, a bounded volume hierarchy (BVH) is built on `tm1` and `tm2` respectively. The BVH on `tm1` is used to iterate over all triangles in `tm1`. Throughout the traversal, the procedure keeps track of a global lower and upper bound on the Hausdorff distance respectively. For each triangle `t` in `tm1`, by traversing the BVH on `tm2`, it is estimated via the global bounds whether `t` can still contribute to the actual Hausdorff distance. From this process, a set of candidate triangles is selected.

The candidate triangles are subsequently subdivided and for each smaller triangle, the BVH on `tm2` is traversed again. This is repeated until the triangle is smaller than the user-given error bound, all vertices of the triangle are projected onto the same triangle in `tm2`, or the triangle's upper bound is lower than the global lower bound. After creation, the subdivided triangles are added to the list of candidate triangles. Thereby, all candidate triangles are processed until a triangle is found in which the Hausdorff distance is realized or in which it is guaranteed to be realized within the user-given error bound.

In the current implementation, the BVH used is an AABB-tree and not the swept sphere volumes as used in the original implementation. This should explain the runtime difference observed with the original implementation.

The function `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` computes the one-sided Hausdorff distance from `tm1` to `tm2`. This component also provides the symmetric distance `CGAL::Polygon_mesh_processing::bounded_error_symmetric_Hausdorff_distance()` and an utility function called `CGAL::Polygon_mesh_processing::is_Hausdorff_distance_larger()` that returns `true` if the Hausdorff distance between two meshes is larger than the user-defined max distance.

### Bounded Hausdorff Distance Example

In the following examples: (a) the distance of a tetrahedron to a remeshed version of itself is computed, (b) the distance of two geometries is computed which is realized strictly in the interior of a triangle of the first geometry, (c) a perturbation of a user-given mesh is compared to the original user-given mesh, (d) two user-given meshes are compared, where the second mesh is gradually moved away from the first one.

File Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/OFF.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Polygon_mesh_processing/distance.h>
#include <CGAL/Polygon_mesh_processing/transform.h>
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = typename Kernel::FT;
using Point_3 = typename Kernel::Point_3;
using Vector_3 = typename Kernel::Vector_3;
using Affine_transformation_3 = typename Kernel::Aff_transformation_3;
using TAG = CGAL::Sequential_tag;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
using Polyhedron = CGAL::Polyhedron_3<Kernel>;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char** argv) {
const double error_bound = 1e-4;
const std::string filepath = (argc > 1 ? argv[1] : CGAL::data_file_path("meshes/blobby.off"));
// We create a tetrahedron, remesh it, and compute the distance.
// The expected distance is error_bound.
std::cout << std::endl << "* remeshing tetrahedron example:" << std::endl;
Surface_mesh mesh1, mesh2;
CGAL::make_tetrahedron(
Point_3(0, 0, 0), Point_3(2, 0, 0),
Point_3(1, 1, 1), Point_3(1, 0, 2), mesh1);
mesh2 = mesh1;
using edge_descriptor = typename boost::graph_traits<Surface_mesh>::edge_descriptor;
Surface_mesh::Property_map<edge_descriptor, bool> is_constrained_map =
mesh2.add_property_map<edge_descriptor, bool>("e:is_constrained", true).first;
const double target_edge_length = 0.05;
PMP::isotropic_remeshing(
mesh2.faces(), target_edge_length, mesh2,
CGAL::parameters::edge_is_constrained_map(is_constrained_map));
std::cout << "* one-sided bounded-error Hausdorff distance: " <<
PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound) << std::endl;
// We load a mesh, save it in two different containers, and
// translate the second mesh by 1 unit. The expected distance is 1.
std::cout << std::endl << "* moving mesh example:" << std::endl;
Surface_mesh surface_mesh;
CGAL::IO::read_OFF(filepath, surface_mesh);
Polyhedron polyhedron;
CGAL::IO::read_OFF(filepath, polyhedron);
PMP::transform(Affine_transformation_3(CGAL::Translation(),
Vector_3(FT(0), FT(0), FT(1))), polyhedron);
std::cout << "* symmetric bounded-error Hausdorff distance: " <<
PMP::bounded_error_symmetric_Hausdorff_distance<TAG>(surface_mesh, polyhedron, error_bound)
<< std::endl << std::endl;
return EXIT_SUCCESS;
}
CGAL::Translation
CGAL::Sequential_tag

# Feature Detection

This package provides methods to detect some features of a polygon mesh.

The function `CGAL::Polygon_mesh_processing::sharp_edges_segmentation()` detects the sharp edges of a polygon mesh and deduces surface patches and vertices incidences. It can be split into three functions : `CGAL::Polygon_mesh_processing::detect_sharp_edges()`, `CGAL::Polygon_mesh_processing::connected_components()` and `CGAL::Polygon_mesh_processing::detect_vertex_incident_patches()`, that respectively detect the sharp edges, compute the patch indices, and give each of `pmesh` vertices the patch indices of its incident faces.

## Feature Detection Example

In the following example, we count how many edges of `pmesh` are incident to two faces whose normals form an angle smaller than 90 degrees, and the number of surface patches that are separated by these edges.

File Polygon_mesh_processing/detect_features_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/detect_features.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <iostream>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/P.off");
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
typedef boost::property_map<Mesh, CGAL::edge_is_feature_t>::type EIFMap;
typedef boost::property_map<Mesh, CGAL::face_patch_id_t<int> >::type PIMap;
typedef boost::property_map<Mesh, CGAL::vertex_incident_patches_t<int> >::type VIMap;
EIFMap eif = get(CGAL::edge_is_feature, mesh);
PIMap pid = get(CGAL::face_patch_id_t<int>(), mesh);
VIMap vip = get(CGAL::vertex_incident_patches_t<int>(), mesh);
std::size_t number_of_patches
= PMP::sharp_edges_segmentation(mesh, 90, eif, pid,
CGAL::parameters::vertex_incident_patches_map(vip));
std::size_t nb_sharp_edges = 0;
for(boost::graph_traits<Mesh>::edge_descriptor e : edges(mesh))
{
if(get(eif, e))
++nb_sharp_edges;
}
std::cout << "This mesh contains " << nb_sharp_edges << " sharp edges" << std::endl;
std::cout << " and " << number_of_patches << " surface patches." << std::endl;
return 0;
}

# Implementation History

A first version of this package was started by Ilker O. Yaz and Sébastien Loriot. Jane Tournois worked on the finalization of the API, code, and documentation.

A prototype of mesh and shape smoothing was developed during the 2017 edition of the Google Summer of Code by Konstantinos Katrioplas, under supervision of Jane Tournois. It was finalized by Mael Rouxel-Labbé and integrated in CGAL 5.0.

Functionalities related to mesh and polygon soup repair have been introduced steadily over multiple versions since CGAL 4.10, in joint work between Sébastien Loriot and Mael Rouxel-Labbé.

The polyhedral envelope containment check was integrated in CGAL 5.3. The implementation makes use of the version of https://github.com/wangbolun300/fast-envelope available on 7th of October 2020. It only uses the high level algorithm of checking that a query is covered by a set of prisms, where each prism is an offset for an input triangle. That is, the implementation in CGAL does not use indirect predicates.

The interpolated corrected curvatures were implemented during GSoC 2022. This was implemented by Hossam Saeed and under the supervision of David Coeurjolly, Jaques-Olivier Lachaud, and Sébastien Loriot. The implementation is based on [11]. DGtal's implementation was also used as a reference during the project.

The curvature-based sizing field version of isotropic remeshing was added by Ivan Pađen during GSoC 2023, under the supervision of Sébastien Loriot and Jane Tournois.

The implementation of the ACVD Remeshing method was initiated by Hossam Saeed during GSoC 2023, under supervision of Sébastien Valette and Sébastien Loriot, who later finalized the work. The implementation is based on [19]. and preceding work. ACVD's implementation was also used as a reference during the project. The `apply_iterative_snap_rounding` option for autorefinement was implemented in 2025, by Léo Valque, based on his work with Sylvain Lazard [13].