Merge pull request #4749 from danston/Polygon_mesh_processing-triangulate_hole_with_cdt2-danston

BGL/include/CGAL/boost/graph/parameters_interface.h

@@ -61,6 +61,7 @@ CGAL_add_named_parameter(geom_traits_t, geom_traits, geom_traits)
 CGAL_add_named_parameter(vertex_incident_patches_t, vertex_incident_patches, vertex_incident_patches_map)
 CGAL_add_named_parameter(density_control_factor_t, density_control_factor, density_control_factor)
 CGAL_add_named_parameter(use_delaunay_triangulation_t, use_delaunay_triangulation, use_delaunay_triangulation)
+CGAL_add_named_parameter(use_2d_constrained_delaunay_triangulation_t, use_2d_constrained_delaunay_triangulation, use_2d_constrained_delaunay_triangulation)
 CGAL_add_named_parameter(fairing_continuity_t, fairing_continuity, fairing_continuity)
 CGAL_add_named_parameter(sparse_linear_solver_t, sparse_linear_solver, sparse_linear_solver)
 CGAL_add_named_parameter(number_of_relaxation_steps_t, number_of_relaxation_steps, number_of_relaxation_steps)

BGL/test/BGL/test_cgal_bgl_named_params.cpp

@@ -63,6 +63,7 @@ void test(const NamedParameters& np)
   assert(get_parameter(np, CGAL::internal_np::vertex_incident_patches).v == 11);
   assert(get_parameter(np, CGAL::internal_np::density_control_factor).v == 12);
   assert(get_parameter(np, CGAL::internal_np::use_delaunay_triangulation).v == 13);
+  assert(get_parameter(np, CGAL::internal_np::use_2d_constrained_delaunay_triangulation).v == 4573);
   assert(get_parameter(np, CGAL::internal_np::fairing_continuity).v == 14);
   assert(get_parameter(np, CGAL::internal_np::sparse_linear_solver).v == 15);
   assert(get_parameter(np, CGAL::internal_np::number_of_relaxation_steps).v == 16);
@@ -165,6 +166,7 @@ void test(const NamedParameters& np)
   check_same_type<11>(get_parameter(np, CGAL::internal_np::vertex_incident_patches));
   check_same_type<12>(get_parameter(np, CGAL::internal_np::density_control_factor));
   check_same_type<13>(get_parameter(np, CGAL::internal_np::use_delaunay_triangulation));
+  check_same_type<4573>(get_parameter(np, CGAL::internal_np::use_2d_constrained_delaunay_triangulation));
   check_same_type<14>(get_parameter(np, CGAL::internal_np::fairing_continuity));
   check_same_type<15>(get_parameter(np, CGAL::internal_np::sparse_linear_solver));
   check_same_type<16>(get_parameter(np, CGAL::internal_np::number_of_relaxation_steps));
@@ -314,6 +316,7 @@ int main()
                          .vertex_incident_patches_map(A<11>(11))
                          .density_control_factor(A<12>(12))
                          .use_delaunay_triangulation(A<13>(13))
+                         .use_2d_constrained_delaunay_triangulation(A<4573>(4573))
                          .fairing_continuity(A<14>(14))
                          .sparse_linear_solver(A<15>(15))
                          .number_of_relaxation_steps(A<16>(16))

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polygon_mesh.h

@@ -94,7 +94,9 @@ triangulate_hole_polygon_mesh(PolygonMesh& pmesh,
             OutputIterator out,
             VertexPointMap vpmap,
             bool use_delaunay_triangulation,
-            const Kernel& k)
+            const Kernel& k,
+            const bool use_cdt,
+            const typename Kernel::FT max_squared_distance)
 {
   typedef Halfedge_around_face_circulator<PolygonMesh>   Hedge_around_face_circulator;
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
@@ -173,6 +175,13 @@ triangulate_hole_polygon_mesh(PolygonMesh& pmesh,
   // fill hole using polyline function, with custom tracer for PolygonMesh
   Tracer_polyhedron<PolygonMesh, OutputIterator>
     tracer(out, pmesh, P_edges);
+
+#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+    if(use_cdt && triangulate_hole_polyline_with_cdt(P, tracer, is_valid, k, max_squared_distance))
+    {
+      return std::make_pair(tracer.out, CGAL::internal::Weight_min_max_dihedral_and_area(0,0));
+    }
+#endif
   CGAL::internal::Weight_min_max_dihedral_and_area weight =
     triangulate_hole_polyline(P, Q, tracer, WC(is_valid),
       use_delaunay_triangulation, k)

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polyline.h

@@ -22,6 +22,16 @@
 #include <CGAL/Delaunay_triangulation_cell_base_3.h>
 #include <CGAL/Triangulation_vertex_base_with_info_3.h>
 #endif
+
+#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+#include <CGAL/centroid.h>
+#include <CGAL/Triangulation_vertex_base_with_info_2.h>
+#include <CGAL/Triangulation_face_base_with_info_2.h>
+#include <CGAL/Constrained_Delaunay_triangulation_2.h>
+#include <CGAL/Triangulation_2_projection_traits_3.h>
+#include <queue>
+#endif
+
 #include <CGAL/utility.h>
 #include <CGAL/iterator.h>
 #include <CGAL/use.h>
@@ -1208,7 +1218,258 @@ public:
   }
 };
 
