CGAL 6.2 - Polygon Mesh Processing
Loading...
Searching...
No Matches
User Manual

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

Note
Starting with CGAL 6.2, the documentation for the "Polygon Mesh Processing" package has been reorganized into several specialized packages. This package retains the core functionalities, while advanced and specialized features have been moved to dedicated packages:

  • Boolean Operations On Meshes: algorithms for Boolean operations on polygon meshes, including clipping, splitting, and slicing with planes, boxes, or other meshes.
  • Meshing and Remeshing of Polygon Meshes: algorithms for meshing and remeshing, such as triangulation, refinement, simplification, optimization, and smoothing.
  • Polygon Mesh Repair: tools for detecting and correcting combinatorial and geometric defects in polygon meshes and polygon soups, including face orientation, hole filling, removal of degeneracies, and boundary stitching.

Introduction

This package provides a comprehensive set of methods and classes for polygon mesh processing, ranging from basic operations on mesh elements to advanced geometry processing algorithms. The implementation is primarily based on the algorithms and references presented 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. In this package, a polygon mesh is considered to have the topology of a 2-manifold. Note that all these requirements are mostly combinatorial, and do not impose any geometric constraints on the shape of the polygons. As such, this definition does not prevent the presence of defects such as self-intersections, degenerate faces or edges, etc.

API

This package follows the BGL API described in CGAL and the Boost Graph Library. It can thus be used either with CGAL::Surface_mesh, CGAL::Polyhedron_3, 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:

  • Predicates : predicates that can be evaluated on the polygon mesh, which includes point location and self-intersection tests.
  • Connected Components : methods to deal with connected components of a polygon mesh (extraction, marks, removal, ...).
  • Computing Normals : normal computation at the vertices and on the faces of a polygon mesh.
  • Feature Detection : methods to detect sharp geometric features on a polygon mesh.
  • Computing Curvatures : computation of curvatures (mean, gaussian, principal) on a polygon mesh.
  • Hausdorff Distance : methods to compute distances between polygon meshes.

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.

Predicates

This package provides several predicates to determine the characteristics of a triangle mesh or a subset of its faces.

Intersections Detection

Intersection tests between triangle meshes and/or polylines can be performed using the function CGAL::Polygon_mesh_processing::do_intersect() . Additionally, the function CGAL::Polygon_mesh_processing::intersecting_meshes() can be used to collect all pairs of intersecting meshes within 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 demonstrates self-intersection detection in the pig.off mesh. Detected self-intersections are illustrated in Figure 75.2.

Example: Polygon_mesh_processing/self_intersections_example.cpp

Show / Hide
#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 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;
}
bool is_triangle_mesh(const FaceGraph &g)
std::string data_file_path(const std::string &filename)

Figure 75.2 Detection of self-intersections in a triangle mesh. Intersecting triangles are displayed in dark grey and red in the right image.


Side of Triangle Mesh

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

A point is considered to be on the bounded side of the mesh if an odd number of surfaces are crossed when moving from the point to infinity.

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

Example: Polygon_mesh_processing/point_inside_example.cpp

Show / Hide
#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 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;
}
double size = max_coordinate(poly);
unsigned int nb_points = 100;
std::vector<Point> points;
points.reserve(nb_points);
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;
}
This class provides an efficient point location functionality with respect to a domain bounded by one...
Definition: Side_of_triangle_mesh.h:76
double to_double(const NT &x)
bool is_empty(const FaceGraph &g)
Bounded_side
ON_BOUNDED_SIDE

Polyhedral Envelope Containment Check

The class CGAL::Polyhedral_envelope provides functors to check whether a query point, segment, or triangle is fully contained within a polyhedral envelope of a triangle mesh or triangle soup.

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

Given a distance \( \delta = \epsilon / \sqrt(3)\), a prism is associated with each triangle by intersecting halfspaces parallel and orthogonal to the triangle and its edges, and additional halfspaces for obtuse angles, with face normals corresponding to angle bisectors. These halfspaces are at distance \( \delta \) and 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.3 The prism for a single triangle (left), the polyhedral envelope (middle), and the Minkowski sum envelope (right) for a triangle mesh.


The polyhedral envelope is guaranteed to be contained within the Minkowski sum envelope. The containment test is exact for the polyhedral envelope and conservative for the Minkowski sum envelope: if a query is inside the polyhedral envelope, it is also inside the Minkowski sum envelope; if outside, its relation to the Minkowski sum envelope is undetermined.

