Authors: Sébastien Loriot, Mael Rouxel-Labbé, and Ilker O. Yaz

Introduction

Geometric models, whether acquired through imprecise tools or imperfect processes, frequently exhibit defects such as inconsistent orientations, degeneracies, gaps, missing data, non-manifold features, or self-intersections.

This CGAL package provides a comprehensive set of functions to detect and address both combinatorial and geometric defects. The first part of this documentation focuses on combinatorial repairs, with particular attention to orienting and repairing polygon soups—a necessary step to obtain a valid combinatorial polygon mesh structure. The second part describes functions for geometric repairs, such as removing nearly degenerate faces, or filling holes.

Orientation

This package offers multiple functions to compute consistent face orientations for set of faces (Section Polygon Soups) and polygon meshes (see Section Polygon Meshes).

Polygon Soups

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

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

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

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

Example: Polygon_mesh_processing/orient_polygon_soup_example.cpp

Show / Hide

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedron_items_with_id_3.h>
 
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
 
#include <CGAL/IO/polygon_soup_io.h>
 
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel          K;
typedef CGAL::Polyhedron_3<K, CGAL::Polyhedron_items_with_id_3>      Mesh;
 
// Optional visitor for orientating a polygon soup to demonstrate usage for some functions.
// inherits from the default class as some functions are not overloaded
struct Visitor : public CGAL::Polygon_mesh_processing::Default_orientation_visitor
{
  void non_manifold_edge(std::size_t id1, std::size_t id2, std::size_t nb_poly)
  {
    std::cout << "The edge " << id1 << ", " << id2 << " is not manifold: " << nb_poly << " incident polygons." << std::endl;
  }
  void non_manifold_vertex(std::size_t id, std::size_t nb_cycles)
  {
     std::cout << "The vertex " << id << " is not manifold: "  << nb_cycles << " connected components of vertices in the link." << std::endl;
  }
  void duplicated_vertex(std::size_t v1, std::size_t v2)
  {
    std::cout << "The vertex " << v1 << " has been duplicated, its new id is " << v2 << "." << std::endl;
  }
  void vertex_id_in_polygon_replaced(std::size_t p_id, std::size_t i1, std::size_t i2)
  {
    std::cout << "In the polygon " << p_id << ", the index " << i1 << " has been replaced by " << i2 << "." << std::endl;
  }
  void polygon_orientation_reversed(std::size_t p_id)
  {
    std::cout << "The polygon " << p_id << " has been reversed." << std::endl;
  }
};
 
int main(int argc, char* argv[])
{
  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/tet-shuffled.off");
 
  std::vector<K::Point_3> points;
  std::vector<std::vector<std::size_t> > polygons;
 
  if(!CGAL::IO::read_polygon_soup(filename, points, polygons) || points.empty())
  {
    std::cerr << "Cannot open file " << std::endl;
    return EXIT_FAILURE;
  }
 
  Visitor visitor;
  CGAL::Polygon_mesh_processing::orient_polygon_soup(points, polygons, CGAL::parameters::visitor(visitor));
 
  Mesh mesh;
  CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh);
 
  // Number the faces because 'orient_to_bound_a_volume' needs a face <--> index map
  int index = 0;
  for(Mesh::Face_iterator fb=mesh.facets_begin(), fe=mesh.facets_end(); fb!=fe; ++fb)
    fb->id() = index++;
 
  if(CGAL::is_closed(mesh))
    CGAL::Polygon_mesh_processing::orient_to_bound_a_volume(mesh);
 
  std::ofstream out("tet-oriented1.off");
  out.precision(17);
  out << mesh;
  out.close();
 
  CGAL::Polygon_mesh_processing::reverse_face_orientations(mesh);
  std::ofstream out2("tet-oriented2.off");
  out2.precision(17);
  out2 << mesh;
  out2.close();
 
  return EXIT_SUCCESS;
}

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

Polygon Meshes

This package provides functions for orienting faces in a closed polygon mesh:

CGAL::Polygon_mesh_processing::orient() orients each connected component of a closed polygon mesh outward or inward.
CGAL::Polygon_mesh_processing::orient_to_bound_a_volume() orients the connected components of a closed polygon mesh so that they bound a volume (see Definitions for a precise definition).
CGAL::Polygon_mesh_processing::is_outward_oriented() checks whether an oriented polygon mesh is oriented such that the normals to all faces point outward from the domain bounded by the mesh.
CGAL::Polygon_mesh_processing::reverse_face_orientations() reverses the orientation of halfedges around faces, thereby reversing the computed normals (see Section Computing Normals).
CGAL::Polygon_mesh_processing::volume_connected_components() provides information about the 3D arrangement of surface connected components in a triangle mesh, with many named parameter options, making it a generalization of is_outward_oriented().
CGAL::Polygon_mesh_processing::duplicate_non_manifold_edges_in_polygon_soup() duplicates points and edges to make a soup orientable, without altering face orientations.
CGAL::Polygon_mesh_processing::orient_triangle_soup_with_reference_triangle_mesh() orients the triangles of a soup according to a reference mesh.
CGAL::Polygon_mesh_processing::merge_reversible_connected_components() merges the connected components of a polygon mesh when possible.

The following example demonstrates how to repair and orient a soup to obtain a mesh from a reference:

Example: Polygon_mesh_processing/orientation_pipeline_example.cpp

Show / Hide

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

Combinatorial Repair

Polygon Soup Repairing

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

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

Stitching

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

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

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

The input mesh should be manifold; otherwise, stitching may not succeed.

Stitching Example

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

Example: Polygon_mesh_processing/stitch_borders_example.cpp

Show / Hide

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

Polygon Mesh Manifoldness

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

Manifoldness Repair Example

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

Example: Polygon_mesh_processing/manifoldness_repair_example.cpp

Show / Hide

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

Duplicated Vertices in Boundary Cycles

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

Geometric Repair

Removal of Almost Degenerate Triangle Faces

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

Self-intersection Resolution (Autorefinement) in Triangle Soups

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

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

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 [2] 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 [3], 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 complex hole boundaries, the generated patch may have self-intersections. After hole filling, the patch can be refined and faired using the meshing functions CGAL::Polygon_mesh_processing::refine() and CGAL::Polygon_mesh_processing::fair() (see Section Chapter_PMPRemeshing).

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

API

This package provides four functions for hole filling:

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

Examples

Triangulate a Polyline

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

Example: Polygon_mesh_processing/triangulate_polyline_example.cpp

Show / Hide

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

Hole Filling From the Border of the Hole

If the input polygon mesh contains one or more holes, they can be filled iteratively by detecting border edges (edges with only one incident non-null face) after each filling step.

Holes are filled sequentially, and the process stops when no border edge remains.

The example below illustrates this process, where holes are iteratively filled, refined, and faired. Optionally, only holes not exceeding a specified diameter or number of edges can be filled. This example assumes the mesh is stored in a CGAL::Surface_mesh data structure. Analogous examples for CGAL::Polyhedron_3 and other classes are available in the code base.

Example: Polygon_mesh_processing/hole_filling_example_SM.cpp

Show / Hide

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

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

An additional parameter, visitor, can be used to track the algorithm's phases, enabling users to implement timeouts or monitor progress.

Example: Polygon_mesh_processing/hole_filling_visitor_example.cpp

Show / Hide

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
 
#include <CGAL/boost/graph/border.h>
#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Real_timer.h>
 
#include <boost/lexical_cast.hpp>
 
#include <iostream>
#include <iterator>
#include <string>
#include <tuple>
#include <vector>
#include <stdexcept>
 
typedef CGAL::Real_timer                                    Timer;
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3                                     Point;
typedef CGAL::Surface_mesh<Point>                           Mesh;
 
typedef boost::graph_traits<Mesh>::vertex_descriptor        vertex_descriptor;
typedef boost::graph_traits<Mesh>::halfedge_descriptor      halfedge_descriptor;
typedef boost::graph_traits<Mesh>::face_descriptor          face_descriptor;
 
namespace PMP = CGAL::Polygon_mesh_processing;
 
bool is_small_hole(halfedge_descriptor h, Mesh & mesh,
                   double max_hole_diam, int max_num_hole_edges)
{
  int num_hole_edges = 0;
  CGAL::Bbox_3 hole_bbox;
  for (halfedge_descriptor hc : CGAL::halfedges_around_face(h, mesh))
  {
    const Point& p = mesh.point(target(hc, mesh));
 
    hole_bbox += p.bbox();
    ++num_hole_edges;
 
    // Exit early, to avoid unnecessary traversal of large holes
    if (num_hole_edges > max_num_hole_edges) return false;
    if (hole_bbox.xmax() - hole_bbox.xmin() > max_hole_diam) return false;
    if (hole_bbox.ymax() - hole_bbox.ymin() > max_hole_diam) return false;
    if (hole_bbox.zmax() - hole_bbox.zmin() > max_hole_diam) return false;
  }
 
  return true;
}
 
