\( \newcommand{\E}{\mathrm{E}} \) \( \newcommand{\A}{\mathrm{A}} \) \( \newcommand{\R}{\mathrm{R}} \) \( \newcommand{\N}{\mathrm{N}} \) \( \newcommand{\Q}{\mathrm{Q}} \) \( \newcommand{\Z}{\mathrm{Z}} \) \( \def\ccSum #1#2#3{ \sum_{#1}^{#2}{#3} } \def\ccProd #1#2#3{ \sum_{#1}^{#2}{#3} }\)
CGAL 4.12.1 - Polygon Mesh Processing
User Manual

Authors
Sébastien Loriot, Jane Tournois, Ilker O. Yaz
neptun_head.jpg

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 [4].

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 for Polygon Mesh Processing describes their usage and provides a list of the parameters that are used in this package.

Outline

The algorithms described in this manual are organized in sections:

  • Meshing : meshing algorithms, including triangulation of non-triangulated meshes, refinement, optimization by fairing, and isotropic remeshing of triangulated surface meshes.
  • 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 Repairing : reparation of polygon meshes and polygon soups.
  • Computing Normals : normal computation at vertices and on faces of 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, ...)

Meshing

A surface patch can be refined by inserting new vertices and flipping edges to get a triangulation. Using a criterion presented in [6], 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 [3]. The visual results of aforementioned steps are depicted by Figure 60.5 (c and d).

API

Meshing

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

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

The incremental triangle-based isotropic remeshing algorithm introduced by Botsch et al [2], [4] 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 60.1. The algorithm has only two parameters : the target edge length for the remeshed surface patch, and the number of iterations of the abovementioned sequence of operations. The bigger this number, the smoother and closer to target edge length the mesh will be.

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.

iso_remeshing.png
Figure 60.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.

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 <fstream>
#include <map>
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef Polyhedron::Vertex_handle Vertex_handle;
// 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++];
Polyhedron::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 char* filename = (argc > 1) ? argv[1] : "data/blobby.off";
std::ifstream input(filename);
Polyhedron poly;
if ( !input || !(input >> poly) || poly.empty()
std::cerr << "Not a valid input file." << std::endl;
return 1;
}
std::vector<Polyhedron::Facet_handle> new_facets;
std::vector<Vertex_handle> new_vertices;
faces(poly),
std::back_inserter(new_facets),
std::back_inserter(new_vertices),
CGAL::Polygon_mesh_processing::parameters::density_control_factor(2.));
std::ofstream refined_off("refined.off");
refined_off << poly;
refined_off.close();
std::cout << "Refinement added " << new_vertices.size() << " vertices." << std::endl;
Polyhedron::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 = CGAL::Polygon_mesh_processing::fair(poly, region);
std::cout << "Fairing : " << (success ? "succeeded" : "failed") << std::endl;
std::ofstream faired_off("faired.off");
faired_off << poly;
faired_off.close();
return 0;
}

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 <boost/foreach.hpp>
#include <fstream>
typedef Kernel::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Surface_mesh;
int main(int argc, char* argv[])
{
const char* filename = (argc > 1) ? argv[1] : "data/P.off";
const char* outfilename = (argc > 2) ? argv[2] : "P_tri.off";
std::ifstream input(filename);
Surface_mesh mesh;
if (!input || !(input >> mesh) || mesh.is_empty())
{
std::cerr << "Not a valid off file." << std::endl;
return 1;
}
// Confirm that all faces are triangles.
BOOST_FOREACH(boost::graph_traits<Surface_mesh>::face_descriptor fit, faces(mesh))
if (next(next(halfedge(fit, mesh), mesh), mesh)
!= prev(halfedge(fit, mesh), mesh))
std::cerr << "Error: non-triangular face left in mesh." << std::endl;
std::ofstream cube_off(outfilename);
cube_off << mesh;
return 0;
}

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 <boost/function_output_iterator.hpp>
#include <fstream>
#include <vector>
typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
typedef boost::graph_traits<Mesh>::edge_descriptor edge_descriptor;
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 char* filename = (argc > 1) ? argv[1] : "data/pig.off";
std::ifstream input(filename);
Mesh mesh;
if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh)) {
std::cerr << "Not a valid input file." << std::endl;
return 1;
}
double target_edge_length = 0.04;
unsigned int nb_iter = 3;
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;
faces(mesh),
target_edge_length,
mesh,
PMP::parameters::number_of_iterations(nb_iter)
.protect_constraints(true)//i.e. protect border, here
);
std::cout << "Remeshing done." << std::endl;
return 0;
}

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.