The algorithm of Wang et al. [10] for polyhedral envelope containment proceeds 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 algorithm checks if the query is completely covered. The details of how to check this covering can be found in the paper.

Polyhedral envelope containment is used by Surface_mesh_simplification::Polyhedral_envelope_filter in the Triangulated Surface Mesh Simplification package to simplify triangle meshes within a given tolerance.

Polyhedral Envelope Examples

The following example demonstrates construction of a polyhedral envelope for a CGAL::Surface_mesh and performing queries.

Example: Polygon_mesh_processing/polyhedral_envelope.cpp

Show / Hide
#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 Kernel::Point_3 Point_3;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
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;
}
This class can be used to check if a query point, segment, or triangle is inside or outside a polyhed...
Definition: Polyhedral_envelope.h:106

As connectivity information is not required, the same check can be performed on a triangle soup.

Example: Polygon_mesh_processing/polyhedral_envelope_of_triangle_soup.cpp

Show / Hide
#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 Kernel::Point_3 Point_3;
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;
}
bool read_OFF(std::istream &is, Graph &g, const NamedParameters &np=parameters::default_values())

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

Example: Polygon_mesh_processing/polyhedral_envelope_mesh_containment.cpp

Show / Hide
#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 Kernel::Point_3 Point_3;
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;
}
void isotropic_remeshing(const FaceRange &faces, SizingFunction &sizing, PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
FT edge_length(typename boost::graph_traits< PolygonMesh >::halfedge_descriptor h, const PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
computes the length of an edge of a given polygon mesh.
Definition: measure.h:109

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.

Connected Components

Collecting Connected Components

Functions are provided to enumerate and store the connected components of a polygon mesh. Connected components may be closed and geometrically separated, or separated by border or user-specified constraint edges.

The main entry point is the function CGAL::Polygon_mesh_processing::connected_components(), which collects all the connected components and fills a property map with the indices of the different connected components.

If a single connected component is to be extracted, the function CGAL::Polygon_mesh_processing::connected_component() collects all the faces that belong to the same connected component as the face that is provided as a parameter.

When a 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.

Modifying Connected Components

It is often useful to retain or remove specific connected components, for example, to discard small noisy components in favor of larger ones.

The functions CGAL::Polygon_mesh_processing::keep_connected_components() and CGAL::Polygon_mesh_processing::remove_connected_components() enable the user to keep or 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).

Finally, it can be useful to quickly remove some connected components based on characteristics of the surface mesh. 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.

Finally, CGAL::Polygon_mesh_processing::split_connected_components() splits the mesh into separate meshes for each connected component.

Connected Components Examples

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 means to demarcate the border of a connected component, and prevents the propagation of a connected component index to cross it.

Example: Polygon_mesh_processing/connected_components_example.cpp

Show / Hide
#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 Kernel::Point_3 Point;
typedef Kernel::Compare_dihedral_angle_3 Compare_dihedral_angle_3;
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;
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;
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;
4,
CGAL::parameters::edge_is_constrained_map(Constraint<Mesh>(mesh, bound)));
std::cerr << "- We keep the two largest components" << std::endl;
2,
CGAL::parameters::edge_is_constrained_map(Constraint<Mesh>(mesh, bound)));
return 0;
}
boost::property_traits< FaceComponentMap >::value_type connected_components(const PolygonMesh &pmesh, FaceComponentMap fcm, const NamedParameters &np=parameters::default_values())
computes for each face the index of the corresponding connected component.
Definition: connected_components.h:199
std::size_t keep_large_connected_components(PolygonMesh &pmesh, const ThresholdValueType threshold_value, const NamedParameters &np=parameters::default_values())
removes connected components whose size is (strictly) smaller than a given threshold value,...
Definition: connected_components.h:535
FaceOutputIterator connected_component(typename boost::graph_traits< PolygonMesh >::face_descriptor seed_face, const PolygonMesh &pmesh, FaceOutputIterator out, const NamedParameters &np=parameters::default_values())
discovers all the faces in the same connected component as seed_face and records them in out.
Definition: connected_components.h:119
std::size_t keep_largest_connected_components(PolygonMesh &pmesh, std::size_t nb_components_to_keep, const NamedParameters &np=parameters::default_values())
removes the small connected components and all isolated vertices.
Definition: connected_components.h:384
Oriented_side opposite(const Oriented_side &o)

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