-/***********************************************************************************
+#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+
+/************************************************************************
+ * Triangulate hole by using a cdt_2
+ ************************************************************************/
+// /!\ points.first == points.last
+
+template <typename Traits>
+bool is_planar_2(
+  const std::vector<typename Traits::Point_3>& points,
+  const typename Traits::Vector_3& avg_normal,
+  const typename Traits::FT max_squared_distance,
+  const Traits& traits) {
+
+  typedef typename Traits::FT FT;
+  typedef typename Traits::Point_3 Point_3;
+  typedef typename Traits::Plane_3 Plane_3;
+  typedef typename Traits::Construct_projected_point_3 Projection_3;
+  typedef typename Traits::Compute_squared_distance_3 Squared_distance_3;
+
+  const Projection_3 projection_3 =
+    traits.construct_projected_point_3_object();
+  const Squared_distance_3 squared_distance_3 =
+    traits.compute_squared_distance_3_object();
+
+  const std::size_t n = points.size() - 1; // the first equals to the last
+  if (n < 3) {
+    return false; // cant be a plane!
+  }
+
+  // Compute centroid.
+  const Point_3 centroid =
+    CGAL::centroid(points.begin(), points.end() - 1);
+
+  // Compute covariance matrix.
+  FT xx = FT(0), yy = FT(0), zz = FT(0);
+  FT xy = FT(0), xz = FT(0), yz = FT(0);
+  for (std::size_t i = 0; i < n; ++i) {
+    const Point_3& p = points[i];
+    const FT dx = p.x() - centroid.x();
+    const FT dy = p.y() - centroid.y();
+    const FT dz = p.z() - centroid.z();
+    xx += dx * dx; yy += dy * dy; zz += dz * dz;
+    xy += dx * dy; xz += dx * dz; yz += dy * dz;
+  }
+
+  // Check the planarity.
+  const FT x = yy * zz - yz * yz;
+  const FT y = xx * zz - xz * xz;
+  const FT z = xx * yy - xy * xy;
+  FT maxv = -FT(1);
+  maxv = (CGAL::max)(maxv, x);
+  maxv = (CGAL::max)(maxv, y);
+  maxv = (CGAL::max)(maxv, z);
+  if (maxv <= FT(0)) {
+    return false; // a good plane does not exist for sure!
+  }
+
+  // Here, avg_squared_distance is a little bit more tolerant than avg_distance^2.
+  CGAL_assertion(avg_normal != typename Traits::Vector_3());
+  const Plane_3 plane = Plane_3(centroid, avg_normal);
+  FT avg_squared_distance = FT(0);
+  for (std::size_t i = 0; i < n; ++i) {
+    const Point_3& p = points[i];
+    const Point_3  q = projection_3(plane, p);
+    avg_squared_distance += CGAL::abs(squared_distance_3(p, q));
+  }
+  avg_squared_distance /= static_cast<FT>(n);
+  // std::cout << "avg squared distance: " << avg_squared_distance << std::endl;
+
+  CGAL_assertion(max_squared_distance >= FT(0));
+  if (avg_squared_distance > max_squared_distance) {
+    return false; // the user distance criteria are not satisfied!
+  }
+
+  // std::cout << "The hole seems to be near planar." << std::endl;
+  return true;
+}
+
+template <
+  typename PointRange, // need size()
+  typename Tracer,
+  typename Validity_checker,
+  typename Traits
+>
+bool
+triangulate_hole_polyline_with_cdt(const PointRange& points,
+                                   Tracer& tracer,
+                                   const Validity_checker& is_valid,
+                                   const Traits& traits,
+                                   const typename Traits::FT max_squared_distance)
+{
+  typedef typename Traits::FT FT;
+  typedef typename Traits::Point_3 Point_3;
+  typedef typename Traits::Vector_3 Vector_3;
+  typedef typename Traits::Collinear_3 Collinear_3;
+
+  // Compute an average normal of the hole.
+  const Collinear_3 collinear_3 =
+    traits.collinear_3_object();
+
+  std::vector<Point_3> P(std::begin(points), std::end(points));
+  CGAL_assertion(P.size() >= 3);
+  if (P.front() != P.back()) {
+    P.push_back(P.front());
+  }
+
+  FT x = FT(0), y = FT(0), z = FT(0);
+  std::size_t num_normals = 0;
+  const Point_3& ref_point = P[0];
+  const std::size_t size = P.size() - 1;
+  for (std::size_t i = 1; i < size - 1; ++i) {
+    const std::size_t ip = i + 1;
+
+    const Point_3& p1 = ref_point; // 3 points, which form a triangle
+    const Point_3& p2 = P[i];
+    const Point_3& p3 = P[ip];
+
+    // Skip in case we have collinear points.
+    if (collinear_3(p1, p2, p3)) {
+      continue;
+    }
+
+    // Computing the normal of a triangle.
+    const Vector_3 n = CGAL::normal(p1, p2, p3);
+    // If it is a positive normal ->
+    if (
+      ( n.x() >  FT(0)                                    ) ||
+      ( n.x() == FT(0) && n.y() >  FT(0)                  ) ||
+      ( n.x() == FT(0) && n.y() == FT(0) && n.z() > FT(0) )) {
+      x += n.x(); y += n.y(); z += n.z();
+    } else { // otherwise invert ->
+      x -= n.x(); y -= n.y(); z -= n.z();
+    }
+    ++num_normals;
+  }
+
+  if (num_normals < 1) {
+    // std::cerr << "WARNING: num normals, cdt 2 falls back to the original solution!" << std::endl;
+    return false;
+  }
+
+  // Setting the final normal.
+  x /= static_cast<FT>(num_normals);
+  y /= static_cast<FT>(num_normals);
+  z /= static_cast<FT>(num_normals);
+  const Vector_3 avg_normal = Vector_3(x, y, z);
+  // std::cout << "avg normal: " << avg_normal << std::endl;
+
+  // Checking the hole planarity.
+  if (!is_planar_2(P, avg_normal, max_squared_distance, traits)) {
+    // std::cerr << "WARNING: planarity, cdt 2 falls back to the original solution!" << std::endl;
+    return false;
+  }
+
+  // Checking the hole simplicity.
+  typedef Triangulation_2_projection_traits_3<Traits> P_traits;
+  const P_traits p_traits(avg_normal);
+  if (!is_simple_2(P.begin(), P.end() - 1, p_traits)) {
+    // std::cerr << "WARNING: simplicity, cdt 2 falls back to the original solution!" << std::endl;
+    return false;
+  }
+
+  Lookup_table_map<int> lambda(static_cast<int>(size), -1);
+
+  // Create and fill the cdt_2.
+  typedef CGAL::Triangulation_vertex_base_with_info_2<std::size_t, P_traits> Vb;
+  typedef CGAL::Triangulation_face_base_with_info_2<bool, P_traits>          Fbi;
+  typedef CGAL::Constrained_triangulation_face_base_2<P_traits, Fbi>         Fb;
+  typedef CGAL::Triangulation_data_structure_2<Vb, Fb>                       TDS;
+  // If the polygon is simple, there should be no intersection.
+  typedef CGAL::No_constraint_intersection_tag                               Itag;
+  typedef CGAL::Constrained_Delaunay_triangulation_2<P_traits, TDS, Itag>    CDT;
+  P_traits cdt_traits(avg_normal);
+  CDT cdt(cdt_traits);
+
+  std::vector< std::pair<Point_3, std::size_t> > points_and_ids;
+  points_and_ids.reserve(size);
+  for (std::size_t i = 0; i < size; ++i) {
+    points_and_ids.push_back(std::make_pair(P[i], i));
+  }
+
+  std::vector<typename CDT::Vertex_handle> vertices(size);
+  cdt.insert(points_and_ids.begin(), points_and_ids.end());
+  for (typename CDT::Vertex_handle v : cdt.finite_vertex_handles()) {
+    vertices[v->info()] = v;
+  }
+
+  for (std::size_t i = 0; i < size; ++i) {
+    const std::size_t ip = (i + 1) % size;
+    if (vertices[i] != vertices[ip]) {
+      cdt.insert_constraint(vertices[i], vertices[ip]);
+    }
+  }
+
+  // Mark external faces.
+  for (typename CDT::All_faces_iterator fit = cdt.all_faces_begin(),
+    end = cdt.all_faces_end(); fit != end; ++fit) {
+    fit->info() = false;
+  }
+
+  std::queue<typename CDT::Face_handle> face_queue;
+  face_queue.push(cdt.infinite_vertex()->face());
+  while (!face_queue.empty()) {
+
+    typename CDT::Face_handle fh = face_queue.front();
+    face_queue.pop();
+    if (fh->info()) {
+      continue;
+    }
+
+    fh->info() = true;
+    for (int i = 0; i < 3; ++i) {
+      if (!cdt.is_constrained(typename CDT::Edge(fh, i))) {
+        face_queue.push(fh->neighbor(i));
+      }
+    }
+  }
+
+  if (cdt.dimension() != 2 || cdt.number_of_vertices() != size) {
+    // std::cerr << "WARNING: dim + num vertices, cdt 2 falls back to the original solution!" << std::endl;
+    return false;
+  }
+
+  // Fill the lambda.
+  for (typename CDT::Finite_faces_iterator fit = cdt.finite_faces_begin(),
+    end = cdt.finite_faces_end(); fit != end; ++fit) {
+    if (!fit->info()) { // if it is not external
+
+      std::vector<int> is(3);
+      for (int i = 0; i < 3; ++i) {
+        is[i] = static_cast<int>(fit->vertex(i)->info());
+      }
+
+      std::sort(is.begin(), is.end());
+      lambda.put(is[0], is[2], is[1]);
+      if (!is_valid(P, is[0], is[1], is[2])) {
+        // std::cerr << "WARNING: validity, cdt 2 falls back to the original solution!" << std::endl;
+        return false;
+      }
+    }
+  }
+
+  // Call the tracer. It correctly orients the patch faces.
+  // std::cout << "CDT is being used!" << std::endl;
+  tracer(lambda, 0, static_cast<int>(size) - 1);
+  return true;
+}
+
+#endif
+
+/*******************************************************************************
  * Internal entry point for both polyline and Polyhedron_3 triangulation functions
  ***********************************************************************************/
 template <

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h