corefine.png
Figure 60.2 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.

bounded_vols.jpg
Figure 60.3 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

bool_op.png
Figure 60.4 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. Note that 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.

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 <fstream>
int main(int argc, char* argv[])
{
const char* filename1 = (argc > 1) ? argv[1] : "data/blobby.off";
const char* filename2 = (argc > 2) ? argv[2] : "data/eight.off";
std::ifstream input(filename1);
Mesh mesh1, mesh2;
if (!input || !(input >> mesh1))
{
std::cerr << "First mesh is not a valid off file." << std::endl;
return 1;
}
input.close();
input.open(filename2);
if (!input || !(input >> mesh2))
{
std::cerr << "Second mesh is not a valid off file." << 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";
std::ofstream output("union.off");
output << out;
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 substract 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/remesh.h>
#include <CGAL/boost/graph/selection.h>
#include <fstream>
typedef boost::graph_traits<Mesh>::edge_descriptor edge_descriptor;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
namespace params = PMP::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 char* filename1 = (argc > 1) ? argv[1] : "data/blobby.off";
const char* filename2 = (argc > 2) ? argv[2] : "data/eight.off";
std::ifstream input(filename1);
Mesh mesh1, mesh2;
if (!input || !(input >> mesh1))
{
std::cerr << "First mesh is not a valid off file." << std::endl;
return 1;
}
input.close();
input.open(filename2);
if (!input || !(input >> mesh2))
{
std::cerr << "Second mesh is not a valid off file." << 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 =
mesh2,
mesh1,
params::all_default(), // default parameters for mesh1
params::all_default(), // default parameters for mesh2
params::edge_is_constrained_map(is_constrained_map));
if (valid_difference)
{
std::cout << "Difference was successfully computed\n";
std::ofstream output("difference.off");
output << mesh1;
}
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);
BOOST_FOREACH(edge_descriptor e, edges(mesh1))
if (is_constrained_map[e])
{
// insert all faces incident to the target vertex
BOOST_FOREACH(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
selected_faces,
0.02,
mesh1,
params::edge_is_constrained_map(is_constrained_map) );
std::ofstream output("difference_remeshed.off");
output << mesh1;
return 0;
}

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 <fstream>
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef Mesh::Property_map<vertex_descriptor,EK::Point_3> Exact_point_map;
typedef Mesh::Property_map<vertex_descriptor,bool> Exact_point_computed;
namespace params = PMP::parameters;
struct Coref_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>::category category;
typedef boost::property_traits<Exact_point_map>::key_type key_type;
// exterior references
Exact_point_computed* exact_point_computed_ptr;
Exact_point_map* exact_point_ptr;
Mesh* mesh_ptr;
Exact_point_computed& exact_point_computed() const
{
CGAL_assertion(exact_point_computed_ptr!=NULL);
return *exact_point_computed_ptr;
}
Exact_point_map& exact_point() const
{
CGAL_assertion(exact_point_ptr!=NULL);
return *exact_point_ptr;
}
Mesh& mesh() const
{
CGAL_assertion(mesh_ptr!=NULL);
return *mesh_ptr;
}
// Converters
Coref_point_map()
: exact_point_computed_ptr(NULL)
, exact_point_ptr(NULL)
, mesh_ptr(NULL)
{}
Coref_point_map(Exact_point_map& ep,
Exact_point_computed& epc,
Mesh& m)
: exact_point_computed_ptr(&epc)
, exact_point_ptr(&ep)
, mesh_ptr(&m)
{}
friend
reference get(const Coref_point_map& map, key_type k)
{
// create exact point if it does not exist
if (!map.exact_point_computed()[k]){
map.exact_point()[k]=map.to_exact(map.mesh().point(k));
map.exact_point_computed()[k]=true;
}
return map.exact_point()[k];
}
friend
void put(const Coref_point_map& map, key_type k, const EK::Point_3& p)
{
map.exact_point_computed()[k]=true;
map.exact_point()[k]=p;
// create the input point from the exact one
map.mesh().point(k)=map.to_input(p);
}
};
int main(int argc, char* argv[])
{
const char* filename1 = (argc > 1) ? argv[1] : "data/blobby.off";
const char* filename2 = (argc > 2) ? argv[2] : "data/eight.off";
std::ifstream input(filename1);
Mesh mesh1, mesh2;
if (!input || !(input >> mesh1))
{
std::cerr << "First mesh is not a valid off file." << std::endl;
return 1;
}
input.close();
input.open(filename2);
if (!input || !(input >> mesh2))
{
std::cerr << "Second mesh is not a valid off file." << std::endl;
return 1;
}
Exact_point_map mesh1_exact_points =
mesh1.add_property_map<vertex_descriptor,EK::Point_3>("e:exact_point").first;
Exact_point_computed mesh1_exact_points_computed =
mesh1.add_property_map<vertex_descriptor,bool>("e:exact_points_computed").first;
Exact_point_map mesh2_exact_points =
mesh2.add_property_map<vertex_descriptor,EK::Point_3>("e:exact_point").first;
Exact_point_computed mesh2_exact_points_computed =
mesh2.add_property_map<vertex_descriptor,bool>("e:exact_points_computed").first;
Coref_point_map mesh1_pm(mesh1_exact_points, mesh1_exact_points_computed, mesh1);
Coref_point_map mesh2_pm(mesh2_exact_points, mesh2_exact_points_computed, mesh2);
mesh2,
mesh1,
params::vertex_point_map(mesh1_pm),
params::vertex_point_map(mesh2_pm),
params::vertex_point_map(mesh1_pm) ) )
{
mesh2,
mesh2,
params::vertex_point_map(mesh1_pm),
params::vertex_point_map(mesh2_pm),
params::vertex_point_map(mesh2_pm) ) )
{
std::cout << "Intersection and union were successfully computed\n";
std::ofstream output("inter_union.off");
output << mesh2;
return 0;
}
std::cout << "Union could not be computed\n";
return 1;
}
std::cout << "Intersection could not be computed\n";
return 1;
}

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 [6] 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 [7], 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.