Example: Polygon_mesh_processing/face_filtered_graph_example.cpp

Show / Hide
#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 Kernel::Point_3 Point;
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;
}

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:

Surface Location Example

The following example demonstrates usage of these functions.

Example: Polygon_mesh_processing/locate_example.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/locate.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits_3.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/Dynamic_property_map.h>
typedef K::FT FT;
typedef K::Point_2 Point_2;
typedef K::Ray_2 Ray_2;
typedef K::Point_3 Point_3;
typedef K::Ray_3 Ray_3;
typedef typename boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<Mesh>::face_descriptor face_descriptor;
namespace CP = CGAL::parameters;
namespace PMP = CGAL::Polygon_mesh_processing;
typedef PMP::Barycentric_coordinates<FT> Barycentric_coordinates;
typedef PMP::Face_location<Mesh, FT> Face_location;
int main(int /*argc*/, char** /*argv*/)
{
// Generate a simple 3D triangle mesh that with vertices on the plane xOy
Mesh tm;
CGAL::make_grid(10, 10, tm);
// Basic usage
Face_location random_location = PMP::random_location_on_mesh<FT>(tm);
const face_descriptor f = random_location.first;
const Barycentric_coordinates& coordinates = random_location.second;
std::cout << "Random location on the mesh: face " << f
<< " and with coordinates [" << coordinates[0] << "; "
<< coordinates[1] << "; "
<< coordinates[2] << "]\n";
std::cout << "It corresponds to point (" << PMP::construct_point(random_location, tm) << ")\n\n";
// Locate a known 3D point in the mesh
const Point_3 query(1.2, 7.4, 0);
Face_location query_location = PMP::locate(query, tm);
std::cout << "Point (" << query << ") is located in face " << query_location.first
<< " with barycentric coordinates [" << query_location.second[0] << "; "
<< query_location.second[1] << "; "
<< query_location.second[2] << "]\n\n";
// Locate a 3D point in the mesh as the intersection of the mesh and a 3D ray.
// The AABB tree can be cached in case many queries are performed (otherwise, it is rebuilt
// on each call, which is expensive).
typedef CGAL::AABB_face_graph_triangle_primitive<Mesh> AABB_face_graph_primitive;
typedef CGAL::AABB_traits_3<K, AABB_face_graph_primitive> AABB_face_graph_traits;
const Ray_3 ray_3(Point_3(4.2, 6.8, 2.4), Point_3(7.2, 2.3, -5.8));
Face_location ray_location = PMP::locate_with_AABB_tree(ray_3, tree, tm);
std::cout << "Intersection of the 3D ray and the mesh is in face " << ray_location.first
<< " with barycentric coordinates [" << ray_location.second[0] << " "
<< ray_location.second[1] << " "
<< ray_location.second[2] << "]\n";
std::cout << "It corresponds to point (" << PMP::construct_point(ray_location, tm) << ")\n";
std::cout << "Is it on the face's border? " << (PMP::is_on_face_border(ray_location, tm) ? "Yes" : "No") << "\n\n";
// -----------------------------------------------------------------------------------------------
// Now, we artificially project the mesh to the natural 2D dimensional plane, with a little translation
// via a custom vertex point property map
typedef CGAL::dynamic_vertex_property_t<Point_2> Point_2_property;
typedef typename boost::property_map<Mesh, Point_2_property>::type Projection_pmap;
Projection_pmap projection_pmap = get(Point_2_property(), tm);
for(vertex_descriptor v : vertices(tm))
{
const Point_3& p = tm.point(v);
put(projection_pmap, v, Point_2(p.x() + 1, p.y())); // simply ignoring the z==0 coordinate and translating along Ox
}
// Locate the same 3D point but in a 2D context
const Point_2 query_2(query.x() + 1, query.y());
Face_location query_location_2 = PMP::locate(query_2, tm, CP::vertex_point_map(projection_pmap));
std::cout << "Point (" << query_2 << ") is located in face " << query_location_2.first
<< " with barycentric coordinates [" << query_location_2.second[0] << "; "
<< query_location_2.second[1] << "; "
<< query_location_2.second[2] << "]\n\n";
// Shoot a 2D ray and locate the intersection with the mesh in 2D
const Ray_2 ray_2(Point_2(-10, -10), Point_2(10, 10));
Face_location ray_location_2 = PMP::locate(ray_2, tm, CP::vertex_point_map(projection_pmap)); // This rebuilds an AABB tree on each call
std::cout << "Intersection of the 2D ray and the mesh is in face " << ray_location_2.first
<< " with barycentric coordinates [" << ray_location_2.second[0] << "; "
<< ray_location_2.second[1] << "; "
<< ray_location_2.second[2] << "]\n";
std::cout << "It corresponds to point (" << PMP::construct_point(ray_location_2, tm, CP::vertex_point_map(projection_pmap)) << ")\n";
if(PMP::is_on_mesh_border(ray_location_2, tm))
std::cout << "It is on the border of the mesh!\n" << std::endl;
return EXIT_SUCCESS;
}
bool triangulate_faces(FaceRange face_range, PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
std::pair< typename boost::graph_traits< TriangleMesh >::face_descriptor, Barycentric_coordinates< FT > > Face_location
If tm is the input triangulated surface mesh and given the pair (f, bc) such that bc is the triplet o...
Definition: locate.h:83
bool is_on_face_border(const Face_location< TriangleMesh, FT > &loc, const TriangleMesh &tm)
Given a location, returns whether the location is on the boundary of the face or not.
Definition: locate.h:790
Face_location< TriangleMesh, FT > locate_with_AABB_tree(const Point &p, const AABB_tree< AABB_traits_3< Geom_traits, AABB_face_graph_triangle_primitive< TriangleMesh, Point3VPM > > > &tree, const TriangleMesh &tm, const NamedParameters &np=parameters::default_values())
returns the face location nearest to the given point, as a location.
Definition: locate.h:1599
std::array< FT, 3 > Barycentric_coordinates
A triplet of coordinates describing the barycentric coordinates of a point with respect to the vertic...
Definition: locate.h:71
Point construct_point(const Face_location< TriangleMesh, FT > &loc, const TriangleMesh &tm, const NamedParameters &np=parameters::default_values())
Given a location in a face, returns the geometric position described by these coordinates,...
Definition: locate.h:591
Face_location< TriangleMesh, FT > locate(const Point &p, const TriangleMesh &tm, const NamedParameters &np=parameters::default_values())
returns the nearest face location to the given point.
Definition: locate.h:1688
bool is_on_mesh_border(const Face_location< TriangleMesh, FT > &loc, const TriangleMesh &tm)
Given a location, returns whether the location is on the border of the mesh or not.
Definition: locate.h:826
void build_AABB_tree(const TriangleMesh &tm, AABB_tree< AABB_traits_3< Geom_traits, CGAL::AABB_face_graph_triangle_primitive< TriangleMesh, Point3VPM > > > &outTree, const NamedParameters &np=parameters::default_values())
creates an AABB tree suitable for use with locate_with_AABB_tree().
Definition: locate.h:1524
boost::graph_traits< Graph >::halfedge_descriptor make_grid(typename boost::graph_traits< Graph >::vertices_size_type i, typename boost::graph_traits< Graph >::vertices_size_type j, Graph &g, const CoordinateFunctor &calculator, bool triangulated=false)

Computing Normals

Methods are provided to compute normals on polygon meshes, either per face or per vertex:

When computing all the normals of faces and vertices, the following functions should be preferred as they factorize some computations:

Computing Normals Example

In the following examples we associate a normal vector to each vertex and to each face of a mesh of type CGAL::Surface_mesh.

Example: Polygon_mesh_processing/compute_normals_example.cpp

Show / Hide
#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 K::Point_3 Point;
typedef K::Vector_3 Vector;
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;
}
void compute_normals(const PolygonMesh &pmesh, VertexNormalMap vertex_normals, FaceNormalMap face_normals, const NamedParameters &np=parameters::default_values())
computes the outward unit vector normal for all vertices and faces of the polygon mesh.
Definition: compute_normal.h:901
const CGAL::Null_vector NULL_VECTOR