struct Stop : std::exception
{
    Stop()
    {}
};
 
struct Progress :
    public PMP::Hole_filling::Default_visitor
{
  Progress(double time_limit)
      : time_limit(time_limit)
  {}
 
  Progress(const Progress&) = delete;
 
  void start_planar_phase() const
  {
    std::cout << "Start planar phase"<< std::endl;
  }
 
  void end_planar_phase(bool success) const
  {
    std::cout << "End planar phase " << (success? "(success)" : "(failed)") << std::endl;
  }
 
  void start_quadratic_phase(std::size_t n)
  {
    timer.start();
    quadratic_i = 0;
    quadratic_n = n;
    quadratic_report = n / 10;
    std::cout << "Start quadratic phase with estimated " << n << " steps" << std::endl;
  }
 
  void quadratic_step()
  {
    if (quadratic_i++ == quadratic_report) {
      std::cout << double(quadratic_i) / double(quadratic_n) * 100 << "%" << std::endl;
      quadratic_report += quadratic_n / 10;
    }
  }
 
  void end_quadratic_phase(bool success) const
  {
    timer.stop();
    std::cout << "End quadratic phase " << timer.time() << " sec. " << (success ? "(success)" : "(failed)") << std::endl;
    timer.reset();
  }
 
 
  void start_cubic_phase(std::size_t n)
  {
    timer.start();
      cubic_n = n;
      cubic_report = n / 10;
      std::cout << "Start cubic phase with " << n << " steps" << std::endl;
  }
 
 
  void cubic_step()
  {
    if (timer.time() > time_limit) {
        std::cout << "Let's stop here" << std::endl;
        throw Stop();
    }
    if (cubic_i++ == cubic_report) {
      std::cout << double(cubic_i) / double(cubic_n) * 100 << "%" << std::endl;
      cubic_report += cubic_n / 10;
    }
  }
 
  void end_cubic_phase() const
  {
    std::cout << "End cubic phase " << timer.time() << " sec. "  << std::endl;
  }
 
  mutable Timer timer;
  double time_limit;
  std::size_t quadratic_n = 0, quadratic_i = 0, quadratic_report = 0;
  std::size_t cubic_n = 0, cubic_i = 0, cubic_report = 0;
};
 
// Incrementally fill the holes that are no larger than given diameter
// and with no more than a given number of edges (if specified).
 
int main(int argc, char* argv[])
{
  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/mech-holes-shark.off");
 
  Mesh mesh;
  if(!PMP::IO::read_polygon_mesh(filename, mesh))
  {
    std::cerr << "Invalid input." << std::endl;
    return 1;
  }
 
  // Both of these must be positive in order to be considered
  double max_hole_diam   = (argc > 2) ? boost::lexical_cast<double>(argv[2]): -1.0;
  int max_num_hole_edges = (argc > 3) ? boost::lexical_cast<int>(argv[3]) : -1;
 
  unsigned int nb_holes = 0;
  std::vector<halfedge_descriptor> border_cycles;
 
  // collect one halfedge per boundary cycle
  CGAL::extract_boundary_cycles(mesh, std::back_inserter(border_cycles));
 
  for(halfedge_descriptor h : border_cycles)
  {
    if(max_hole_diam > 0 && max_num_hole_edges > 0 &&
       !is_small_hole(h, mesh, max_hole_diam, max_num_hole_edges))
      continue;
 
    Progress progress(10.0);
    bool success = false;
 
    try {
        success = std::get<0>(PMP::triangulate_refine_and_fair_hole(mesh,
            h,
            CGAL::parameters::visitor(std::ref(progress)).use_delaunay_triangulation(true)));
    }
    catch (const Stop&) {
        std::cout << "We stopped with a timeout" << std::endl;
    }
    std::cout << "  Is fairing successful: " << success << std::endl;
    ++nb_holes;
  }
 
  std::cout << std::endl;
  std::cout << nb_holes << " holes have been filled" << std::endl;
 
  CGAL::IO::write_polygon_mesh("filled_SM.off", mesh, CGAL::parameters::stream_precision(17));
  std::cout << "Mesh written to: filled_SM.off" << std::endl;
 
  return 0;
}

Performance

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

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

The following running times were observed:

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

Implementation History

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