mech_hole_horz.jpg
Figure 60.5 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:

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>
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
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;
polyline_collinear,
back_inserter(patch_will_be_empty));
CGAL_assertion(patch_will_be_empty.empty());
return 0;
}

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 to get a faired mesh with no hole.


File Polygon_mesh_processing/hole_filling_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <boost/foreach.hpp>
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef Polyhedron::Halfedge_handle Halfedge_handle;
typedef Polyhedron::Facet_handle Facet_handle;
typedef Polyhedron::Vertex_handle Vertex_handle;
int main(int argc, char* argv[])
{
const char* filename = (argc > 1) ? argv[1] : "data/mech-holes-shark.off";
std::ifstream input(filename);
Polyhedron poly;
if ( !input || !(input >> poly) || poly.empty() ) {
std::cerr << "Not a valid off file." << std::endl;
return 1;
}
// Incrementally fill the holes
unsigned int nb_holes = 0;
BOOST_FOREACH(Halfedge_handle h, halfedges(poly))
{
if(h->is_border())
{
std::vector<Facet_handle> patch_facets;
std::vector<Vertex_handle> patch_vertices;
bool success = CGAL::cpp11::get<0>(
poly,
h,
std::back_inserter(patch_facets),
std::back_inserter(patch_vertices),
CGAL::Polygon_mesh_processing::parameters::vertex_point_map(get(CGAL::vertex_point, poly)).
geom_traits(Kernel())) );
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 << " Fairing : " << (success ? "succeeded" : "failed") << std::endl;
++nb_holes;
}
}
std::cout << std::endl;
std::cout << nb_holes << " holes have been filled" << std::endl;
std::ofstream out("filled.off");
out.precision(17);
out << poly << std::endl;
return 0;
}

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