Computing Curvatures

This package provides methods to compute curvatures on polygon meshes based on "Interpolated Corrected Curvatures on Polyhedral Surfaces" [6]. 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. Furthermore, they give accurate results on meshes with noise on vertex positions, 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 [7] and [6]. They introduce a new way to compute curvatures on polygon meshes. The main idea in [7] 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. [6] 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 scalar product, and \( \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 with respect to mesh data structure and can be used with CGAL::Surface_mesh, CGAL::Polyhedron_3, or any polygon mesh structure meeting the FaceGraph concept requirements.

Curvatures are computed for all vertices using CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures(), with named parameters to select which curvatures (and directions) to compute. An overload is available for computing curvatures at a single vertex.

Results

Figure 75.4 illustrates various curvature measures on a triangular mesh.

(a)(b)

Figure 75.4 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.5 Varying the integration ball radius yields a scale space of curvature measures, useful for handling noise in the input mesh. The second row illustrates mean curvature 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.6 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 computed curvatures, as shown in the examples. For each property map, a curvature value is associated with each vertex.

Interpolated Corrected Curvatures on a Surface Mesh Example

The following example demonstrates computation of curvatures at vertices and storage in property maps provided by CGAL::Surface_mesh.

Example: Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp

Show / Hide
#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 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);
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;
}
void interpolated_corrected_curvatures(const PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
computes the interpolated corrected curvatures across the mesh pmesh.
Definition: interpolated_corrected_curvatures.h:1083
bool read_polygon_mesh(const std::string &fname, Graph &g, const NamedParameters &np=parameters::default_values())
a struct for storing principal curvatures and directions.
Definition: interpolated_corrected_curvatures.h:43

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 CGAL::Polyhedron_3 does not provide storage for the curvatures.

Example: Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp

Show / Hide
#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 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 =
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++;
}
}