@@ -16,6 +16,7 @@
 
 #include <CGAL/disable_warnings.h>
 
+#include <CGAL/bounding_box.h>
 #include <CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polygon_mesh.h>
 #include <CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polyline.h>
 #include <CGAL/Polygon_mesh_processing/refine.h>
@@ -42,10 +43,13 @@ namespace Polygon_mesh_processing {
   /*!
   \ingroup hole_filling_grp
   triangulates a hole in a polygon mesh.
-  The hole must not contain any non-manifold vertex,
-  nor self-intersections.
-  The patch generated does not introduce non-manifold edges nor degenerate triangles.
-  If a hole cannot be triangulated, `pmesh` is not modified and nothing is recorded in `out`.
+
+  Depending on the choice of the underlying algorithm different preconditions apply.
+  When using the 2D constrained Delaunay triangulation, the border edges of the hole
+  must not intersect the surface. Otherwise, additionally, the boundary
+  of the hole must not contain any non-manifold vertex. The patch generated does not
+  introduce non-manifold edges nor degenerate triangles. If a hole cannot be triangulated,
+  `pmesh` is not modified and nothing is recorded in `out`.
 
   @tparam PolygonMesh a model of `MutableFaceGraph`
   @tparam OutputIterator a model of `OutputIterator`
@@ -81,6 +85,27 @@ namespace Polygon_mesh_processing {
       \cgalParamExtra{If no valid triangulation can be found in this search space, the algorithm
                       falls back to the non-Delaunay triangulations search space to find a solution.}
     \cgalParamNEnd
+
+    \cgalParamNBegin{use_2d_constrained_delaunay_triangulation}
+      \cgalParamDescription{If `true`, the points of the boundary of the hole are used
+                            to estimate a fitting plane and a 2D constrained Delaunay triangulation
+                            is then used to fill the hole projected in the fitting plane.}
+      \cgalParamType{Boolean}
+      \cgalParamDefault{`true`}
+      \cgalParamExtra{If the boundary of the hole is not planar (according to the
+                      parameter `threshold_distance`) or if no valid 2D triangulation
+                      can be found, the algorithm falls back to the method using
+                      the 3D Delaunay triangulation. This parameter is a good choice for near planar holes.}
+    \cgalParamNEnd
+
+    \cgalParamNBegin{threshold_distance}
+      \cgalParamDescription{The maximum distance between the vertices of
+                            the hole boundary and the least squares plane fitted to this boundary.}
+      \cgalParamType{double}
+      \cgalParamDefault{one quarter of the height of the bounding box of the hole}
+      \cgalParamExtra{This parameter is used only in conjunction with
+                      the parameter `use_2d_constrained_delaunay_triangulation`.}
+    \cgalParamNEnd
   \cgalNamedParamsEnd
 
   @return `out`
@@ -115,13 +140,45 @@ namespace Polygon_mesh_processing {
 #endif
 
     CGAL_precondition(face(border_halfedge, pmesh) == boost::graph_traits<PolygonMesh>::null_face());
+    bool use_cdt =
+    #ifdef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+        false;
+#else
+        choose_parameter(get_parameter(np, internal_np::use_2d_constrained_delaunay_triangulation), false);
+#endif
 
-    return internal::triangulate_hole_polygon_mesh(pmesh,
+    typename GeomTraits::FT max_squared_distance = typename GeomTraits::FT(-1);
+    if (use_cdt) {
+
+      std::vector<typename GeomTraits::Point_3> points;
+      typedef Halfedge_around_face_circulator<PolygonMesh> Hedge_around_face_circulator;
+      const auto vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), get_property_map(vertex_point, pmesh));
+      Hedge_around_face_circulator circ(border_halfedge, pmesh), done(circ);
+      do {
+        points.push_back(get(vpmap, target(*circ, pmesh)));
+      } while (++circ != done);
+
+      const typename GeomTraits::Iso_cuboid_3 bbox = CGAL::bounding_box(points.begin(), points.end());
+      typename GeomTraits::FT default_squared_distance = CGAL::abs(CGAL::squared_distance(bbox.vertex(0), bbox.vertex(5)));
+      default_squared_distance /= typename GeomTraits::FT(16); // one quarter of the bbox height
+
+      const typename GeomTraits::FT threshold_distance = choose_parameter(
+        get_parameter(np, internal_np::threshold_distance), typename GeomTraits::FT(-1));
+      max_squared_distance = default_squared_distance;
+      if (threshold_distance >= typename GeomTraits::FT(0))
+        max_squared_distance = threshold_distance * threshold_distance;
+      CGAL_assertion(max_squared_distance >= typename GeomTraits::FT(0));
+    }
+
+    return internal::triangulate_hole_polygon_mesh(
+      pmesh,
       border_halfedge,
       out,
       choose_parameter(get_parameter(np, internal_np::vertex_point), get_property_map(vertex_point, pmesh)),
       use_dt3,
-      choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits))).first;
+      choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits)),
+      use_cdt,
+      max_squared_distance).first;
   }
 
   template<typename PolygonMesh, typename OutputIterator>
@@ -192,6 +249,27 @@ namespace Polygon_mesh_processing {
                       falls back to the non-Delaunay triangulations search space to find a solution.}
     \cgalParamNEnd
 
+    \cgalParamNBegin{use_2d_constrained_delaunay_triangulation}
+      \cgalParamDescription{If `true`, the points of the boundary of the hole are used
+                            to estimate a fitting plane and a 2D constrained Delaunay triangulation
+                            is then used to fill the hole projected in the fitting plane.}
+      \cgalParamType{Boolean}
+      \cgalParamDefault{`true`}
+      \cgalParamExtra{If the boundary of the hole is not planar (according to the
+                      parameter `threshold_distance`) or if no valid 2D triangulation
+                      can be found, the algorithm falls back to the method using
+                      the 3D Delaunay triangulation. This parameter is a good choice for near planar holes.}
+    \cgalParamNEnd
+
+    \cgalParamNBegin{threshold_distance}
+      \cgalParamDescription{The maximum distance between the vertices of
+                            the hole boundary and the least squares plane fitted to this boundary.}
+      \cgalParamType{double}
+      \cgalParamDefault{one quarter of the height of the bounding box of the hole}
+      \cgalParamExtra{This parameter is used only in conjunction with
+                      the parameter `use_2d_constrained_delaunay_triangulation`.}
+    \cgalParamNEnd
+
     \cgalParamNBegin{density_control_factor}
       \cgalParamDescription{factor to control density of the ouput mesh,
                             where larger values cause denser refinements, as in `refine()`}
@@ -283,6 +361,27 @@ namespace Polygon_mesh_processing {
                       falls back to the non-Delaunay triangulations search space to find a solution.}
     \cgalParamNEnd
 
+    \cgalParamNBegin{use_2d_constrained_delaunay_triangulation}
+      \cgalParamDescription{If `true`, the points of the boundary of the hole are used
+                            to estimate a fitting plane and a 2D constrained Delaunay triangulation
+                            is then used to fill the hole projected in the fitting plane.}
+      \cgalParamType{Boolean}
+      \cgalParamDefault{`true`}
+      \cgalParamExtra{If the boundary of the hole is not planar (according to the
+                      parameter `threshold_distance`) or if no valid 2D triangulation
+                      can be found, the algorithm falls back to the method using
+                      the 3D Delaunay triangulation. This parameter is a good choice for near planar holes.}
+    \cgalParamNEnd
+
+    \cgalParamNBegin{threshold_distance}
+      \cgalParamDescription{The maximum distance between the vertices of
+                            the hole boundary and the least squares plane fitted to this boundary.}
+      \cgalParamType{double}
+      \cgalParamDefault{one quarter of the height of the bounding box of the hole}
+      \cgalParamExtra{This parameter is used only in conjunction with
+                      the parameter `use_2d_constrained_delaunay_triangulation`.}
+    \cgalParamNEnd
+
     \cgalParamNBegin{density_control_factor}
       \cgalParamDescription{factor to control density of the ouput mesh,
                             where larger values cause denser refinements, as in `refine()`}
@@ -409,6 +508,27 @@ namespace Polygon_mesh_processing {
       \cgalParamExtra{If no valid triangulation can be found in this search space, the algorithm
                       falls back to the non-Delaunay triangulations search space to find a solution.}
     \cgalParamNEnd
+
+    \cgalParamNBegin{use_2d_constrained_delaunay_triangulation}
+      \cgalParamDescription{If `true`, the points of the boundary of the hole are used
+                            to estimate a fitting plane and a 2D constrained Delaunay triangulation
+                            is then used to fill the hole projected in the fitting plane.}
+      \cgalParamType{Boolean}
+      \cgalParamDefault{`true`}
+      \cgalParamExtra{If the boundary of the hole is not planar (according to the
+                      parameter `threshold_distance`) or if no valid 2D triangulation
+                      can be found, the algorithm falls back to the method using
+                      the 3D Delaunay triangulation. This parameter is a good choice for near planar holes.}
+    \cgalParamNEnd
+
+    \cgalParamNBegin{threshold_distance}
+      \cgalParamDescription{The maximum distance between the vertices of
+                            the hole boundary and the least squares plane fitted to this boundary.}
+      \cgalParamType{double}
+      \cgalParamDefault{one quarter of the height of the bounding box of the hole}
+      \cgalParamExtra{This parameter is used only in conjunction with
+                      the parameter `use_2d_constrained_delaunay_triangulation`.}
+    \cgalParamNEnd
   \cgalNamedParamsEnd
 
   \todo handle islands
@@ -426,7 +546,13 @@ namespace Polygon_mesh_processing {
     using parameters::choose_parameter;
     using parameters::get_parameter;
 
-    bool use_dt3 =
+    bool use_cdt =
+#ifdef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+      false;
+#else
+      choose_parameter(get_parameter(np, internal_np::use_2d_constrained_delaunay_triangulation), true);
+#endif
+bool use_dt3 =
 #ifdef CGAL_HOLE_FILLING_DO_NOT_USE_DT3
       false;
 #else
@@ -448,7 +574,32 @@ namespace Polygon_mesh_processing {
     typedef typename PointRange1::iterator InIterator;
     typedef typename std::iterator_traits<InIterator>::value_type Point;
     typedef typename CGAL::Kernel_traits<Point>::Kernel Kernel;
+#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+    struct Always_valid{
+      bool operator()(const std::vector<Point>&, int,int,int)const
+      {return true;}
+    };
+    Always_valid is_valid;
 
+    const typename Kernel::Iso_cuboid_3 bbox = CGAL::bounding_box(points.begin(), points.end());
+    typename Kernel::FT default_squared_distance = CGAL::abs(CGAL::squared_distance(bbox.vertex(0), bbox.vertex(5)));
+    default_squared_distance /= typename Kernel::FT(16); // one quarter of the bbox height
+
+    const typename Kernel::FT threshold_distance = choose_parameter(
+      get_parameter(np, internal_np::threshold_distance), typename Kernel::FT(-1));
+    typename Kernel::FT max_squared_distance = default_squared_distance;
+    if (threshold_distance >= typename Kernel::FT(0))
+      max_squared_distance = threshold_distance * threshold_distance;
+    CGAL_assertion(max_squared_distance >= typename Kernel::FT(0));
+
+    if(!use_cdt ||
+       !triangulate_hole_polyline_with_cdt(
+         points,
+         tracer,
+         is_valid,
+         choose_parameter<Kernel>(get_parameter(np, internal_np::geom_traits)),
+         max_squared_distance))
+#endif
     triangulate_hole_polyline(points, third_points, tracer, WC(),
                               use_dt3,
                               choose_parameter<Kernel>(get_parameter(np, internal_np::geom_traits)));

Polygon_mesh_processing/package_info/Polygon_mesh_processing/dependencies

@@ -21,6 +21,7 @@ Modular_arithmetic
 Number_types
 Polygon
 Polygon_mesh_processing
+Principal_component_analysis_LGPL
 Profiling_tools
 Property_map
 Random_numbers

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -73,8 +73,7 @@ create_single_source_cgal_program("remeshing_test.cpp")
 create_single_source_cgal_program("measures_test.cpp")
 create_single_source_cgal_program("triangulate_faces_test.cpp")
 create_single_source_cgal_program("triangulate_faces_hole_filling_dt3_test.cpp")
-create_single_source_cgal_program(
-  "triangulate_faces_hole_filling_all_search_test.cpp")
+create_single_source_cgal_program("triangulate_faces_hole_filling_all_search_test.cpp")
 create_single_source_cgal_program("test_pmp_remove_border_edge.cpp")
 create_single_source_cgal_program("test_pmp_distance.cpp")
 create_single_source_cgal_program("test_corefinement_and_constraints.cpp")
@@ -104,8 +103,8 @@ create_single_source_cgal_program("test_pmp_manifoldness.cpp")
 create_single_source_cgal_program("test_mesh_smoothing.cpp")
 create_single_source_cgal_program("test_remove_caps_needles.cpp")
 create_single_source_cgal_program("test_simplify_polylines_pmp.cpp")
+create_single_source_cgal_program("triangulate_hole_with_cdt_2_test.cpp")
 create_single_source_cgal_program("test_pmp_polyhedral_envelope.cpp")
-
 #  create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp")
 
 if(TARGET CGAL::TBB_support)

Polygon_mesh_processing/test/Polygon_mesh_processing/data/w_horizontal_hole.off

@@ -0,0 +1,110 @@
+OFF
+54 54 0
+519 0 100
+734 -817 0
+734 -817 100
+775 -636 0
+775 -636 100
+944 0 0
+944 0 100
+1131 0 0
+1131 0 100
+1463 -1062 0
+1463 -1062 100
+1288 -1062 0
+1288 -1062 100
+1106 -448 0
+1106 -448 100
+1045 -243 0
+1045 -243 100
+992 -446 0
+992 -446 100
+833 -1062 0
+833 -1062 100
+648 -1062 0
+648 -1062 100
+479 -440 0
+479 -440 100
+468 -400 0
+468 -400 100
+459 -366 0
+459 -366 100
+451 -334 0
+451 -334 100
+445 -307 0
+445 -307 100
+439 -284 0
+439 -284 100
+434 -263 0
+434 -263 100
+430 -247 0
+430 -247 100
+426 -234 0
+426 -234 100
+424 -225 0
+424 -225 100
+424 -221 0
+424 -221 100
+361 -449 0
+361 -449 100
+192 -1062 0
+192 -1062 100
+6 -1062 0
+6 -1062 100
+331 0 0
+331 0 100
+519 0 0
+3 2 1 0
+3 53 0 1
+3 4 3 2
+3 1 2 3
+3 6 5 4
+3 3 4 5
+3 8 7 6
+3 5 6 7
+3 10 9 8
+3 7 8 9
+3 12 11 10
+3 9 10 11
+3 14 13 12
+3 11 12 13
+3 16 15 14
+3 13 14 15
+3 18 17 16
+3 15 16 17
+3 20 19 18
+3 17 18 19
+3 22 21 20
+3 19 20 21
+3 24 23 22
+3 21 22 23
+3 26 25 24
+3 23 24 25
+3 28 27 26
+3 25 26 27
+3 30 29 28
+3 27 28 29
+3 32 31 30
+3 29 30 31
+3 34 33 32
+3 31 32 33
+3 36 35 34
+3 33 34 35
+3 38 37 36
+3 35 36 37
+3 40 39 38
+3 37 38 39
+3 42 41 40
+3 39 40 41
+3 44 43 42
+3 41 42 43
+3 46 45 44
+3 43 44 45
+3 48 47 46
+3 45 46 47
+3 50 49 48
+3 47 48 49
+3 52 51 50
+3 49 50 51
+3 0 53 52
+3 51 52 53

Polygon_mesh_processing/test/Polygon_mesh_processing/data/w_orthogonal_hole.off

@@ -0,0 +1,156 @@
+OFF
+54 100 0
+519 0 100 
+734 -817 0 
+734 -817 100 
+775 -636 0 
+775 -636 100 
+944 0 0 
+944 0 100 
+1131 0 0 
+1131 0 100 
+1463 -1062 0 
+1463 -1062 100 
+1288 -1062 0 
+1288 -1062 100 
+1106 -448 0 
+1106 -448 100 
+1045 -243 0 
+1045 -243 100 
+992 -446 0 
+992 -446 100 
+833 -1062 0 
+833 -1062 100 
+648 -1062 0 
+648 -1062 100 
+479 -440 0 
+479 -440 100 
+468 -400 0 
+468 -400 100 
+459 -366 0 
+459 -366 100 
+451 -334 0 
+451 -334 100 
+445 -307 0 
+445 -307 100 
+439 -284 0 
+439 -284 100 
+434 -263 0 
+434 -263 100 
+430 -247 0 
+430 -247 100 
+426 -234 0 
+426 -234 100 
+424 -225 0 
+424 -225 100 
+424 -221 0 
+424 -221 100 
+361 -449 0 
+361 -449 100 
+192 -1062 0 
+192 -1062 100 
+6 -1062 0 
+6 -1062 100 
+331 0 0 
+331 0 100 
+519 0 0 
+3 2 1 0
+3 53 0 1
+3 4 3 2
+3 1 2 3
+3 6 5 4
+3 3 4 5
+3 8 7 6
+3 5 6 7
+3 12 11 10
+3 9 10 11
+3 16 15 14
+3 13 14 15
+3 18 17 16
+3 15 16 17
+3 20 19 18
+3 17 18 19
+3 22 21 20
+3 19 20 21
+3 24 23 22
+3 21 22 23
+3 26 25 24
+3 23 24 25
+3 28 27 26
+3 25 26 27
+3 30 29 28
+3 27 28 29
+3 32 31 30
+3 29 30 31
+3 34 33 32
+3 31 32 33
+3 36 35 34
+3 33 34 35
+3 38 37 36
+3 35 36 37
+3 40 39 38
+3 37 38 39
+3 42 41 40
+3 39 40 41
+3 44 43 42
+3 41 42 43
+3 46 45 44
+3 43 44 45
+3 48 47 46
+3 45 46 47
+3 50 49 48
+3 47 48 49
+3 52 51 50
+3 49 50 51
+3 0 53 52
+3 51 52 53
+3 50 48 46
+3 36 34 0
+3 52 44 0
+3 52 50 46
+3 38 36 0
+3 52 46 44
+3 44 42 40
+3 44 40 38
+3 44 38 0
+3 24 2 0
+3 26 24 0
+3 28 26 0
+3 30 28 0
+3 18 4 20
+3 32 30 0
+3 2 22 20
+3 14 10 8
+3 12 10 14
+3 14 8 16
+3 16 8 6
+3 24 22 2
+3 4 2 20
+3 16 6 18
+3 18 6 4
+3 34 32 0
+3 21 23 1
+3 11 13 9
+3 19 21 1
+3 17 5 15
+3 13 7 9
+3 5 7 15
+3 15 7 13
+3 3 5 17
+3 3 17 19
+3 1 3 19
+3 37 53 35
+3 35 53 33
+3 39 41 43
+3 39 43 37
+3 51 53 43
+3 43 53 37
+3 49 51 45
+3 49 45 47
+3 51 43 45
+3 33 53 31
+3 31 53 29
+3 29 53 27
+3 27 53 25
+3 25 53 23
+3 23 53 1

Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp

@@ -95,7 +95,7 @@ CGAL::internal::Weight_min_max_dihedral_and_area
 }
 
 
-void test_triangulate_hole_weight(const char* file_name, bool use_DT, std::size_t nb_remaining_holes) {
+void test_triangulate_hole_weight(const char* file_name, bool use_DT, std::size_t nb_remaining_holes) { //don't test with cdt as we are testing the weights and there is no weight in the cdt version
   typedef CGAL::internal::Weight_min_max_dihedral_and_area Weight;
 
   std::cout << "test_triangulate_hole_weight + useDT: " << use_DT << std::endl;
@@ -107,7 +107,7 @@ void test_triangulate_hole_weight(const char* file_name, bool use_DT, std::size_
   for(std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it) {
     std::vector<Facet_handle> patch;
     Weight w_algo = CGAL::Polygon_mesh_processing::internal::triangulate_hole_polygon_mesh(
-      poly, *it, back_inserter(patch), get(CGAL::vertex_point, poly), use_DT, Kernel()).second;
+      poly, *it, back_inserter(patch), get(CGAL::vertex_point, poly), use_DT, Kernel(), false, 0).second;
     if(patch.empty()) { continue; }
     Weight w_test = calculate_weight_for_patch(poly, patch.begin(), patch.end());
 
@@ -127,7 +127,7 @@ void test_triangulate_hole_weight(const char* file_name, bool use_DT, std::size_
 }
 /******************************************************************/
 
-void test_triangulate_hole(const char* file_name) {
+void test_triangulate_hole(const char* file_name, bool use_cdt) {
   std::cout << "test_triangulate_hole:" << std::endl;
   std::cout << "  File: "<< file_name  << std::endl;
   Polyhedron poly;
@@ -136,7 +136,8 @@ void test_triangulate_hole(const char* file_name) {
 
   for(std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it) {
     std::vector<Facet_handle> patch;
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch));
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch),
+                                                    CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
     if(patch.empty()) {
       std::cerr << "  Error: empty patch created." << std::endl;
       assert(false);
@@ -151,7 +152,7 @@ void test_triangulate_hole(const char* file_name) {
   std::cout << "  Done!" << std::endl;
 }
 
-void test_triangulate_hole_should_be_no_output(const char* file_name) {
+void test_triangulate_hole_should_be_no_output(const char* file_name, bool use_cdt) {
   std::cout << "test_triangulate_hole_should_be_no_output:" << std::endl;
   std::cout << "  File: "<< file_name  << std::endl;
   Polyhedron poly;
@@ -161,7 +162,8 @@ void test_triangulate_hole_should_be_no_output(const char* file_name) {
   for(std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it) {
     std::vector<Facet_handle> patch;
     CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch),
-      CGAL::Polygon_mesh_processing::parameters::use_delaunay_triangulation(false));
+      CGAL::Polygon_mesh_processing::parameters::use_delaunay_triangulation(false)
+        .use_2d_constrained_delaunay_triangulation(use_cdt));
     if(!patch.empty()) {
       std::cerr << "  Error: patch should be empty" << std::endl;
       assert(false);
@@ -178,7 +180,7 @@ void test_triangulate_hole_should_be_no_output(const char* file_name) {
   std::cout << "  Done!" << std::endl;
 }
 
-void test_triangulate_and_refine_hole(const char* file_name) {
+void test_triangulate_and_refine_hole(const char* file_name, bool use_cdt) {
   std::cout << "test_triangulate_and_refine_hole:" << std::endl;
   std::cout << "  File: "<< file_name  << std::endl;
   Polyhedron poly;
@@ -189,7 +191,8 @@ void test_triangulate_and_refine_hole(const char* file_name) {
     std::vector<Facet_handle> patch_facets;
     std::vector<Vertex_handle> patch_vertices;
     CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly, *it,
-      back_inserter(patch_facets), back_inserter(patch_vertices));
+      back_inserter(patch_facets), back_inserter(patch_vertices),
+                                                               CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
 
     if(patch_facets.empty()) {
       std::cerr << "  Error: empty patch created." << std::endl;
@@ -205,7 +208,7 @@ void test_triangulate_and_refine_hole(const char* file_name) {
   std::cout << "  Done!" << std::endl;
 }
 
-void test_triangulate_refine_and_fair_hole(const char* file_name) {
+void test_triangulate_refine_and_fair_hole(const char* file_name, bool use_cdt) {
   std::cout << "test_triangulate_refine_and_fair_hole:" << std::endl;
   std::cout << "  File: "<< file_name  << std::endl;
   Polyhedron poly;
@@ -216,7 +219,8 @@ void test_triangulate_refine_and_fair_hole(const char* file_name) {
     std::vector<Facet_handle> patch_facets;
     std::vector<Vertex_handle> patch_vertices;
     CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole(poly,
-      *it, back_inserter(patch_facets), back_inserter(patch_vertices));
+      *it, back_inserter(patch_facets), back_inserter(patch_vertices),
+      CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
 
     if(patch_facets.empty()) {
       std::cerr << "  Error: empty patch created." << std::endl;
@@ -232,7 +236,7 @@ void test_triangulate_refine_and_fair_hole(const char* file_name) {
   std::cout << "  Done!" << std::endl;
 }
 
-void test_ouput_iterators_triangulate_hole(const char* file_name) {
+void test_ouput_iterators_triangulate_hole(const char* file_name, bool use_cdt) {
   std::cout << "test_ouput_iterators_triangulate_hole:" << std::endl;
   std::cout << "  File: "<< file_name  << std::endl;
 
@@ -245,7 +249,8 @@ void test_ouput_iterators_triangulate_hole(const char* file_name) {
   std::vector<Halfedge_handle>::iterator it_2 = border_reps_2.begin();
   for(std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it, ++it_2) {
     std::vector<Facet_handle> patch;
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch));
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch),
+              CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
 
     std::vector<Facet_handle> patch_2 = patch;
     Facet_handle* output_it =
@@ -254,13 +259,14 @@ void test_ouput_iterators_triangulate_hole(const char* file_name) {
     if(patch.size() != (std::size_t)(output_it - &*patch_2.begin())) {
       std::cerr << "  Error: returned facet output iterator is not valid!" << std::endl;
       std::cerr << "  " << patch.size() << " vs " << (output_it - &*patch_2.begin()) << std::endl;
+
       assert(false);
     }
   }
   std::cout << "  Done!" << std::endl;
 }
 
-void test_ouput_iterators_triangulate_and_refine_hole(const char* file_name) {
+void test_ouput_iterators_triangulate_and_refine_hole(const char* file_name, bool use_cdt) {
   std::cout << "test_ouput_iterators_triangulate_and_refine_hole:" << std::endl;
   std::cout << "  File: "<< file_name  << std::endl;
 
@@ -275,7 +281,8 @@ void test_ouput_iterators_triangulate_and_refine_hole(const char* file_name) {
     std::vector<Facet_handle> patch_facets;
     std::vector<Vertex_handle> patch_vertices;
     CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly,
-      *it, back_inserter(patch_facets), back_inserter(patch_vertices));
+      *it, back_inserter(patch_facets), back_inserter(patch_vertices),
+      CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
     // create enough space to hold outputs
     std::vector<Facet_handle> patch_facets_2 = patch_facets;
     std::vector<Vertex_handle> patch_vertices_2 = patch_vertices;
@@ -283,7 +290,8 @@ void test_ouput_iterators_triangulate_and_refine_hole(const char* file_name) {
 
     std::pair<Facet_handle*, Vertex_handle*> output_its =
       CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly_2,
-        *it_2, &*patch_facets_2.begin(), &*patch_vertices_2.begin());
+        *it_2, &*patch_facets_2.begin(), &*patch_vertices_2.begin(),
+        CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
 
     if(patch_facets.size() != (std::size_t) (output_its.first - &*patch_facets_2.begin())) {
       std::cout << "  Error: returned facet output iterator is not valid!" << std::endl;
@@ -320,7 +328,7 @@ void test_triangulate_refine_and_fair_hole_compile() {
   (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices),
   CGAL::Polygon_mesh_processing::parameters::weight_calculator(
     CGAL::internal::Uniform_weight_fairing<Polyhedron>(poly)).
-    sparse_linear_solver(Default_solver()));
+    sparse_linear_solver(Default_solver()).use_2d_constrained_delaunay_triangulation(false));
 
   // default solver
   read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
@@ -361,26 +369,30 @@ int main()
   std::cerr.precision(17);
 
   generate_elephant_with_hole();
-
   std::vector<std::string> input_files;
   input_files.push_back("elephant_triangle_hole.off");
   input_files.push_back("elephant_quad_hole.off");
   input_files.push_back("data/mech-holes-shark.off");
-
-  for(std::vector<std::string>::iterator it = input_files.begin(); it != input_files.end(); ++it)
-  {
-    test_triangulate_hole(it->c_str());
-    test_triangulate_and_refine_hole(it->c_str());
-    test_triangulate_refine_and_fair_hole(it->c_str());
-    test_ouput_iterators_triangulate_and_refine_hole(it->c_str());
-    test_ouput_iterators_triangulate_hole(it->c_str());
+  // std::cerr.precision(15);
+  for(std::vector<std::string>::iterator it = input_files.begin(); it != input_files.end(); ++it) {
+    test_triangulate_hole(it->c_str(), true);
+    test_triangulate_hole(it->c_str(), false);
+    test_triangulate_and_refine_hole(it->c_str(), true);
+    test_triangulate_and_refine_hole(it->c_str(), false);
+    test_triangulate_refine_and_fair_hole(it->c_str(), true);
+    test_triangulate_refine_and_fair_hole(it->c_str(), false);
+    test_ouput_iterators_triangulate_and_refine_hole(it->c_str(), true);
+    test_ouput_iterators_triangulate_and_refine_hole(it->c_str(), false);
+    test_ouput_iterators_triangulate_hole(it->c_str(), true);
+    test_ouput_iterators_triangulate_hole(it->c_str(), false);
     test_triangulate_hole_weight(it->c_str(), true, 0);
     test_triangulate_hole_weight(it->c_str(), false, 0);
     std::cout << "------------------------------------------------" << std::endl;
   }
-
-  test_triangulate_hole_should_be_no_output("data/non_manifold_vertex.off");
-  test_triangulate_hole_should_be_no_output("data/two_tris_collinear.off");
+  test_triangulate_hole_should_be_no_output("data/non_manifold_vertex.off", true);
+  test_triangulate_hole_should_be_no_output("data/non_manifold_vertex.off", false);
+  test_triangulate_hole_should_be_no_output("data/two_tris_collinear.off", true);
+  test_triangulate_hole_should_be_no_output("data/two_tris_collinear.off", false);
 
   test_triangulate_refine_and_fair_hole_compile();
   std::cout << "All Done!" << std::endl;

Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_polyline_test.cpp

@@ -150,7 +150,8 @@ void test_1(const char* file_name, bool use_DT, bool save_output) {
   std::vector<boost::tuple<int, int, int> > tris;
   CGAL::Polygon_mesh_processing::triangulate_hole_polyline(
     points, std::back_inserter(tris),
-    CGAL::Polygon_mesh_processing::parameters::use_delaunay_triangulation(use_DT));
+    CGAL::Polygon_mesh_processing::parameters::use_delaunay_triangulation(use_DT)
+        .use_2d_constrained_delaunay_triangulation(true));
 
   check_triangles(points, tris);
   check_constructed_polyhedron(file_name, &tris, &points, save_output);
@@ -168,7 +169,8 @@ void test_2(const char* file_name, bool use_DT, bool save_output) {
   std::vector<boost::tuple<int, int, int> > tris;
   CGAL::Polygon_mesh_processing::triangulate_hole_polyline(
     points, extras, std::back_inserter(tris),
-    CGAL::Polygon_mesh_processing::parameters::use_delaunay_triangulation(use_DT));
+    CGAL::Polygon_mesh_processing::parameters::use_delaunay_triangulation(use_DT)
+        .use_2d_constrained_delaunay_triangulation(true));
 
   check_triangles(points, tris);
   check_constructed_polyhedron(file_name, &tris, &points, save_output);
@@ -185,7 +187,8 @@ void test_should_be_no_output(const char* file_name, bool use_DT) {
   std::vector<boost::tuple<int, int, int> > tris;
   CGAL::Polygon_mesh_processing::triangulate_hole_polyline(
     points, std::back_inserter(tris),
-    CGAL::Polygon_mesh_processing::parameters::use_delaunay_triangulation(use_DT));
+    CGAL::Polygon_mesh_processing::parameters::use_delaunay_triangulation(use_DT)
+        .use_2d_constrained_delaunay_triangulation(true));
 
   if(!tris.empty()) {
     std::cerr << "  Error: patch should be empty" << std::endl;

Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_with_cdt_2_test.cpp

@@ -0,0 +1,175 @@
+#include <set>
+#include <vector>
+#include <fstream>
+#include <cassert>
+#include <string>
+
+#define CGAL_NO_CDT_2_WARNING
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
+#include <CGAL/Polygon_mesh_processing/orientation.h>
+
+template<
+class PolygonMesh,
+class Halfedge_handle>
+void detect_borders(
+  PolygonMesh& pmesh,
+  std::vector<Halfedge_handle>& borders) {
+
+  typedef CGAL::Halfedge_around_face_circulator<PolygonMesh>
+    Halfedge_around_facet_circulator;
+
+  borders.clear();
+  std::set<Halfedge_handle> border_map;
+  for (Halfedge_handle h : halfedges(pmesh)) {
+    if (
+      face(h, pmesh) == boost::graph_traits<PolygonMesh>::null_face() &&
+      border_map.find(h) == border_map.end()) {
+
+      borders.push_back(h);
+      Halfedge_around_facet_circulator hf_around_facet(h, pmesh);
+      Halfedge_around_facet_circulator done(hf_around_facet);
+
+      do {
+        assert(border_map.insert(*hf_around_facet).second); // is insertion ok?
+      } while (++hf_around_facet != done);
+    }
+  }
+}
+
+// This test is inspired by the issue: https://github.com/CGAL/cgal/issues/4464.
+template<
+typename PolygonMesh,
+typename GeomTraits>
+void test_triangulate_hole_with_cdt_2(
+  const std::string kernel_name,
+  int argc, char **argv,
+  const std::string file_name,
+  const std::size_t num_borders,
+  const std::size_t num_patch_faces,
+  const bool verbose) {
+
+  typedef typename boost::graph_traits<PolygonMesh>::face_descriptor Face_handle;
+  typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor Halfedge_handle;
+
+  // Reading the file.
+  if (verbose) {
+    std::cout << "test with the " << kernel_name << " kernel:" << std::endl;
+  }
+  PolygonMesh pmesh;
+  std::string path = "data/" + file_name + ".off";
+  std::ifstream in(path.c_str(), std::ios_base::in);
+  CGAL::set_ascii_mode(in);
+  CGAL::read_OFF(in, pmesh);
+  in.close();
+  if (verbose) {
+    std::cout << "* finished reading the file" << std::endl;
+  }
+
+  // Detecting the hole borders.
+  std::vector<Halfedge_handle> borders;
+  detect_borders(pmesh, borders);
+  if (verbose) {
+    std::cout << "* number of detected borders: " <<
+      borders.size() << std::endl;
+  }
+  assert(borders.size() == num_borders);
+
+  // Triangulating the holes.
+  std::vector<Face_handle> patch_faces;
+  for (const Halfedge_handle& h : borders) {
+    patch_faces.clear();
+    CGAL::Polygon_mesh_processing::triangulate_hole(
+      pmesh,
+      h,
+      std::back_inserter(patch_faces),
+      CGAL::Polygon_mesh_processing::parameters::vertex_point_map(
+        get(CGAL::vertex_point, pmesh)).
+        use_2d_constrained_delaunay_triangulation(true).
+        geom_traits(GeomTraits()));
+
+    if (verbose) {
+      std::cout << "* number of faces in the constructed patch: " <<
+        patch_faces.size() << std::endl;
+    }
+    assert(patch_faces.size() == num_patch_faces);
+  }
+  assert(pmesh.is_valid() && is_closed(pmesh));
+  assert(CGAL::Polygon_mesh_processing::is_outward_oriented(pmesh,
+    CGAL::parameters::all_default()));
+
+  // Writing the file.
+  if (verbose) {
+    path = "";
+    if (argc > 1) path = std::string(argv[1]);
+    path += "4464_" + file_name + ".off";
+    std::ofstream out(path.c_str(), std::ios_base::out);
+    CGAL::set_ascii_mode(out);
+    CGAL::write_OFF(out, pmesh);
+    out.close();
+    std::cout << "* finished writing the file" << std::endl;
+  }
+}
+
+int main(int argc, char **argv) {
+
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel EI;
+  typedef CGAL::Exact_predicates_exact_constructions_kernel EE;
+
+  typedef CGAL::Surface_mesh<EI::Point_3> Surface_mesh_EI;
+  typedef CGAL::Polyhedron_3<EE> Polyhedron_3_EE;
+
+  // Checking on a data file with two planar, simple, and horizontal holes.
+  test_triangulate_hole_with_cdt_2<Surface_mesh_EI, EI>(
+    "exact_inexact", argc, argv, "w_horizontal_hole", 2, 25, false);
+  test_triangulate_hole_with_cdt_2<Polyhedron_3_EE, EE>(
+    "exact_exact", argc, argv, "w_horizontal_hole", 2, 25, false);
+  std::cout <<
+    "test_triangulate_hole_with_cdt_2: horizontal planar hole SUCCESS" << std::endl;
+
+  // Checking on a data file with two planar, simple, and orthogonal holes.
+  test_triangulate_hole_with_cdt_2<Surface_mesh_EI, EI>(
+    "exact_inexact", argc, argv, "w_orthogonal_hole", 2, 2, false);
+  test_triangulate_hole_with_cdt_2<Polyhedron_3_EE, EE>(
+    "exact_exact", argc, argv, "w_orthogonal_hole", 2, 2, false);
+  std::cout <<
+    "test_triangulate_hole_with_cdt_2: orthogonal planar hole SUCCESS" << std::endl;
+
+  // Checking on a data file with two planar, simple, and horizontal holes.
+  test_triangulate_hole_with_cdt_2<Surface_mesh_EI, EI>(
+    "exact_inexact", argc, argv, "elephant_flat_hole", 1, 17, false);
+  test_triangulate_hole_with_cdt_2<Polyhedron_3_EE, EE>(
+    "exact_exact", argc, argv, "elephant_flat_hole", 1, 17, false);
+  std::cout <<
+    "test_triangulate_hole_with_cdt_2: near flat hole SUCCESS" << std::endl;
+
+  // Checking on a data file with two planar, simple, and horizontal holes.
+  test_triangulate_hole_with_cdt_2<Surface_mesh_EI, EI>(
+    "exact_inexact", argc, argv, "elephant_concave_hole", 1, 24, false);
+  test_triangulate_hole_with_cdt_2<Polyhedron_3_EE, EE>(
+    "exact_exact", argc, argv, "elephant_concave_hole", 1, 24, false);
+  std::cout <<
+    "test_triangulate_hole_with_cdt_2: concave hole SUCCESS" << std::endl;
+
+  // Checking on a data file with one simple but not a planar hole.
+  test_triangulate_hole_with_cdt_2<Surface_mesh_EI, EI>(
+    "exact_inexact", argc, argv, "elephant_curved_hole", 1, 19, false);
+  test_triangulate_hole_with_cdt_2<Polyhedron_3_EE, EE>(
+    "exact_exact", argc, argv, "elephant_curved_hole", 1, 19, false);
+  std::cout <<
+    "test_triangulate_hole_with_cdt_2: curved hole SUCCESS" << std::endl;
+
+  // Checking on a data file with one hole that is neither simple nor planar.
+  test_triangulate_hole_with_cdt_2<Surface_mesh_EI, EI>(
+    "exact_inexact", argc, argv, "elephant_complex_hole", 1, 29, false);
+  test_triangulate_hole_with_cdt_2<Polyhedron_3_EE, EE>(
+    "exact_exact", argc, argv, "elephant_complex_hole", 1, 29, false);
+  std::cout <<
+    "test_triangulate_hole_with_cdt_2: complex hole SUCCESS" << std::endl;
+
+  return EXIT_SUCCESS;
+}

Triangulation_2/include/CGAL/Triangulation_2/internal/Triangulation_2_projection_traits_base_3.h

@@ -289,6 +289,33 @@ public:
   }
 }; // end class Compare_along_axis
 
+template <class Traits>
+class Less_xy_along_axis
+{
+  // private members
+  typedef typename Traits::Vector_3 Vector_3;
+  typedef typename Traits::Point_2 Point;
+  Vector_3 base1, base2;
+public:
+  Less_xy_along_axis(const Vector_3& base1, const Vector_3& base2) : base1(base1), base2(base2)
+  {
+    CGAL_PROFILER("Construct Less_xy_along_axis")
+    CGAL_TIME_PROFILER("Construct Less_xy_along_axis")
+  }
+
+  typedef bool result_type;
+
+  bool operator() (const Point &p, const Point &q) const {
+
+    Compare_along_axis<Traits> cx(base1);
+    Comparison_result crx = cx(p, q);
+    if (crx == SMALLER) { return true; }
+    if (crx == LARGER) { return false; }
+    Less_along_axis<Traits> ly(base2);
+    return ly(p, q);
+  }
+}; // end class Less_xy_along_axis
+
 } // end namespace TriangulationProjectionTraitsCartesianFunctors
 
 
@@ -353,6 +380,8 @@ public:
     Less_along_axis<Self>                                    Less_x_2;
   typedef TriangulationProjectionTraitsCartesianFunctors::
     Less_along_axis<Self>                                    Less_y_2;
+  typedef TriangulationProjectionTraitsCartesianFunctors::
+    Less_xy_along_axis<Self>                                 Less_xy_2;
 
   typedef TriangulationProjectionTraitsCartesianFunctors::
     Projected_orientation_with_normal_3<Self>                Orientation_2;
@@ -393,6 +422,12 @@ public:
     return Less_y_2(this->base2());
   }
 
+  Less_xy_2
+  less_xy_2_object() const
+  {
+    return Less_xy_2(this->base1(), this->base2());
+  }
+
   Compare_x_2
   compare_x_2_object() const
   {