Performance

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

elephants-with-holes.png
Figure 60.7 The elephant on the left/right has a hole with 963/7657 vertices.

This takes time

# 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 packages provides several predicates to be evaluated with respect to a triangle mesh.

Self Intersections

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

selfintersections.jpg
Figure 60.8 Detecting self-intersections on a triangle mesh. The intersecting triangles are displayed in dark grey on the right image.

Self Intersections Example


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 <fstream>
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
int main(int argc, char* argv[])
{
const char* filename = (argc > 1) ? argv[1] : "data/pig.off";
std::ifstream input(filename);
Mesh mesh;
if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh))
{
std::cerr << "Not a valid input file." << std::endl;
return 1;
}
bool intersecting = PMP::does_self_intersect(mesh,
PMP::parameters::vertex_point_map(get(CGAL::vertex_point, mesh)));
std::cout
<< (intersecting ? "There are self-intersections." : "There is no self-intersection.")
<< std::endl;
std::vector<std::pair<face_descriptor, face_descriptor> > intersected_tris;
PMP::self_intersections(mesh, std::back_inserter(intersected_tris));
std::cout << intersected_tris.size() << " pairs of triangles intersect." << std::endl;
return 0;
}

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/point_generators_3.h>
#include <CGAL/Side_of_triangle_mesh.h>
#include <vector>
#include <fstream>
#include <limits>
#include <boost/foreach.hpp>
typedef K::Point_3 Point;
typedef CGAL::Polyhedron_3<K> Polyhedron;
double max_coordinate(const Polyhedron& poly)
{
double max_coord = -std::numeric_limits<double>::infinity();
BOOST_FOREACH(Polyhedron::Vertex_handle v, vertices(poly))
{
Point p = v->point();
max_coord = (std::max)(max_coord, p.x());
max_coord = (std::max)(max_coord, p.y());
max_coord = (std::max)(max_coord, p.z());
}
return max_coord;
}
int main(int argc, char* argv[])
{
const char* filename = (argc > 1) ? argv[1] : "data/eight.off";
std::ifstream input(filename);
Polyhedron poly;
if (!input || !(input >> poly) || poly.empty()
{
std::cerr << "Not a valid input file." << std::endl;
return 1;
}
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;
}

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.

Orientation

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

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 Polygon Soup Example puts these functions at work on a polygon soup.

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).

Combinatorial Repairing

Stitching

It happens that a polygon 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 such a polygon mesh consists in two main steps. First, border edges that are similar but duplicated are detected and paired. Then, they are "stitched" together so that the edges and vertices duplicates are removed from the mesh, and each of these remaining edges is incident to exactly two faces.

The function CGAL::Polygon_mesh_processing::stitch_borders() performs such repairing operation. 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 <iostream>
#include <fstream>
typedef CGAL::Polyhedron_3<K> Polyhedron;
int main(int argc, char* argv[])
{
const char* filename = (argc > 1) ? argv[1] : "data/full_border_quads.off";
std::ifstream input(filename);
Polyhedron mesh;
if (!input || !(input >> mesh) || mesh.is_empty()) {
std::cerr << "Not a valid off file." << 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;
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;
std::ofstream output("mesh_stitched.off");
output << std::setprecision(17) << mesh;
return 0;
}