Interpolated Corrected Curvatures on a Vertex Example

The following example demonstrates computation of curvatures at a specific vertex.

Example: Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp

Show / Hide
#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 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;
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;
}
GT::FT min_curvature
min curvature magnitude
Definition: interpolated_corrected_curvatures.h:46
GT::FT max_curvature
max curvature magnitude
Definition: interpolated_corrected_curvatures.h:49

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. [8]. These curvatures are computed at each vertex of the mesh, and are based on the angles of the incident triangles. The functions are:

Feature Detection

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

The function CGAL::Polygon_mesh_processing::sharp_edges_segmentation() detects sharp edges and deduces surface patches and vertex incidences. It is composed of three functions:

These functions respectively detect sharp edges, compute patch indices, and assign patch indices to each vertex based on incident faces.

The following example counts the number of edges incident to two faces whose normals form an angle less than 90 degrees, and the number of surface patches separated by these edges.

Example: Polygon_mesh_processing/detect_features_example.cpp

Show / Hide
#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 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;
}
boost::graph_traits< PolygonMesh >::faces_size_type sharp_edges_segmentation(const PolygonMesh &pmesh, FT angle_in_deg, EdgeIsFeatureMap edge_is_feature_map, PatchIdMap patch_id_map, const NamedParameters &np=parameters::default_values())
This function calls successively CGAL::Polygon_mesh_processing::detect_sharp_edges(),...
Definition: detect_features.h:460

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 [5]. 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.7).

Figure 75.7 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, lower and upper bounds of the Hausdorff distance to the point set are computed. Triangles are refined until the difference between bounds is below 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.

Example: Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp

Show / Hide
#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 K::Point_3 Point;
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;
std::cout << "Approximated Hausdorff distance: "
<TAG>(tm1, tm2, CGAL::parameters::number_of_points_per_area_unit(4000))
<< std::endl;
return 0;
}
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:1176
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 as follows:

// computes the approximation error of the reconstruction
double max_dist =
(output_mesh,
CGAL::make_range (boost::make_transform_iterator
boost::make_transform_iterator
4000);
std::cout << "Max distance to point_set: " << max_dist << std::endl;
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:1255
Iterator_range< T > make_range(const T &b, const T &e)

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 [9]. 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 a 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.

Example: Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp

Show / Hide
#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 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 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;
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;
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;
}
void transform(const Transformation &transformation, PolygonMesh &mesh, const NamedParameters &np=parameters::default_values())
applies a transformation to every vertex of a PolygonMesh.
Definition: transform.h:50

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.

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.

Interpolated corrected curvatures were implemented during GSoC 2022 by Hossam Saeed, supervised by David Coeurjolly, Jacques-Olivier Lachaud, and Sébastien Loriot. The implementation is based on [6]. DGtal's implementation was also referenced during development.