Merge remote-tracking branch 'cgal/master' into SLS-Weighted_skeleton-GF

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

@@ -110,22 +110,75 @@ public:
   typedef typename boost::property_traits<PMap>::value_type type;
 };
 
+
+template <typename PolygonMesh,
+          typename VPM_from_NP>
+struct GetVertexPointMap_impl
+{
+  typedef VPM_from_NP type;
+  typedef VPM_from_NP const_type;
+
+  template<class NamedParameters>
+  static const_type
+  get_const_map(const NamedParameters& np, const PolygonMesh&)
+  {
+    return parameters::get_parameter(np, internal_np::vertex_point);
+  }
+
+  template<class NamedParameters>
+  static type
+  get_map(const NamedParameters& np, PolygonMesh&)
+  {
+    return parameters::get_parameter(np, internal_np::vertex_point);
+  }
+};
+
+template <typename PolygonMesh>
+struct GetVertexPointMap_impl<PolygonMesh, internal_np::Param_not_found>
+{
+  typedef typename property_map_selector<PolygonMesh, boost::vertex_point_t>::const_type const_type;
+  typedef typename property_map_selector<PolygonMesh, boost::vertex_point_t>::type type;
+
+  template<class NamedParameters>
+  static const_type
+  get_const_map(const NamedParameters& /* np */, const PolygonMesh& pm)
+  {
+    return get_const_property_map(boost::vertex_point, pm);
+  }
+
+  template<class NamedParameters>
+  static type
+  get_map(const NamedParameters& /* np */, PolygonMesh& pm)
+  {
+    return get_property_map(boost::vertex_point, pm);
+  }
+};
+
 template <typename PolygonMesh,
           typename NamedParameters = parameters::Default_named_parameters>
 class GetVertexPointMap
 {
-  typedef typename property_map_selector<PolygonMesh, boost::vertex_point_t>::const_type
+  typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_point_t,
-  DefaultVPMap_const;
+                                                       NamedParameters,
-  typedef typename property_map_selector<PolygonMesh, boost::vertex_point_t>::type
+                                                       internal_np::Param_not_found>::type VPM_from_NP;
-  DefaultVPMap;
+
+  typedef GetVertexPointMap_impl<PolygonMesh, VPM_from_NP> Impl;
 
 public:
-  typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_point_t,
+  typedef typename Impl::type type;
-                                                       NamedParameters,
+  typedef typename Impl::const_type const_type;
-                                                       DefaultVPMap>::type type;
+
-  typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_point_t,
+  static const_type
-                                                       NamedParameters,
+  get_const_map(const NamedParameters& np, const PolygonMesh& pm)
-                                                       DefaultVPMap_const>::type const_type;
+  {
+    return Impl::get_const_map(np, pm);
+  }
+
+  static type
+  get_map(const NamedParameters& np, PolygonMesh& pm)
+  {
+    return Impl::get_map(np, pm);
+  }
 };
 
 template<typename PolygonMesh, typename NamedParameters>
@@ -138,10 +191,15 @@ public:
   typedef typename CGAL::Kernel_traits<Point>::Kernel Kernel;
 };
 
-template <typename PolygonMesh,
+
-          typename NamedParametersGT = parameters::Default_named_parameters,
+template<typename PolygonMesh, class GT, class NamedParametersVPM>
-          typename NamedParametersVPM = NamedParametersGT>
+struct GetGeomTraits_impl
-class GetGeomTraits
+{
+  typedef GT type;
+};
+
+template<typename PolygonMesh, class NamedParametersVPM>
+struct GetGeomTraits_impl<PolygonMesh, internal_np::Param_not_found, NamedParametersVPM>
 {
   typedef typename CGAL::graph_has_property<PolygonMesh, boost::vertex_point_t>::type Has_internal_pmap;
 
@@ -154,12 +212,20 @@ class GetGeomTraits
   typedef typename boost::mpl::if_c<Has_internal_pmap::value ||
                                     !std::is_same<internal_np::Param_not_found, NP_vpm>::value,
                                     typename GetK<PolygonMesh, NamedParametersVPM>::Kernel,
-                                    Fake_GT>::type DefaultKernel;
+                                    Fake_GT>::type type;
+};
 
-public:
+template <typename PolygonMesh,
+          typename NamedParametersGT = parameters::Default_named_parameters,
+          typename NamedParametersVPM = NamedParametersGT>
+struct GetGeomTraits
+{
   typedef typename internal_np::Lookup_named_param_def<internal_np::geom_traits_t,
                                                        NamedParametersGT,
-                                                       DefaultKernel>::type type;
+                                                       internal_np::Param_not_found>::type GT_from_NP;
+  typedef typename GetGeomTraits_impl<PolygonMesh,
+                                          GT_from_NP,
+                                          NamedParametersVPM>::type type;
 };
 
 // Define the following structs:
@@ -278,6 +344,21 @@ public:
   typedef typename CGAL::Identity_property_map<const Dummy_point> const_type;
 };
 
+template <typename PointRange, typename NamedParameters>
+struct GetPolygonSoupGeomTraits
+{
+  typedef typename internal_np::Lookup_named_param_def <
+                     internal_np::geom_traits_t,
+                     NamedParameters,
+                     typename CGAL::Kernel_traits<
+                       typename boost::property_traits<
+                         typename GetPointMap<PointRange, NamedParameters>::type
+                       >::value_type
+                     >::type
+                   > ::type                                                         type;
+};
+
+
 template <class PointRange, class NamedParameters, class PointMap = Default, class NormalMap = Default>
 struct Point_set_processing_3_np_helper
 {

Installation/CHANGES.md

@@ -40,11 +40,15 @@ CGAL tetrahedral Delaunay refinement algorithm.
 
 -   Added the function `CGAL::Polygon_mesh_processing::remove_almost_degenerate_faces()` to remove badly shaped triangles faces in a mesh.
 
+-   Added the function `CGAL::Polygon_mesh_processing::triangulate_polygons()`, which allows users to triangulate polygon soups.
+
 -   Added the functions `CGAL::Polygon_mesh_processing::remesh_planar_patches()` and
     `CGAL::Polygon_mesh_processing::remesh_almost_coplanar_patches()` to retriangulate patches of coplanar faces in a mesh.
 
 -   Added a named parameter to `CGAL::Polygon_mesh_processing::smooth_shape()` to disable scaling to compensate volume loss.
 
+-   Added the functions `CGAL::Polygon_mesh_processing::does_triangle_soup_self_intersect()` and `CGAL::Polygon_mesh_processing::triangle_soup_self_intersections()` to identify and report self-intersections in a triangle soup, similarly to existing functions on triangle meshes.
+
 ### [3D Simplicial Mesh Data Structure](https://doc.cgal.org/5.6/Manual/packages.html#PkgSMDS3) (new package)
 
 -   This new package wraps all the existing code that deals with a `MeshComplex_3InTriangulation_3` to describe 3D simplicial meshes, and makes the data structure independent from the tetrahedral mesh generation package.
@@ -117,6 +121,7 @@ extracted from labeled images.
 
 ### [2D Minkowski Sums](https://doc.cgal.org/5.6/Manual/packages.html#PkgMinkowskiSum2)
 -   Fixed a bug that made holes in the Minkowski sum disappear
+-   Fixed hole filtering. The code used to erroneously remove holes from the container of holes of polygons with holes that did not affect the minkowsi sum. Now it simply adds those that do affect it.
 
 [Release 5.5](https://github.com/CGAL/cgal/releases/tag/v5.5)
 -----------

Installation/cmake/modules/FindGMP.cmake

@@ -26,22 +26,25 @@ if( NOT GMP_in_cache )
             NAMES gmp.h
             HINTS ENV GMP_INC_DIR
                   ENV GMP_DIR
+                  $ENV{GMP_DIR}/include
                   ${CGAL_INSTALLATION_PACKAGE_DIR}/auxiliary/gmp/include
             PATH_SUFFIXES include
   	        DOC "The directory containing the GMP header files"
            )
 
-  find_library(GMP_LIBRARY_RELEASE NAMES gmp libgmp-10 mpir
+  find_library(GMP_LIBRARY_RELEASE NAMES gmp libgmp-10 gmp-10 mpir
     HINTS ENV GMP_LIB_DIR
           ENV GMP_DIR
+          $ENV{GMP_DIR}/lib
           ${CGAL_INSTALLATION_PACKAGE_DIR}/auxiliary/gmp/lib
     PATH_SUFFIXES lib
     DOC "Path to the Release GMP library"
     )
 
-  find_library(GMP_LIBRARY_DEBUG NAMES gmpd gmp libgmp-10 mpir
+  find_library(GMP_LIBRARY_DEBUG NAMES gmpd gmp libgmp-10 gmp-10 mpir
     HINTS ENV GMP_LIB_DIR
           ENV GMP_DIR
+          $ENV{GMP_DIR}/include
           ${CGAL_INSTALLATION_PACKAGE_DIR}/auxiliary/gmp/lib
     PATH_SUFFIXES lib
     DOC "Path to the Debug GMP library"

Installation/cmake/modules/FindMPFR.cmake

@@ -25,6 +25,7 @@ if (NOT MPFR_in_cache)
             NAMES mpfr.h
             HINTS ENV MPFR_INC_DIR
                   ENV MPFR_DIR
+                  $ENV{MPFR_DIR}/include
                   ${CGAL_INSTALLATION_PACKAGE_DIR}/auxiliary/gmp/include
             PATH_SUFFIXES include
   	        DOC "The directory containing the MPFR header files"
@@ -33,6 +34,7 @@ if (NOT MPFR_in_cache)
   find_library(MPFR_LIBRARIES NAMES mpfr libmpfr-4 libmpfr-1
     HINTS ENV MPFR_LIB_DIR
           ENV MPFR_DIR
+          $ENV{MPFR_DIR}/lib
           ${CGAL_INSTALLATION_PACKAGE_DIR}/auxiliary/gmp/lib
     PATH_SUFFIXES lib
     DOC "Path to the MPFR library"

Installation/include/CGAL/config.h

@@ -597,4 +597,61 @@ inline std::string data_file_path(const std::string& filename)
 
 } // end namespace CGAL
 
+
+#if BOOST_VERSION < 107900
+
+// Workaround for an accidental enable if of Eigen::Matrix in the
+// boost::multiprecision::cpp_int constructor for some versions of
+// boost
+
+namespace Eigen{
+  template <class A, int B, int C, int D, int E, int F>
+  class Matrix;
+  template <class A, int B, class C>
+  class Ref;
+
+  template <class A, class B, int C>
+  class Product;
+
+  template<typename BinaryOp, typename Lhs, typename Rhs>  class CwiseBinaryOp;
+
+}
+
+namespace boost {
+    namespace multiprecision {
+        namespace detail {
+            template <typename T>
+            struct is_byte_container;
+
+
+            template <class A, int B, int C, int D, int E, int F>
+            struct is_byte_container< Eigen::Matrix<A, B, C, D, E, F>>
+            {
+                static const bool value = false;
+            };
+
+            template <class A, int B, class C>
+            struct is_byte_container< Eigen::Ref<A, B, C>>
+            {
+                static const bool value = false;
+            };
+
+            template <class A, class B, int C>
+            struct is_byte_container< Eigen::Product<A, B, C>>
+            {
+                static const bool value = false;
+            };
+
+            template <class A, class B, class C>
+            struct is_byte_container< Eigen::CwiseBinaryOp<A, B, C>>
+            {
+                static const bool value = false;
+            };
+
+        }
+    }
+}
+
+#endif // BOOST_VERSION < 107900
+
 #endif // CGAL_CONFIG_H

Installation/test/Installation/test_gmp_mpfr_dll.cpp

@@ -8,10 +8,8 @@ int main() {
 #define GMP_SONAME "libgmp-10"
 #define MPFR_SONAME "libmpfr-4"
 #define GMP_SONAME_BACKUP "gmp"
+#define GMP_SONAME_BACKUP_2 "gmp-10"
 #define MPFR_SONAME_BACKUP "mpfr-6"
-#define GMP_MAJOR 5
-#define MPFR_MAJOR 3
-
 
 #include <iostream>
 #include <cassert>
@@ -35,12 +33,21 @@ bool get_version_info(const LPCTSTR name,
     std::cerr << name << " is not loaded!\n";
     return false;
   }
+  else
+    std::cerr << name << " is loaded.\n";
+
   char fileName[_MAX_PATH];
   DWORD size = GetModuleFileName(g_dllHandle, fileName, _MAX_PATH);
   fileName[size] = NULL;
   std::cerr << "Query FileVersion of \"" << fileName << "\"\n";
   DWORD handle = 0;
   size = GetFileVersionInfoSize(fileName, &handle);
+
+  DWORD err = GetLastError();
+  if (size == 0) {
+    std::cerr << "GetFileVersionInfoSize failed with error " << err << std::endl;
+  }
+
   BYTE* versionInfo = new BYTE[size];
   if (!GetFileVersionInfo(fileName, handle, size, versionInfo))
   {
@@ -66,7 +73,9 @@ int main() {
   int major, minor, patch, build;
   if(!get_version_info(GMP_SONAME, major, minor, patch, build)) {
     if(!get_version_info(GMP_SONAME_BACKUP, major, minor, patch, build)) {
-      return 1;
+      if (!get_version_info(GMP_SONAME_BACKUP_2, major, minor, patch, build)) {
+        return 1;
+      }
     }
   }
 

Minkowski_sum_2/include/CGAL/Minkowski_sum_2/Hole_filter_2.h

@@ -26,16 +26,12 @@ namespace CGAL {
  * to the Minkowski sum boundary.
  */
 template <typename Kernel_, typename Container_>
-class Hole_filter_2
+class Hole_filter_2 {
-{
 private:
-  typedef Kernel_               Kernel;
+  using Kernel = Kernel_;
-  typedef Container_            Container;
+  using Container = Container_;
-
+  using Polygon_2 = CGAL::Polygon_2<Kernel, Container>;
-  typedef CGAL::Polygon_2<Kernel, Container>            Polygon_2;
+  using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2<Kernel, Container>;
-  typedef CGAL::Polygon_with_holes_2<Kernel, Container> Polygon_with_holes_2;
-  typedef typename Polygon_with_holes_2::Hole_iterator  Hole_iterator;
-  typedef std::vector<Hole_iterator>                    Hole_iterator_vector;
 
 public:
   /*! Filter out holes of a polygon with holes.
@@ -45,29 +41,17 @@ public:
    */
   void operator()(const Polygon_with_holes_2& pgn1,
                   const Polygon_2& pgn2,
-                  Polygon_with_holes_2& filtered_pgn1) const
+                  Polygon_with_holes_2& filtered_pgn1) const {
-  {
-    filtered_pgn1 = pgn1;
-
-    Hole_iterator_vector to_erase;
     Bbox_2 boundary_bbox = pgn2.bbox();
-
+    filtered_pgn1 = Polygon_with_holes_2(pgn1.outer_boundary());
-    Hole_iterator it = filtered_pgn1.holes_begin();
+    for (const auto& h : pgn1.holes()) {
-    while (it != filtered_pgn1.holes_end()) {
+      Bbox_2 hole_bbox = h.bbox();
-      Bbox_2 hole_bbox = (*it).bbox();
+      if ((hole_bbox.ymax()-hole_bbox.ymin() >=
-
+           boundary_bbox.ymax()-boundary_bbox.ymin()) &&
-      if ((hole_bbox.ymax()-hole_bbox.ymin() <
+          (hole_bbox.xmax()-hole_bbox.xmin() >=
-           boundary_bbox.ymax()-boundary_bbox.ymin()) ||
-          (hole_bbox.xmax()-hole_bbox.xmin() <
            boundary_bbox.xmax()-boundary_bbox.xmin()))
-      {
+        filtered_pgn1.add_hole(h);
-        to_erase.push_back(it);
-      }
-      ++it;
     }
-
-    typename Hole_iterator_vector::iterator it2 = to_erase.begin();
-    while (it2 != to_erase.end()) filtered_pgn1.erase_hole(*it2++);
   }
 
   /*! Filter out holes of a polygon with holes.

Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt

@@ -119,6 +119,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 - `CGAL::Polygon_mesh_processing::fair()`
 - `CGAL::Polygon_mesh_processing::triangulate_face()`
 - `CGAL::Polygon_mesh_processing::triangulate_faces()`
+- `CGAL::Polygon_mesh_processing::triangulate_polygons()`
 - \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::isotropic_remeshing()` \endlink
 - \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::surface_Delaunay_remeshing()` \endlink
 - \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::split_long_edges()` \endlink
@@ -152,6 +153,8 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 \cgalCRPSection{Intersection Functions}
 - `CGAL::Polygon_mesh_processing::does_self_intersect()`
 - `CGAL::Polygon_mesh_processing::self_intersections()`
+- `CGAL::Polygon_mesh_processing::does_triangle_soup_self_intersect()`
+- `CGAL::Polygon_mesh_processing::triangle_soup_self_intersections()`
 - \link PMP_intersection_grp `CGAL::Polygon_mesh_processing::do_intersect()` \endlink
 - `CGAL::Polygon_mesh_processing::intersecting_meshes()`
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp

@@ -41,7 +41,7 @@ public:
 };
 
 
-struct Visitor
+struct Visitor : public CGAL::Polygon_mesh_processing::Triangulate_faces::Default_visitor<Surface_mesh>
 {
    typedef std::unordered_map<face_descriptor,face_descriptor> Container;
 

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

@@ -22,6 +22,7 @@
 #endif
 #include <CGAL/boost/graph/iterator.h>
 #include <CGAL/boost/graph/Euler_operations.h>
+#include <CGAL/use.h>
 #include <vector>
 
 namespace CGAL {
@@ -107,6 +108,11 @@ triangulate_hole_polygon_mesh(PolygonMesh& pmesh,
             Visitor& visitor,
             const typename Kernel::FT max_squared_distance)
 {
+#ifdef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+  CGAL_USE(use_cdt);
+  CGAL_USE(max_squared_distance);
+#endif
+
   typedef Halfedge_around_face_circulator<PolygonMesh>   Hedge_around_face_circulator;
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
   typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;

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

@@ -15,6 +15,11 @@
 
 #include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
 
+#ifdef CGAL_TRIANGULATE_FACES_DO_NOT_USE_CDT2
+# ifndef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+#   define CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+# endif
+#endif
 
 #include <CGAL/value_type_traits.h>
 #ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3
@@ -1442,12 +1447,20 @@ triangulate_hole_polyline_with_cdt(const PointRange& points,
     vertices[v->info()] = v;
   }
 
-  for (std::size_t i = 0; i < size; ++i) {
+  try
-    const std::size_t ip = (i + 1) % size;
+  {
-    if (vertices[i] != vertices[ip]) {
+    for (std::size_t i = 0; i < size; ++i) {
-      cdt.insert_constraint(vertices[i], vertices[ip]);
+      const std::size_t ip = (i + 1) % size;
+      if (vertices[i] != vertices[ip]) {
+        cdt.insert_constraint(vertices[i], vertices[ip]);
+      }
     }
   }
+  catch(const typename CDT::Intersection_of_constraints_exception&)
+  {
+    visitor.end_planar_phase(false);
+    return false;
+  }
 
   // Mark external faces.
   for (typename CDT::All_faces_iterator fit = cdt.all_faces_begin(),

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h

@@ -1125,7 +1125,7 @@ struct Polygon_soup_fixer<PointRange, PolygonRange, std::array<PID, N> >
 ///
 /// \tparam PointRange a model of the concepts `SequenceContainer` and `Swappable`
 ///                    and whose value type is the point type.
-/// \tparam PolygonRange a model of the concept `SequenceContainer`.
+/// \tparam PolygonRange a model of the concept `SequenceContainer`
 ///                      whose `value_type` is itself a model of the concepts `SequenceContainer`,
 ///                      `Swappable`, and `ReversibleContainer` whose `value_type` is `std::size_t`.
 /// \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h

@@ -43,6 +43,7 @@
 #endif
 
 #include <boost/iterator/function_output_iterator.hpp>
+#include <boost/range/irange.hpp>
 
 #include <exception>
 #include <sstream>
@@ -54,6 +55,127 @@ namespace CGAL {
 namespace Polygon_mesh_processing {
 namespace internal {
 
+template <class TM>
+struct Triangle_mesh_and_triangle_soup_wrapper
+{
+  typedef typename boost::graph_traits<TM>::face_descriptor face_descriptor;
+  typedef typename boost::graph_traits<TM>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<TM>::halfedge_descriptor halfedge_descriptor; // private
+
+  static void get_face_vertices(face_descriptor fd, std::array<vertex_descriptor,3>& vh, const TM& tm)
+  {
+    CGAL_assertion(boost::graph_traits<TM>::null_face() != fd);
+    halfedge_descriptor h = halfedge(fd, tm);
+    vh[0]=source(h, tm);
+    vh[1]=target(h, tm);
+    vh[2]=target(next(h, tm), tm);
+  }
+
+  static bool faces_have_a_shared_edge(face_descriptor f, face_descriptor g, std::array<vertex_descriptor, 4>& vh, const TM& tm)
+  {
+    CGAL_assertion(boost::graph_traits<TM>::null_face() != f);
+    CGAL_assertion(boost::graph_traits<TM>::null_face() != g);
+    halfedge_descriptor h=halfedge(f, tm);
+    for(unsigned int i=0; i<3; ++i)
+    {
+      halfedge_descriptor opp_h = opposite(h, tm);
+      if(face(opp_h, tm) == g)
+      {
+        vh[0]=source(h, tm);
+        vh[1]=target(h, tm);
+        vh[2]=target(next(h, tm), tm);
+        vh[3]=target(next(opp_h, tm), tm);
+        return true;
+      }
+      h = next(h, tm);
+    }
+    return false;
+  }
+
+  static bool is_pure_triangle(const TM& tm)
+  {
+    return is_triangle_mesh(tm);
+  }
+};
+
+template <class PointRange, class TriangleRange>
+struct Triangle_mesh_and_triangle_soup_wrapper<
+    std::pair<const PointRange&,
+              const TriangleRange&>>
+{
+  typedef std::size_t face_descriptor;
+  typedef std::size_t vertex_descriptor;
+
+  typedef std::pair<const PointRange&, const TriangleRange& > Soup;
+
+  static void get_face_vertices(face_descriptor fd, std::array<vertex_descriptor,3>& vh, const Soup& soup)
+  {
+    const auto& face = soup.second[fd];
+    vh[0]=face[0];
+    vh[1]=face[1];
+    vh[2]=face[2];
+  }
+
+  static bool faces_have_a_shared_edge(face_descriptor fd, face_descriptor gd, std::array<vertex_descriptor, 4>& vh, const Soup& soup)
+  {
+    const auto& f = soup.second[fd];
+    const auto& g = soup.second[gd];
+
+    for(unsigned int i=0; i<2; ++i) // no need to check f[2] if neither f[0] nor f[1] are shared
+    {
+      for(unsigned int j=0; j<3; ++j)
+      {
+        if (f[i]==g[j])
+        {
+          vh[0]=f[i];
+          vh[1]=f[i+1];
+          vh[2]=f[(i+2)%3];
+
+          if (vh[1]==g[(j+1)%3])
+          {
+            vh[3]=g[(j+2)%3];
+            return true;
+          }
+          if (vh[1]==g[(j+2)%3])
+          {
+            vh[3]=g[(j+1)%3];
+            return true;
+          }
+
+          if (i==0)
+          {
+            vh[1]=f[i];
+            vh[2]=f[(i+1)%3];
+            vh[0]=f[(i+2)%3];
+            if (vh[0]==g[(j+1)%3])
+            {
+              vh[3]=g[(j+2)%3];
+              return true;
+            }
+            if (vh[0]==g[(j+2)%3])
+            {
+              vh[3]=g[(j+1)%3];
+              return true;
+            }
+          }
+
+          return false;
+        }
+      }
+    }
+
+    return false;
+  }
+
+  static bool is_pure_triangle(const Soup& soup)
+  {
+    for (const typename std::iterator_traits<typename TriangleRange::const_iterator>::value_type& t : soup.second)
+      if (t.size()!=3)
+        return false;
+    return true;
+  }
+};
+
 template<typename Output_iterator>
 struct Throw_at_count_reached_functor {
 
@@ -82,59 +204,49 @@ struct Throw_at_count_reached_functor {
 
 // Checks for 'real' intersections, i.e. not simply a shared vertex or edge
 template <class GT, class TM, class VPM>
-bool do_faces_intersect(typename boost::graph_traits<TM>::halfedge_descriptor h,
+bool do_faces_intersect(typename Triangle_mesh_and_triangle_soup_wrapper<TM>::face_descriptor fh,
-                        typename boost::graph_traits<TM>::halfedge_descriptor g,
+                        typename Triangle_mesh_and_triangle_soup_wrapper<TM>::face_descriptor fg,
                         const TM& tmesh,
                         const VPM vpmap,
                         const typename GT::Construct_segment_3& construct_segment,
                         const typename GT::Construct_triangle_3& construct_triangle,
                         const typename GT::Do_intersect_3& do_intersect)
 {
-  typedef typename boost::graph_traits<TM>::vertex_descriptor                 vertex_descriptor;
+  typedef Triangle_mesh_and_triangle_soup_wrapper<TM>       Wrapper;
-  typedef typename boost::graph_traits<TM>::halfedge_descriptor               halfedge_descriptor;
+  typedef typename Wrapper::vertex_descriptor     vertex_descriptor;
 
-  typedef typename GT::Segment_3                                              Segment;
+  typedef typename GT::Segment_3                            Segment;
-  typedef typename GT::Triangle_3                                             Triangle;
+  typedef typename GT::Triangle_3                          Triangle;
 
-  CGAL_assertion(!is_border(h, tmesh));
-  CGAL_assertion(!is_border(g, tmesh));
 
-  vertex_descriptor hv[3], gv[3];
+  std::array<vertex_descriptor, 3> hv, gv;
-  hv[0] = target(h, tmesh);
+  Wrapper::get_face_vertices(fh, hv, tmesh);
-  hv[1] = target(next(h, tmesh), tmesh);
+  Wrapper::get_face_vertices(fg, gv, tmesh);
-  hv[2] = source(h, tmesh);
-
-  gv[0] = target(g, tmesh);
-  gv[1] = target(next(g, tmesh), tmesh);
-  gv[2] = source(g, tmesh);
 
   // check for shared edge
-  for(unsigned int i=0; i<3; ++i)
+  std::array<vertex_descriptor, 4> verts;
+  if (Wrapper::faces_have_a_shared_edge(fh, fg, verts, tmesh))
   {
-    halfedge_descriptor opp_h = opposite(h, tmesh);
+    if (verts[2]==verts[3]) return false; // only for a soup of triangles
-    if(face(opp_h, tmesh) == face(g, tmesh))
-    {
-      // there is an intersection if the four points are coplanar and the triangles overlap
-      if(CGAL::coplanar(get(vpmap, hv[i]),
-                        get(vpmap, hv[(i+1)%3]),
-                        get(vpmap, hv[(i+2)%3]),
-                        get(vpmap, target(next(opp_h, tmesh), tmesh))) &&
-         CGAL::coplanar_orientation(get(vpmap, hv[(i+2)%3]),
-                                    get(vpmap, hv[i]),
-                                    get(vpmap, hv[(i+1)%3]),
-                                    get(vpmap, target(next(opp_h, tmesh), tmesh)))
-           == CGAL::POSITIVE)
-      {
-        return true;
-      }
-      else
-      {
-        // there is a shared edge but no intersection
-        return false;
-      }
-    }
 
-    h = next(h, tmesh);
+    // there is an intersection if the four points are coplanar and the triangles overlap
+    if(CGAL::coplanar(get(vpmap, verts[0]),
+                      get(vpmap, verts[1]),
+                      get(vpmap, verts[2]),
+                      get(vpmap, verts[3])) &&
+       CGAL::coplanar_orientation(get(vpmap, verts[0]),
+                                  get(vpmap, verts[1]),
+                                  get(vpmap, verts[2]),
+                                  get(vpmap, verts[3]))
+         == CGAL::POSITIVE)
+    {
+      return true;
+    }
+    else
+    {
+      // there is a shared edge but no intersection
+      return false;
+    }
   }
 
   // check for shared vertex --> maybe intersection, maybe not
@@ -188,8 +300,6 @@ template <class Box, class TM, class VPM, class GT,
           class OutputIterator>
 struct Strict_intersect_faces // "strict" as in "not sharing a subface"
 {
-  typedef typename boost::graph_traits<TM>::halfedge_descriptor halfedge_descriptor;
-
   mutable OutputIterator m_iterator;
   const TM& m_tmesh;
   const VPM m_vpmap;
@@ -209,10 +319,7 @@ struct Strict_intersect_faces // "strict" as in "not sharing a subface"
 
   void operator()(const Box* b, const Box* c) const
   {
-    const halfedge_descriptor h = halfedge(b->info(), m_tmesh);
+    if(do_faces_intersect<GT>(b->info(), c->info(), m_tmesh, m_vpmap, m_construct_segment, m_construct_triangle, m_do_intersect))
-    const halfedge_descriptor g = halfedge(c->info(), m_tmesh);
-
-    if(do_faces_intersect<GT>(h, g, m_tmesh, m_vpmap, m_construct_segment, m_construct_triangle, m_do_intersect))
       *m_iterator++ = std::make_pair(b->info(), c->info());
   }
 };
@@ -229,15 +336,17 @@ self_intersections_impl(const FaceRange& face_range,
                         const bool throw_on_SI,
                         const NamedParameters& np)
 {
-  CGAL_precondition(CGAL::is_triangle_mesh(tmesh));
+  typedef TriangleMesh                                                                   TM;
+  typedef Triangle_mesh_and_triangle_soup_wrapper<TM>                                    Wrapper;
+
+  CGAL_precondition(Wrapper::is_pure_triangle(tmesh));
 
   using CGAL::parameters::choose_parameter;
   using CGAL::parameters::get_parameter;
   using CGAL::parameters::is_default_parameter;
 
-  typedef TriangleMesh                                                                   TM;
+  typedef typename Wrapper::face_descriptor                                              face_descriptor;
-  typedef typename boost::graph_traits<TM>::halfedge_descriptor                          halfedge_descriptor;
+  typedef typename Wrapper::vertex_descriptor                                            vertex_descriptor;
-  typedef typename boost::graph_traits<TM>::face_descriptor                              face_descriptor;
 
   typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS                                  Box_policy;
   typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, face_descriptor, Box_policy> Box;
@@ -245,9 +354,9 @@ self_intersections_impl(const FaceRange& face_range,
   typedef typename GetGeomTraits<TM, NamedParameters>::type                              GT;
   GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
 
-  typedef typename GetVertexPointMap<TM, NamedParameters>::const_type                    VPM;
+  typedef GetVertexPointMap<TM, NamedParameters>                                  VPM_helper;
-  VPM vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
+  typedef typename VPM_helper::const_type                                                VPM;
-                               get_const_property_map(boost::vertex_point, tmesh));
+  VPM vpmap = VPM_helper::get_const_map(np, tmesh);
 
   const bool do_limit = !(is_default_parameter<NamedParameters, internal_np::maximum_number_t>::value);
   const unsigned int maximum_number = choose_parameter(get_parameter(np, internal_np::maximum_number), 0);
@@ -268,11 +377,13 @@ self_intersections_impl(const FaceRange& face_range,
   // This loop is very cheap, so there is hardly anything to gain from parallelizing it
   for(face_descriptor f : face_range)
   {
-    halfedge_descriptor h = halfedge(f, tmesh);
+    std::array<vertex_descriptor, 3> vh;
+    Wrapper::get_face_vertices(f, vh, tmesh);
+
     typename boost::property_traits<VPM>::reference
-      p = get(vpmap, target(h,tmesh)),
+      p = get(vpmap, vh[0]),
-      q = get(vpmap, target(next(h, tmesh), tmesh)),
+      q = get(vpmap, vh[1]),
-      r = get(vpmap, target(prev(h, tmesh), tmesh));
+      r = get(vpmap, vh[2]);
 
     // tiny fixme: if f is degenerate, we might still have a real intersection between f
     // and another face f', but right now we are not creating a box for f and thus not returning those
@@ -342,7 +453,7 @@ self_intersections_impl(const FaceRange& face_range,
         std::atomic<unsigned int> atomic_counter(counter);
         Throw_functor throwing_count_functor(atomic_counter, maximum_number, std::back_inserter(face_pairs));
         Throwing_after_count_output_iterator count_filter(throwing_count_functor);
-        Filtered_intersecting_faces_filter limited_callback(tmesh, vpmap, gt, count_filter);
+         Filtered_intersecting_faces_filter limited_callback(tmesh, vpmap, gt, count_filter);
         CGAL::box_self_intersection_d<ConcurrencyTag>(box_ptr.begin(), box_ptr.end(), limited_callback, cutoff);
       }
       catch(const CGAL::internal::Throw_at_output_exception&)
@@ -408,9 +519,9 @@ self_intersections_impl(const FaceRange& face_range,
  *
  * collects intersections between a subset of faces of a triangulated surface mesh.
  * Two faces are said to intersect if the corresponding triangles intersect
- * and the intersection is not an edge nor a vertex incident to both faces.
+ * and the intersection is neither an edge nor a vertex incident to both faces.
  *
- * This function depends on the package \ref PkgBoxIntersectionD
+ * This function depends on the package \ref PkgBoxIntersectionD.
  *
  * @pre `CGAL::is_triangle_mesh(tmesh)`
  *
@@ -475,9 +586,9 @@ self_intersections(const FaceRange& face_range,
  *
  * collects intersections between all the faces of a triangulated surface mesh.
  * Two faces are said to intersect if the corresponding triangles intersect
- * and the intersection is not an edge nor a vertex incident to both faces.
+ * and the intersection is neither an edge nor a vertex incident to both faces.
  *
- * This function depends on the package \ref PkgBoxIntersectionD
+ * This function depends on the package \ref PkgBoxIntersectionD.
  *
  * @pre `CGAL::is_triangle_mesh(tmesh)`
  *
@@ -490,9 +601,9 @@ self_intersections(const FaceRange& face_range,
  *
  * @param tmesh the triangulated surface mesh to be checked
  * @param out output iterator to be filled with all pairs of non-adjacent faces that intersect.
-              In case `tmesh` contains some degenerate faces, for each degenerate face `f` a pair `(f,f)`
+ *            In case `tmesh` contains some degenerate faces, for each degenerate face `f` a pair `(f,f)`
-              will be put in `out` before any other self intersection between non-degenerate faces.
+ *            will be put in `out` before any other self intersection between non-degenerate faces. <br>
-              These are the only pairs where degenerate faces will be reported.
+ *            Note that these are the only pairs where degenerate faces will be reported.
  * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
  *
  * \cgalNamedParamsBegin
@@ -541,7 +652,7 @@ self_intersections(const TriangleMesh& tmesh,
  *
  * \brief tests if a set of faces of a triangulated surface mesh self-intersects.
  *
- * This function depends on the package \ref PkgBoxIntersectionD
+ * This function depends on the package \ref PkgBoxIntersectionD.
  *
  * @pre `CGAL::is_triangle_mesh(tmesh)`
  *
@@ -584,8 +695,6 @@ bool does_self_intersect(const FaceRange& face_range,
                          const TriangleMesh& tmesh,
                          const NamedParameters& np = parameters::default_values())
 {
-  CGAL_precondition(CGAL::is_triangle_mesh(tmesh));
-
   try
   {
     CGAL::Emptyset_iterator unused_out;
@@ -614,7 +723,7 @@ bool does_self_intersect(const FaceRange& face_range,
  *
  * \brief tests if a triangulated surface mesh self-intersects.
  *
- * This function depends on the package \ref PkgBoxIntersectionD
+ * This function depends on the package \ref PkgBoxIntersectionD.
  *
  * @pre `CGAL::is_triangle_mesh(tmesh)`
  *
@@ -656,6 +765,219 @@ bool does_self_intersect(const TriangleMesh& tmesh,
   return does_self_intersect<ConcurrencyTag>(faces(tmesh), tmesh, np);
 }
 
+
+#ifndef DOXYGEN_RUNNING
+template <class PointRange, class VPM>
+struct Property_map_for_soup
+{
+  typedef std::size_t key_type;
+  typedef typename boost::property_traits<VPM>::value_type value_type;
+  //typedef typename boost::property_traits<VPM>::category category;
+  typedef boost::readable_property_map_tag category;
+  typedef typename boost::property_traits<VPM>::reference reference;
+
+  const PointRange& points;
+  VPM vpm;
+
+  Property_map_for_soup(const PointRange& points, VPM vpm)
+    : points(points)
+    , vpm(vpm)
+  {}
+
+  inline friend
+  reference get(const Property_map_for_soup<PointRange, VPM>& map, key_type k)
+  {
+    return get(map.vpm, map.points[k]);
+  }
+};
+#endif
+
+/**
+ * \ingroup PMP_intersection_grp
+ *
+ * collects intersections between all the triangles in a triangle soup.
+ *
+ * Two triangles of the soup are said to intersect if the corresponding geometric triangles intersect
+ * and the intersection is neither an edge nor a vertex of both triangles
+ * (with the same point ids, ignoring the orientation for an edge).
+ *
+ * This function depends on the package \ref PkgBoxIntersectionD.
+ *
+ * @tparam ConcurrencyTag enables sequential versus parallel algorithm.
+ *                        Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`.
+ * @tparam PointRange a model of the concept `RandomAccessContainer`
+ *         whose value type is the point type
+ * @tparam TriangleRange a model of the concept `RandomAccessContainer` whose
+ *         value type is a model of the concept `RandomAccessContainer` whose value type is `std::size_t`
+ * @tparam TriangleIdPairOutputIterator a model of `OutputIterator` holding objects of type
+ *   `std::pair<std::size_t,std::size_t>`
+ * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+ *
+ * @param points points of the soup of triangles
+ * @param triangles each element in the range describes a triangle using the indices of the points in `points`
+ * @param out output iterator to be filled with all pairs of ids of triangles intersecting (the id of a triangle is its position in `triangles`).
+ *            In case the triangle soup contains some degenerate triangles, for each degenerate triangle `t` with id `i` a pair `(i,i)`
+ *            will be put in `out` before any other self intersection between non-degenerate triangles.<br>
+ *            Note that these are the only pairs where degenerate triangles will be reported.
+ * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{maximum_number}
+ *     \cgalParamDescription{the maximum number of self intersections that will be detected and returned by the function.}
+ *     \cgalParamType{unsigned int}
+ *     \cgalParamDefault{No limit.}
+ *     \cgalParamExtra{In parallel mode, the number of returned self-intersections is at least `maximum_number`
+ *     (and not exactly that number) as no strong synchronization is put on threads for performance reasons.}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{point_map}
+ *     \cgalParamDescription{a property map associating points to the elements of the range `points`}
+ *     \cgalParamType{a model of `ReadablePropertyMap` whose value type is a point type from a \cgal `Kernel`.}
+ *     \cgalParamDefault{`CGAL::Identity_property_map`}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{geom_traits}
+ *     \cgalParamDescription{an instance of a geometric traits class}
+ *     \cgalParamType{a class model of `PMPSelfIntersectionTraits`}
+ *     \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
+ *     \cgalParamExtra{The geometric traits class must be compatible with the point type of the point map.}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ *
+ * @return `out`
+ *
+ * @sa `does_triangle_soup_self_intersect()`
+ * @sa `self_intersections()`
+ * @sa `does_self_intersect()`
+ */
+template <class ConcurrencyTag = Sequential_tag,
+          class PointRange,
+          class TriangleRange,
+          class TriangleIdPairOutputIterator,
+          class CGAL_NP_TEMPLATE_PARAMETERS>
+TriangleIdPairOutputIterator
+triangle_soup_self_intersections(const PointRange& points,
+                                 const TriangleRange& triangles,
+                                 TriangleIdPairOutputIterator out,
+                                 const CGAL_NP_CLASS& np = parameters::default_values())
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+  using parameters::is_default_parameter;
+
+  typedef typename CGAL::GetPointMap<PointRange, CGAL_NP_CLASS>::const_type Point_map_base;
+  Point_map_base pm_base = choose_parameter<Point_map_base>(get_parameter(np, internal_np::point_map));
+  typedef Property_map_for_soup<PointRange, Point_map_base> Point_map;
+  typedef typename GetPolygonSoupGeomTraits<PointRange, CGAL_NP_CLASS>::type GT;
+  GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
+
+  const bool do_limit = !(is_default_parameter<CGAL_NP_CLASS, internal_np::maximum_number_t>::value);
+  if (do_limit)
+  {
+    return self_intersections<ConcurrencyTag>(boost::irange<std::size_t>(0, triangles.size()),
+                                              std::make_pair(std::cref(points), std::cref(triangles)),
+                                              out,
+                                              parameters::vertex_point_map(Point_map(points,pm_base)).
+                                              geom_traits(gt).
+                                              maximum_number(choose_parameter(get_parameter(np, internal_np::maximum_number), 0)));
+  }
+
+  return self_intersections<ConcurrencyTag>(boost::irange<std::size_t>(0, triangles.size()),
+                                            std::make_pair(std::cref(points), std::cref(triangles)),
+                                            out,
+                                            parameters::vertex_point_map(Point_map(points,pm_base)).
+                                            geom_traits(gt));
+}
+
+/**
+ * \ingroup PMP_intersection_grp
+ *
+ * \brief tests if a triangle soup self-intersects.
+ *
+ * A triangle soup self-intersects if at least two triangles of the soup intersect.
+ * Two triangles of the soup are said to intersect if the corresponding geometric triangles intersect
+ * and the intersection is neither an edge nor a vertex of both triangles
+ * (with the same point ids, ignoring the orientation for an edge).
+ *
+ * This function depends on the package \ref PkgBoxIntersectionD.
+ *
+ * @tparam ConcurrencyTag enables sequential versus parallel algorithm.
+ *                        Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`.
+ * @tparam PointRange a model of the concept `RandomAccessContainer`
+ *         whose value type is the point type
+ * @tparam TriangleRange a model of the concept `RandomAccessContainer` whose
+ *         value type is a model of the concept `RandomAccessContainer` whose value type is `std::size_t`
+ * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+ *
+ * @param points points of the soup of triangles
+ * @param triangles each element in the range describes a triangle using the indices of the points in `points`
+ * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{point_map}
+ *     \cgalParamDescription{a property map associating points to the elements of the range `points`}
+ *     \cgalParamType{a model of `ReadablePropertyMap` whose value type is a point type from a \cgal `Kernel`.}
+ *     \cgalParamDefault{`CGAL::Identity_property_map`}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{geom_traits}
+ *     \cgalParamDescription{an instance of a geometric traits class}
+ *     \cgalParamType{a class model of `PMPSelfIntersectionTraits`}
+ *     \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
+ *     \cgalParamExtra{The geometric traits class must be compatible with the point type of the point map.}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ *
+ * @return `true` if the triangle soup self-intersects, and `false` otherwise.
+ *
+ * @sa `triangle_soup_self_intersections()`
+ * @sa `self_intersections()`
+ * @sa `does_self_intersect()`
+ */
+template <class ConcurrencyTag = Sequential_tag,
+          class PointRange,
+          class TriangleRange,
+          class CGAL_NP_TEMPLATE_PARAMETERS>
+bool does_triangle_soup_self_intersect(const PointRange& points,
+                                       const TriangleRange& triangles,
+                                       const CGAL_NP_CLASS& np = parameters::default_values())
+{
+  try
+  {
+    using parameters::choose_parameter;
+    using parameters::get_parameter;
+
+    CGAL::Emptyset_iterator unused_out;
+    typedef typename CGAL::GetPointMap<PointRange, CGAL_NP_CLASS>::const_type Point_map_base;
+    Point_map_base pm_base = choose_parameter<Point_map_base>(get_parameter(np, internal_np::point_map));
+    typedef Property_map_for_soup<PointRange, Point_map_base> Point_map;
+    typedef typename GetPolygonSoupGeomTraits<PointRange, CGAL_NP_CLASS>::type GT;
+    GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
+
+    internal::self_intersections_impl<ConcurrencyTag>(boost::irange<std::size_t>(0, triangles.size()),
+                                                      std::make_pair(std::cref(points), std::cref(triangles)),
+                                                      unused_out, true /*throw*/,
+                                                      parameters::vertex_point_map(Point_map(points,pm_base))
+                                                                 .geom_traits(gt));
+  }
+  catch (const CGAL::internal::Throw_at_output_exception&)
+  {
+    return true;
+  }
+  #if defined(CGAL_LINKED_WITH_TBB) && TBB_USE_CAPTURED_EXCEPTION
+  catch (const tbb::captured_exception& e)
+  {
+    const char* ti1 = e.name();
+    const char* ti2 = typeid(const CGAL::internal::Throw_at_output_exception&).name();
+    const std::string tn1(ti1);
+    const std::string tn2(ti2);
+    if (tn1 == tn2) return true;
+    else throw;
+  }
+  #endif
+  return false;
+}
+
 }// namespace Polygon_mesh_processing
 }// namespace CGAL
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_faces.h

@@ -15,38 +15,27 @@
 
 #include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
 
-#include <CGAL/disable_warnings.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
 
+#include <CGAL/array.h>
 #include <CGAL/boost/graph/helpers.h>
 #include <CGAL/boost/graph/Euler_operations.h>
-
-#ifndef CGAL_TRIANGULATE_FACES_DO_NOT_USE_CDT2
-#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/Projection_traits_3.h>
-#else
-#include <CGAL/use.h>
-#endif
-
-#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
-#include <CGAL/Polygon_mesh_processing/compute_normal.h>
-#include <CGAL/Named_function_parameters.h>
 #include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Named_function_parameters.h>
 
-#include <boost/range/size.hpp>
+#include <boost/range/value_type.hpp>
 
+#include <algorithm>
+#include <iterator>
+#include <map>
 #include <queue>
-#include <vector>
 #include <utility>
-#include <CGAL/array.h>
+#include <vector>
 
 namespace CGAL {
-
 namespace Polygon_mesh_processing {
+namespace Triangulate_faces {
 
-namespace Triangulate_faces
-{
 /** \ingroup PMP_meshing_grp
 *   %Default new face visitor model of `PMPTriangulateFaceVisitor`.
 *   All its functions have an empty body. This class can be used as a
@@ -54,7 +43,9 @@ namespace Triangulate_faces
 *   overridden.
 */
 template<class PolygonMesh>
-struct Default_visitor {
+struct Default_visitor
+  : public Hole_filling::Default_visitor
+{
   typedef boost::graph_traits<PolygonMesh> GT;
   typedef typename GT::face_descriptor face_descriptor;
 
@@ -63,245 +54,31 @@ struct Default_visitor {
   void after_subface_created(face_descriptor /*f_new*/) {}
 };
 
-} //end namespace Triangulate_faces
+} // namespace Triangulate_faces
 
 namespace internal {
 
-template <class PM
+template <typename PolygonMesh>
-          , typename VertexPointMap
+class Triangulate_polygon_mesh_modifier
-          , typename Kernel
-          , typename Visitor>
-class Triangulate_modifier
 {
-  typedef Kernel Traits;
+  using vertex_descriptor = typename boost::graph_traits<PolygonMesh>::vertex_descriptor;
+  using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
+  using face_descriptor = typename boost::graph_traits<PolygonMesh>::face_descriptor;
 
-  typedef typename boost::graph_traits<PM>::vertex_descriptor vertex_descriptor;
+private:
-  typedef typename boost::graph_traits<PM>::halfedge_descriptor halfedge_descriptor;
+  template <typename VPM,
-  typedef typename boost::graph_traits<PM>::face_descriptor face_descriptor;
+            typename Visitor,
-  typedef typename boost::graph_traits<PM>::edge_descriptor edge_descriptor;
+            typename NamedParameters>
-  typedef typename Kernel::Point_3 Point;
+  bool triangulate_face_with_hole_filling(face_descriptor f,
-
+                                          PolygonMesh& pmesh,
-  struct Face_info {
+                                          const VPM vpm,
-    typename boost::graph_traits<PM>::halfedge_descriptor e[3];
+                                          Visitor visitor,
-    bool is_external;
+                                          const NamedParameters& np)
-  };
-
-  typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
-  VertexPointMap _vpmap;
-  Traits _traits;
-
-public:
-  Triangulate_modifier(VertexPointMap vpmap, const Traits& traits = Traits())
-    : _vpmap(vpmap), _traits(traits)
-  {
-  }
-
-  template <class Face_handle>
-  bool is_external(Face_handle fh) const {
-    return fh->info().is_external;
-  }
-
-  bool triangulate_face(face_descriptor f, PM& pmesh, bool use_cdt, Visitor visitor)
-  {
-    typedef typename Traits::FT FT;
-
-    typename Traits::Vector_3 normal =
-      Polygon_mesh_processing::compute_face_normal(
-        f, pmesh, CGAL::parameters::geom_traits(_traits)
-                                   .vertex_point_map(_vpmap));
-
-    if(normal == typename Traits::Vector_3(0,0,0))
-      return false;
-
-    std::size_t original_size = CGAL::halfedges_around_face(halfedge(f, pmesh), pmesh).size();
-    if(original_size == 4)
-    {
-      halfedge_descriptor v0, v1, v2, v3;
-      v0 = halfedge(f, pmesh);
-      Point_ref p0 = get(_vpmap, target(v0, pmesh));
-      v1 = next(v0, pmesh);
-      Point_ref p1 = get(_vpmap, target(v1, pmesh));
-      v2 = next(v1, pmesh);
-      Point_ref p2 = get(_vpmap, target(v2, pmesh));
-      v3 = next(v2, pmesh);
-      Point_ref p3 = get(_vpmap, target(v3, pmesh));
-
-      /* Chooses the diagonal that will split the quad in two triangles that maximize
-       * the scalar product of of the un-normalized normals of the two triangles.
-       * The lengths of the un-normalized normals (computed using cross-products of two vectors)
-       *  are proportional to the area of the triangles.
-       * Maximize the scalar product of the two normals will avoid skinny triangles,
-       * and will also taken into account the cosine of the angle between the two normals.
-       * In particular, if the two triangles are oriented in different directions,
-       * the scalar product will be negative.
-       */
-      FT p1p3 = CGAL::cross_product(p2-p1,p3-p2) * CGAL::cross_product(p0-p3,p1-p0);
-      FT p0p2 = CGAL::cross_product(p1-p0,p1-p2) * CGAL::cross_product(p3-p2,p3-p0);
-      visitor.before_subface_creations(f);
-      halfedge_descriptor res = (p0p2>p1p3)
-                              ?  CGAL::Euler::split_face(v0, v2, pmesh)
-                              :  CGAL::Euler::split_face(v1, v3, pmesh);
-
-      visitor.after_subface_created(face(res,pmesh));
-      visitor.after_subface_created(face(opposite(res,pmesh),pmesh));
-
-      visitor.after_subface_creations();
-    }
-    else
-    {
-#ifndef CGAL_TRIANGULATE_FACES_DO_NOT_USE_CDT2
-      if (use_cdt)
-      {
-        typedef CGAL::Projection_traits_3<Traits>   P_traits;
-        typedef CGAL::Triangulation_vertex_base_with_info_2<halfedge_descriptor,
-                                                            P_traits>        Vb;
-        typedef CGAL::Triangulation_face_base_with_info_2<Face_info,
-                                                          P_traits>          Fb1;
-        typedef CGAL::Constrained_triangulation_face_base_2<P_traits, Fb1>   Fb;
-        typedef CGAL::Triangulation_data_structure_2<Vb,Fb>                  TDS;
-        typedef CGAL::Exact_intersections_tag                                Itag;
-        typedef CGAL::Constrained_Delaunay_triangulation_2<P_traits,
-                                                           TDS,
-                                                           Itag>             CDT;
-        P_traits cdt_traits(normal);
-        CDT cdt(cdt_traits);
-        return triangulate_face_with_CDT(f, pmesh, cdt, visitor);
-      }
-#else
-      CGAL_USE(use_cdt);
-#endif
-      return triangulate_face_with_hole_filling(f, pmesh, visitor);
-    }
-    return true;
-  }
-
-  template<class CDT>
-  bool triangulate_face_with_CDT(face_descriptor f, PM& pmesh, CDT& cdt, Visitor visitor)
-  {
-    std::size_t original_size = CGAL::halfedges_around_face(halfedge(f, pmesh), pmesh).size();
-
-    // Halfedge_around_facet_circulator
-    typedef typename CDT::Vertex_handle Tr_Vertex_handle;
-    halfedge_descriptor start = halfedge(f, pmesh);
-    halfedge_descriptor h = start;
-    Tr_Vertex_handle previous, first;
-    do
-    {
-      Tr_Vertex_handle vh = cdt.insert(get(_vpmap, target(h, pmesh)));
-      if (first == Tr_Vertex_handle()) {
-        first = vh;
-      }
-      vh->info() = h;
-      if(previous != Tr_Vertex_handle() && previous != vh) {
-        cdt.insert_constraint(previous, vh);
-      }
-      previous = vh;
-      h = next(h, pmesh);
-
-    } while( h != start );
-    cdt.insert_constraint(previous, first);
-
-    // sets mark is_external
-    for(typename CDT::All_faces_iterator fit = cdt.all_faces_begin(),
-          end = cdt.all_faces_end();
-        fit != end; ++fit)
-    {
-      fit->info().is_external = 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().is_external)
-        continue;
-
-      fh->info().is_external = 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() != original_size)
-      return false;
-
-
-    // then modify the polyhedron
-    visitor.before_subface_creations(f);
-    // make_hole. (see comment in function body)
-    this->make_hole(halfedge(f, pmesh), pmesh);
-
-    for(typename CDT::Finite_edges_iterator eit = cdt.finite_edges_begin(),
-          end = cdt.finite_edges_end();
-        eit != end; ++eit)
-    {
-      typename CDT::Face_handle fh = eit->first;
-      const int index = eit->second;
-      typename CDT::Face_handle opposite_fh = fh->neighbor(eit->second);
-      const int opposite_index = opposite_fh->index(fh);
-
-      const Tr_Vertex_handle va = fh->vertex(cdt. cw(index));
-      const Tr_Vertex_handle vb = fh->vertex(cdt.ccw(index));
-
-      if( ! (is_external(fh) && is_external(opposite_fh))//not both fh are external
-          && ! cdt.is_constrained(*eit) )                  //and edge is not constrained
-      {
-        // strictly internal edge
-        halfedge_descriptor hnew = halfedge(add_edge(pmesh), pmesh),
-          hnewopp = opposite(hnew, pmesh);
-
-        fh->info().e[index] = hnew;
-        opposite_fh->info().e[opposite_index] = hnewopp;
-
-        set_target(hnew,    target(va->info(), pmesh), pmesh);
-        set_target(hnewopp, target(vb->info(), pmesh), pmesh);
-      }
-      if( cdt.is_constrained(*eit) ) //edge is constrained
-      {
-        if(!is_external(fh)) {
-          fh->info().e[index] = va->info();
-        }
-        if(!is_external(opposite_fh)) {
-          opposite_fh->info().e[opposite_index] = vb->info();
-        }
-      }
-    }
-    for(typename CDT::Finite_faces_iterator fit = cdt.finite_faces_begin(),
-          end = cdt.finite_faces_end();
-        fit != end; ++fit)
-    {
-      if(!is_external(fit))
-      {
-        halfedge_descriptor h0 = fit->info().e[0];
-        halfedge_descriptor h1 = fit->info().e[1];
-        halfedge_descriptor h2 = fit->info().e[2];
-        CGAL_assertion(h0 != halfedge_descriptor());
-        CGAL_assertion(h1 != halfedge_descriptor());
-        CGAL_assertion(h2 != halfedge_descriptor());
-
-        set_next(h0, h1, pmesh);
-        set_next(h1, h2, pmesh);
-        set_next(h2, h0, pmesh);
-
-        Euler::fill_hole(h0, pmesh);
-        visitor.after_subface_created(face(h0, pmesh));
-      }
-    }
-    visitor.after_subface_creations();
-    return true;
-  }
-
-  bool triangulate_face_with_hole_filling(face_descriptor f, PM& pmesh, Visitor visitor)
   {
     namespace PMP = CGAL::Polygon_mesh_processing;
 
+    using Point = typename boost::property_traits<VPM>::value_type;
+
     // gather halfedges around the face
     std::vector<Point> hole_points;
     std::vector<vertex_descriptor> border_vertices;
@@ -309,60 +86,60 @@ public:
     for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, pmesh), pmesh))
     {
       vertex_descriptor v = source(h, pmesh);
-      hole_points.push_back( get(_vpmap, v) );
+      hole_points.push_back(get(vpm, v));
       border_vertices.push_back(v);
     }
 
     // use hole filling
     typedef CGAL::Triple<int, int, int> Face_indices;
     std::vector<Face_indices> patch;
-    PMP::triangulate_hole_polyline(hole_points, std::back_inserter(patch),
+    PMP::triangulate_hole_polyline(hole_points, std::back_inserter(patch), np);
-                                   parameters::geom_traits(_traits));
 
     if(patch.empty())
       return false;
 
     // triangulate the hole
-    std::map< std::pair<int, int> , halfedge_descriptor > halfedge_map;
+    std::map<std::pair<int, int>, halfedge_descriptor > halfedge_map;
-    int i=0;
+    int i = 0;
     for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, pmesh), pmesh))
     {
       int j = std::size_t(i+1) == hole_points.size() ? 0 : i+1;
-      halfedge_map[ std::make_pair(i, j) ] = h;
+      halfedge_map[std::make_pair(i, j)] = h;
       ++i;
     }
 
     visitor.before_subface_creations(f);
+
     bool first = true;
     std::vector<halfedge_descriptor> hedges;
     hedges.reserve(4);
     for(const Face_indices& triangle : patch)
     {
-      if (first)
+      if(first)
-        first=false;
+        first = false;
       else
-        f=add_face(pmesh);
+        f = add_face(pmesh);
+
       visitor.after_subface_created(f);
 
-      std::array<int, 4> indices =
+      std::array<int, 4> indices = make_array(triangle.first,
-        make_array( triangle.first,
+                                              triangle.second,
-                    triangle.second,
+                                              triangle.third,
-                    triangle.third,
+                                              triangle.first);
-                    triangle.first );
       for (int i=0; i<3; ++i)
       {
         typename std::map< std::pair<int, int> , halfedge_descriptor >::iterator insert_res =
-          halfedge_map.insert(
+          halfedge_map.emplace(std::make_pair(indices[i], indices[i+1]),
-            std::make_pair( std::make_pair(indices[i], indices[i+1]),
+                               boost::graph_traits<PolygonMesh>::null_halfedge()).first;
-                            boost::graph_traits<PM>::null_halfedge() ) ).first;
+        if(insert_res->second == boost::graph_traits<PolygonMesh>::null_halfedge())
-        if (insert_res->second == boost::graph_traits<PM>::null_halfedge())
         {
           halfedge_descriptor nh = halfedge(add_edge(pmesh), pmesh);
-          insert_res->second=nh;
+          insert_res->second = nh;
-          halfedge_map[std::make_pair(indices[i+1], indices[i])]=opposite(nh, pmesh);
+          halfedge_map[std::make_pair(indices[i+1], indices[i])] = opposite(nh, pmesh);
         }
         hedges.push_back(insert_res->second);
       }
+
       hedges.push_back(hedges.front());
       for(int i=0; i<3;++i)
       {
@@ -370,59 +147,93 @@ public:
         set_face(hedges[i], f, pmesh);
         set_target(hedges[i], border_vertices[indices[i+1]], pmesh);
       }
+
       set_halfedge(f, hedges[0], pmesh);
       hedges.clear();
     }
+
     visitor.after_subface_creations();
+
     return true;
   }
 
-  template<typename FaceRange>
+public:
-  bool operator()(FaceRange face_range, PM& pmesh, bool use_cdt, Visitor visitor)
+  template <typename NamedParameters>
+  bool operator()(face_descriptor f,
+                  PolygonMesh& pmesh,
+                  const NamedParameters& np)
   {
-   bool result = true;
+    using Traits = typename GetGeomTraits<PolygonMesh, NamedParameters>::type;
-    // One need to store facet handles into a vector, because the list of
+    using VPM = typename GetVertexPointMap<PolygonMesh, NamedParameters>::type;
-    // facets of the polyhedron will be modified during the loop, and
-    // that invalidates the range [facets_begin(), facets_end()[.
-    std::vector<face_descriptor> facets;
-    facets.reserve(std::distance(boost::begin(face_range), boost::end(face_range)));
 
-    //only consider non-triangular faces
+    using FT = typename Traits::FT;
-    for(face_descriptor fit : face_range)
+    using Point_ref = typename boost::property_traits<VPM>::reference;
-      if ( next( next( halfedge(fit, pmesh), pmesh), pmesh)
-        !=       prev( halfedge(fit, pmesh), pmesh) )
-        facets.push_back(fit);
 
-    // Iterates on the vector of face descriptors
+    using Visitor = typename internal_np::Lookup_named_param_def<
-    for(face_descriptor f : facets)
+                                            internal_np::visitor_t,
+                                            NamedParameters,
+                                            Triangulate_faces::Default_visitor<PolygonMesh> // default
+                                          >::type;
+
+    using parameters::choose_parameter;
+    using parameters::get_parameter;
+
+    CGAL_precondition(is_valid_face_descriptor(f, pmesh));
+
+    Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
+    VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                               get_property_map(vertex_point, pmesh));
+    Visitor visitor = choose_parameter<Visitor>(get_parameter(np, internal_np::visitor),
+                                                Triangulate_faces::Default_visitor<PolygonMesh>());
+
+    typename Traits::Construct_cross_product_vector_3 cross_product =
+      traits.construct_cross_product_vector_3_object();
+
+    typename boost::graph_traits<PolygonMesh>::degree_size_type  original_size = degree(f, pmesh);
+    if(original_size <= 3)
+      return true;
+
+    if(original_size == 4)
     {
-      if(!this->triangulate_face(f, pmesh, use_cdt, visitor))
+      halfedge_descriptor v0, v1, v2, v3;
-       result = false;
+      v0 = halfedge(f, pmesh);
+      Point_ref p0 = get(vpm, target(v0, pmesh));
+      v1 = next(v0, pmesh);
+      Point_ref p1 = get(vpm, target(v1, pmesh));
+      v2 = next(v1, pmesh);
+      Point_ref p2 = get(vpm, target(v2, pmesh));
+      v3 = next(v2, pmesh);
+      Point_ref p3 = get(vpm, target(v3, pmesh));
+
+      /* Chooses the diagonal that will split the quad in two triangles that maximize
+       * the scalar product of the un-normalized normals of the two triangles.
+       * The lengths of the un-normalized normals (computed using cross-products of two vectors)
+       *  are proportional to the area of the triangles.
+       * Maximize the scalar product of the two normals will avoid skinny triangles,
+       * and will also taken into account the cosine of the angle between the two normals.
+       * In particular, if the two triangles are oriented in different directions,
+       * the scalar product will be negative.
+       */
+      visitor.before_subface_creations(f);
+
+      const FT p1p3 = cross_product(p2-p1, p3-p2) * cross_product(p0-p3, p1-p0);
+      const FT p0p2 = cross_product(p1-p0, p1-p2) * cross_product(p3-p2, p3-p0);
+      halfedge_descriptor res = (p0p2>p1p3) ?  CGAL::Euler::split_face(v0, v2, pmesh)
+                                            :  CGAL::Euler::split_face(v1, v3, pmesh);
+
+      visitor.after_subface_created(face(res, pmesh));
+      visitor.after_subface_created(face(opposite(res, pmesh), pmesh));
+
+      visitor.after_subface_creations();
+
+      return true;
     }
-    return result;
+
+    return triangulate_face_with_hole_filling(f, pmesh, vpm, visitor, np);
   }
+}; // class Triangulate_polygon_mesh_modifier
 
-  void make_hole(halfedge_descriptor h, PM& pmesh)
+} // namespace internal
-  {
-    //we are not using Euler::make_hole because it has a precondition
-    //that the hole is not made on the boundary of the mesh
-    //here we allow making a hole on the boundary, and the pair(s) of
-    //halfedges that become border-border are fixed by the connectivity
-    //setting made in operator()
-    CGAL_assertion(!is_border(h, pmesh));
-    face_descriptor fd = face(h, pmesh);
-
-    for(halfedge_descriptor hd : halfedges_around_face(h, pmesh))
-    {
-      CGAL::internal::set_border(hd, pmesh);
-    }
-    remove_face(fd, pmesh);
-  }
-
-
-}; // end class Triangulate_modifier
-
-}//end namespace internal
 
 /**
 * \ingroup PMP_meshing_grp
@@ -455,50 +266,30 @@ public:
 *
 *   \cgalParamNBegin{visitor}
 *     \cgalParamDescription{a visitor that enables to track how faces are triangulated into subfaces}
-*     \cgalParamType{a class model of `PMPTriangulateFaceVisitor`}
+*     \cgalParamType{a class model of `PMPTriangulateFaceVisitor` and `PMPHolefillingVisitor`}
 *     \cgalParamDefault{`Triangulate_faces::Default_visitor<PolygonMesh>`}
 *     \cgalParamExtra{Note that the visitor will be copied, so
 *                     it must not have any data member that does not have a reference-like type.}
 *   \cgalParamNEnd
 * \cgalNamedParamsEnd
 *
+* This function calls `CGAL::Polygon_mesh_processing::triangulate_hole_polyline()`.
+* Refer to its documentation for its named parameters.
+*
+* @pre The face `f` is not degenerate.
+*
 * @return `true` if the face has been triangulated.
 *
 * @see `triangulate_faces()`
 */
-template<typename PolygonMesh, typename NamedParameters = parameters::Default_named_parameters>
+template <typename PolygonMesh,
+          typename NamedParameters = parameters::Default_named_parameters>
 bool triangulate_face(typename boost::graph_traits<PolygonMesh>::face_descriptor f,
                       PolygonMesh& pmesh,
                       const NamedParameters& np = parameters::default_values())
 {
-  using parameters::choose_parameter;
+  internal::Triangulate_polygon_mesh_modifier<PolygonMesh> modifier;
-  using parameters::get_parameter;
+  return modifier(f, pmesh, np);
-
-  CGAL_precondition(is_valid_face_descriptor(f, pmesh));
-
-  //VertexPointMap
-  typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::type VPMap;
-  VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
-                             get_property_map(vertex_point, pmesh));
-
-  //Kernel
-  typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type Kernel;
-  Kernel traits = choose_parameter<Kernel>(get_parameter(np, internal_np::geom_traits));
-
-  //Option
-  bool use_cdt = choose_parameter(get_parameter(np, internal_np::use_delaunay_triangulation), true);
-
-  typedef typename internal_np::Lookup_named_param_def<
-    internal_np::visitor_t,
-    NamedParameters,
-    Triangulate_faces::Default_visitor<PolygonMesh>//default
-  >::type Visitor;
-  Visitor visitor = choose_parameter<Visitor>(
-                             get_parameter(np, internal_np::visitor),
-                             Triangulate_faces::Default_visitor<PolygonMesh>());
-
-  internal::Triangulate_modifier<PolygonMesh, VPMap, Kernel, Visitor> modifier(vpmap, traits);
-  return modifier.triangulate_face(f, pmesh, use_cdt, visitor);
 }
 
 /**
@@ -534,48 +325,47 @@ bool triangulate_face(typename boost::graph_traits<PolygonMesh>::face_descriptor
 *
 *   \cgalParamNBegin{visitor}
 *     \cgalParamDescription{a visitor that enables to track how faces are triangulated into subfaces}
-*     \cgalParamType{a class model of `PMPTriangulateFaceVisitor`}
+*     \cgalParamType{a class model of `PMPTriangulateFaceVisitor` and `PMPHolefillingVisitor`}
 *     \cgalParamDefault{`Triangulate_faces::Default_visitor<PolygonMesh>`}
 *     \cgalParamExtra{Note that the visitor will be copied, so
 *                     it must not have any data member that does not have a reference-like type.}
 *   \cgalParamNEnd
 * \cgalNamedParamsEnd
 *
+* This function calls `CGAL::Polygon_mesh_processing::triangulate_hole_polyline()` for each face.
+* Refer to its documentation for its named parameters.
+*
+* @pre No face within `face_range` is degenerate.
+*
 * @return `true` if all the faces have been triangulated.
 *
 * @see `triangulate_face()`
+* @see `triangulate_polygons()`
 */
-template <typename FaceRange, typename PolygonMesh, typename NamedParameters = parameters::Default_named_parameters>
+template <typename FaceRange,
+          typename PolygonMesh,
+          typename NamedParameters = parameters::Default_named_parameters>
 bool triangulate_faces(FaceRange face_range,
                        PolygonMesh& pmesh,
                        const NamedParameters& np = parameters::default_values())
 {
-  using parameters::choose_parameter;
+  using face_descriptor = typename boost::graph_traits<PolygonMesh>::face_descriptor;
-  using parameters::get_parameter;
 
-  //VertexPointMap
+  bool result = true;
-  typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::type VPMap;
-  VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
-                             get_property_map(vertex_point, pmesh));
 
-  //Kernel
+  // One needs to store the facets into a vector, because the list of
-  typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type Kernel;
+  // facets of the polyhedron will be modified during the loop, and
-  Kernel traits = choose_parameter<Kernel>(get_parameter(np, internal_np::geom_traits));
+  // that invalidates the range [facets_begin(), facets_end()[.
+  std::vector<face_descriptor> facets(std::begin(face_range), std::end(face_range));
 
-  //Option
+  internal::Triangulate_polygon_mesh_modifier<PolygonMesh> modifier;
-  bool use_cdt = choose_parameter(get_parameter(np, internal_np::use_delaunay_triangulation), true);
+  for(face_descriptor f : facets)
+  {
+    if(!modifier(f, pmesh, np))
+      result = false;
+  }
 
-  typedef typename internal_np::Lookup_named_param_def<
+  return result;
-    internal_np::visitor_t,
-    NamedParameters,
-    Triangulate_faces::Default_visitor<PolygonMesh>//default
-  >::type Visitor;
-  Visitor visitor = choose_parameter<Visitor>(
-                                  get_parameter(np, internal_np::visitor),
-                                  Triangulate_faces::Default_visitor<PolygonMesh>());
-
-  internal::Triangulate_modifier<PolygonMesh, VPMap, Kernel, Visitor> modifier(vpmap, traits);
-  return modifier(face_range, pmesh, use_cdt, visitor);
 }
 
 /**
@@ -608,16 +398,22 @@ bool triangulate_faces(FaceRange face_range,
 *
 *   \cgalParamNBegin{visitor}
 *     \cgalParamDescription{a visitor that enables to track how faces are triangulated into subfaces}
-*     \cgalParamType{a class model of `PMPTriangulateFaceVisitor`}
+*     \cgalParamType{a class model of `PMPTriangulateFaceVisitor` and `PMPHolefillingVisitor`}
 *     \cgalParamDefault{`Triangulate_faces::Default_visitor<PolygonMesh>`}
 *     \cgalParamExtra{Note that the visitor will be copied, so
 *                     it must not have any data member that does not have a reference-like type.}
 *   \cgalParamNEnd
 * \cgalNamedParamsEnd
 *
+* This function calls `CGAL::Polygon_mesh_processing::triangulate_hole_polyline()` on all the faces of the polygon mesh.
+* Refer to its documentation for its named parameters.
+*
+* @pre No face of `pmesh` is degenerate.
+*
 * @return `true` if all the faces have been triangulated.
 *
 * @see `triangulate_face()`
+* @see `triangulate_polygons()`
 */
 template <typename PolygonMesh, typename NamedParameters = parameters::Default_named_parameters>
 bool triangulate_faces(PolygonMesh& pmesh,
@@ -626,10 +422,255 @@ bool triangulate_faces(PolygonMesh& pmesh,
   return triangulate_faces(faces(pmesh), pmesh, np);
 }
 
-} // end namespace Polygon_mesh_processing
+////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// Polygon Soup
 
-} // end namespace CGAL
+namespace Triangulate_polygons {
 
-#include <CGAL/enable_warnings.h>
+/** \ingroup PMP_meshing_grp
+*   %Default new polygon visitor model of `PMPTriangulateFaceVisitor`.
+*   All its functions have an empty body. This class can be used as a
+*   base class if only some of the functions of the concept require to be
+*   overridden.
+*/
+struct Default_visitor
+  : public Hole_filling::Default_visitor
+{
+  template <typename Polygon>
+  void before_subface_creations(const Polygon& /*f_old*/) {}
+
+  template <typename Polygon>
+  void after_subface_created(const Polygon& /*f_new*/) {}
+
+  void after_subface_creations() {}
+};
+
+} // namespace Triangulate_polygons
+
+namespace internal {
+
+class Triangulate_polygon_soup_modifier
+{
+private:
+  template<typename Polygon,
+           typename PointRange,
+           typename PolygonRange,
+           typename PMap,
+           typename Visitor,
+           typename NamedParameters>
+  bool triangulate_polygon_with_hole_filling(const Polygon& polygon,
+                                             const PointRange& points,
+                                             PolygonRange& triangulated_polygons, // output
+                                             PMap pm,
+                                             Visitor visitor,
+                                             const NamedParameters& np)
+  {
+    namespace PMP = CGAL::Polygon_mesh_processing;
+
+    using Point = typename boost::property_traits<PMap>::value_type;
+
+    // gather halfedges around the face
+    std::vector<Point> hole_points;
+    std::vector<std::size_t> hole_points_indices;
+
+    for(std::size_t i : polygon)
+    {
+      hole_points.push_back(get(pm, points[i]));
+      hole_points_indices.push_back(i);
+    }
+
+    // use hole filling
+    typedef CGAL::Triple<int, int, int> Face_indices;
+    std::vector<Face_indices> patch;
+    PMP::triangulate_hole_polyline(hole_points, std::back_inserter(patch), np);
+
+    if(patch.empty())
+      return false;
+
+    visitor.before_subface_creations(polygon);
+
+    for(const Face_indices& triangle : patch)
+    {
+      triangulated_polygons.push_back({hole_points_indices[triangle.first],
+                                       hole_points_indices[triangle.second],
+                                       hole_points_indices[triangle.third]});
+      visitor.after_subface_created(triangulated_polygons.back());
+    }
+
+    visitor.after_subface_creations();
+    return true;
+  }
+
+public:
+  template <typename Polygon,
+            typename PointRange,
+            typename PolygonRange,
+            typename NamedParameters>
+  bool operator()(const Polygon& polygon,
+                  const PointRange& points,
+                  PolygonRange& triangulated_polygons,
+                  const NamedParameters& np)
+  {
+    // PointMap
+    using PMap = typename GetPointMap<PointRange, NamedParameters>::const_type;
+    using Point_ref = typename boost::property_traits<PMap>::reference;
+
+    // Kernel
+    using Point = typename boost::property_traits<PMap>::value_type;
+    using Def_Kernel = typename CGAL::Kernel_traits<Point>::Kernel;
+    using Traits = typename internal_np::Lookup_named_param_def<
+                              internal_np::geom_traits_t,
+                              NamedParameters,
+                              Def_Kernel>::type;
+    using FT = typename Traits::FT;
+
+    // Visitor
+    using Visitor = typename internal_np::Lookup_named_param_def<
+                                            internal_np::visitor_t,
+                                            NamedParameters,
+                                            Triangulate_polygons::Default_visitor // default
+                                          >::type;
+
+    using parameters::choose_parameter;
+    using parameters::get_parameter;
+
+    PMap pm = choose_parameter<PMap>(get_parameter(np, internal_np::point_map));
+    Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
+    Visitor visitor = choose_parameter<Visitor>(get_parameter(np, internal_np::visitor),
+                                                Triangulate_polygons::Default_visitor());
+
+    typename Traits::Construct_cross_product_vector_3 cross_product =
+      traits.construct_cross_product_vector_3_object();
+
+    const std::size_t original_size = polygon.size();
+    if(original_size == 4)
+    {
+      Point_ref p0 = get(pm, points[polygon[0]]);
+      Point_ref p1 = get(pm, points[polygon[1]]);
+      Point_ref p2 = get(pm, points[polygon[2]]);
+      Point_ref p3 = get(pm, points[polygon[3]]);
+
+      /* Chooses the diagonal that will split the quad in two triangles that maximize
+       * the scalar product of the un-normalized normals of the two triangles.
+       * The lengths of the un-normalized normals (computed using cross-products of two vectors)
+       *  are proportional to the area of the triangles.
+       * Maximize the scalar product of the two normals will avoid skinny triangles,
+       * and will also taken into account the cosine of the angle between the two normals.
+       * In particular, if the two triangles are oriented in different directions,
+       * the scalar product will be negative.
+       */
+      visitor.before_subface_creations(polygon);
+
+      const FT p1p3 = cross_product(p2-p1, p3-p2) * cross_product(p0-p3, p1-p0);
+      const FT p0p2 = cross_product(p1-p0, p1-p2) * cross_product(p3-p2, p3-p0);
+      if(p0p2 > p1p3)
+      {
+        triangulated_polygons.push_back({polygon[0], polygon[1], polygon[2]});
+        triangulated_polygons.push_back({polygon[0], polygon[2], polygon[3]});
+      }
+      else
+      {
+        triangulated_polygons.push_back({polygon[0], polygon[1], polygon[3]});
+        triangulated_polygons.push_back({polygon[1], polygon[2], polygon[3]});
+      }
+
+      visitor.after_subface_created(triangulated_polygons[triangulated_polygons.size()-2]);
+      visitor.after_subface_created(triangulated_polygons[triangulated_polygons.size()-1]);
+
+      visitor.after_subface_creations();
+
+      return true;
+    }
+
+    return triangulate_polygon_with_hole_filling(polygon, points, triangulated_polygons, pm, visitor, np);
+  }
+}; // class Triangulate_polygon_soup_modifier
+
+} // namespace internal
+
+/**
+* \ingroup PMP_meshing_grp
+*
+* triangulates all polygons of a polygon soup. This function depends on the package \ref PkgTriangulation2.
+*
+* @tparam PointRange a model of `ConstRange`. The value type of its iterator is the point type.
+* @tparam PolygonRange a model of the concepts `SequenceContainer` and `Swappable`,
+*                      whose `value_type` is itself a model of the concept `SequenceContainer`
+*                      whose `value_type` is `std::size_t`.
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* @param points the point geometry of the soup to be triangulated
+* @param polygons the polygons to be triangulated
+* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{point_map}
+*     \cgalParamDescription{a property map associating points to the elements of the point set `points`}
+*     \cgalParamType{a model of `ReadablePropertyMap` whose key type is the value type
+*                    of the iterator of `PointRange` and whose value type is `geom_traits::Point_3`}
+*     \cgalParamDefault{`CGAL::Identity_property_map<geom_traits::Point_3>`}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{an instance of a geometric traits class}
+*     \cgalParamType{a class model of `Kernel`}
+*     \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
+*     \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{visitor}
+*     \cgalParamDescription{a visitor that enables to track how polygons are divided into triangles}
+*     \cgalParamType{a class model of `PMPTriangulateFaceVisitor` and `PMPHolefillingVisitor`}
+*     \cgalParamDefault{`Triangulate_polygons::Default_visitor`}
+*     \cgalParamExtra{Note that the visitor will be copied, so
+*                     it must not have any data member that does not have a reference-like type.}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+* This function calls `CGAL::Polygon_mesh_processing::triangulate_hole_polyline()` for each polygon.
+* Refer to its documentation for its named parameters.
+*
+* @pre No polygon within `polygons` is degenerate.
+*
+* @return `true` if all the polygons have been triangulated.
+*
+* @see `triangulate_faces()`
+*/
+template <typename PointRange,
+          typename PolygonRange,
+          typename NamedParameters = parameters::Default_named_parameters>
+bool triangulate_polygons(const PointRange& points,
+                          PolygonRange& polygons,
+                          const NamedParameters& np = parameters::default_values())
+{
+  using Polygon = typename boost::range_value<PolygonRange>::type;
+
+  PolygonRange triangulated_polygons;
+  triangulated_polygons.reserve(polygons.size());
+
+  bool success = true;
+
+  internal::Triangulate_polygon_soup_modifier modifier;
+  for(const Polygon& polygon : polygons)
+  {
+    if(polygon.size() <= 3)
+    {
+      triangulated_polygons.push_back(polygon);
+      continue;
+    }
+
+    if(!modifier(polygon, points, triangulated_polygons, np))
+      success = false;
+  }
+
+  std::swap(polygons, triangulated_polygons);
+
+  return success;
+}
+
+} // namespace Polygon_mesh_processing
+} // namespace CGAL
 
 #endif // CGAL_POLYGON_MESH_PROCESSING_TRIANGULATE_FACES_H

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h

@@ -729,11 +729,8 @@ namespace Polygon_mesh_processing {
     using parameters::get_parameter;
     using parameters::get_parameter_reference;
 
-    bool use_cdt =
+#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
-#ifdef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+    bool use_cdt = choose_parameter(get_parameter(np, internal_np::use_2d_constrained_delaunay_triangulation), false);
-      false;
-#else
-      choose_parameter(get_parameter(np, internal_np::use_2d_constrained_delaunay_triangulation), false);
 #endif
 bool use_dt3 =
 #ifdef CGAL_HOLE_FILLING_DO_NOT_USE_DT3

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -21,6 +21,7 @@ create_single_source_cgal_program("point_inside_surface_mesh_test.cpp")
 create_single_source_cgal_program("polygon_mesh_slicer_test.cpp")
 create_single_source_cgal_program("self_intersection_polyhedron_test.cpp")
 create_single_source_cgal_program("self_intersection_surface_mesh_test.cpp")
+create_single_source_cgal_program("self_intersection_triangle_soup_test.cpp")
 create_single_source_cgal_program("pmp_do_intersect_test.cpp")
 create_single_source_cgal_program("test_is_polygon_soup_a_polygon_mesh.cpp")
 create_single_source_cgal_program("test_stitching.cpp")

Polygon_mesh_processing/test/Polygon_mesh_processing/self_intersection_triangle_soup_test.cpp

@@ -0,0 +1,234 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+#include <CGAL/IO/polygon_soup_io.h>
+#include <CGAL/tags.h>
+#include <CGAL/Timer.h>
+
+#include <cstdlib>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel     EPICK;
+typedef CGAL::Exact_predicates_exact_constructions_kernel       EPECK;
+
+namespace PMP = ::CGAL::Polygon_mesh_processing;
+namespace CP = ::CGAL::parameters;
+
+template <typename K>
+int test_self_intersections(const std::string filename,
+                            const bool expected)
+{
+  std::vector<typename K::Point_3> points;
+  std::vector< std::array<std::size_t, 3> > triangles;
+
+  bool read_ok = CGAL::IO::read_polygon_soup(filename, points, triangles);
+
+  if ( !read_ok ) {
+    std::cerr << "Error: cannot read file: " << filename << std::endl;
+    return 1;
+  }
+
+  std::cout << "Reading file: " << filename << std::endl;
+
+  CGAL::Timer timer;
+  timer.start();
+
+  std::vector<std::pair<std::size_t, std::size_t> > intersected_tris;
+
+  PMP::triangle_soup_self_intersections<CGAL::Parallel_if_available_tag>(points, triangles, std::back_inserter(intersected_tris));
+  bool intersecting_1 = !intersected_tris.empty();
+
+  std::cout << "self_intersections test took " << timer.time() << " sec." << std::endl;
+  std::cout << intersected_tris.size() << " pairs of triangles are intersecting." << std::endl;
+
+  timer.reset();
+
+  bool intersecting_2 =
+    PMP::does_triangle_soup_self_intersect<CGAL::Parallel_if_available_tag>(points, triangles);
+
+  std::cout << "does_self_intersect test took " << timer.time() << " sec." << std::endl;
+  std::cout << (intersecting_2 ? "There is a self-intersection." :
+                                 "There are no self-intersections.") << std::endl;
+
+  assert(intersecting_1 == intersecting_2);
+  assert(intersecting_1 == expected);
+
+  std::cout << filename << " passed the tests." << std::endl << std::endl;
+
+  return 0;
+}
+
+template <typename K>
+int test_limited_self_intersections(const std::string& filename)
+{
+  std::vector<typename K::Point_3> points;
+  std::vector< std::array<std::size_t, 3> > triangles;
+
+  bool read_ok = CGAL::IO::read_polygon_soup(filename, points, triangles);
+
+  if ( !read_ok ) {
+    std::cerr << "Error: cannot read file: " << filename << std::endl;
+    return 1;
+  }
+
+  CGAL::Timer timer;
+  timer.start();
+
+  std::vector<std::pair<std::size_t, std::size_t> > intersected_tris;
+
+#ifdef CGAL_LINKED_WITH_TBB
+  PMP::triangle_soup_self_intersections()<CGAL::Parallel_if_available_tag>(
+        points, triangles,
+        std::back_inserter(intersected_tris), CGAL::parameters::maximum_number(40));
+  std::cout << "self_intersections test for 40 SI took " << timer.time() << " sec." << std::endl;
+  std::cout << "Found " << intersected_tris.size() << " SIs." << std::endl;
+  if(intersected_tris.size() < 40)
+  {
+    std::cerr<<"Not enough intersections found in parallel."<<std::endl;
+    return 1;
+  }
+  intersected_tris.clear();
+  timer.reset();
+#endif
+  PMP::triangle_soup_self_intersections<CGAL::Sequential_tag>(
+    points, triangles,
+    std::back_inserter(intersected_tris), CGAL::parameters::maximum_number(40));
+  std::cout << "self_intersections test for 40 SI took " << timer.time() << " sec." << std::endl;
+  timer.reset();
+  if(intersected_tris.size() != 40)
+  {
+    std::cerr<<"Too many intersections found in sequential"<<std::endl;
+    return 1;
+  }
+  return 0;
+}
+
+int main(int argc, char** argv)
+{
+  // If file(s) are provided, the associated expected result must also be provided.
+  // Note that this expected value is a Boolean that is passed in command line
+  // with either 'true' or 'false' (and not integers), that is for example:
+  // > self_intersection_surface_mesh_test data/U.off false
+
+  // First test ----------------------------------------------------------------
+  bool expected = false;
+  std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off");
+  if(argc > 1) {
+    assert(argc > 2);
+    std::stringstream ss(argv[2]);
+    ss >> std::boolalpha >> expected;
+    assert(!ss.fail()); // make sure that argv[2] is either 'true' or 'false'
+  }
+
+  std::cout << "First test (EPICK):" << std::endl;
+  int r = test_self_intersections<EPICK>(filename, expected);
+
+  std::cout << "First test (EPECK):" << std::endl;
+  r += test_self_intersections<EPECK>(filename, expected);
+
+  // Second test ---------------------------------------------------------------
+  expected = true;
+  filename = (argc > 3) ? argv[3] : CGAL::data_file_path("meshes/mannequin-devil.off");
+  if(argc > 3) {
+    assert(argc > 4);
+    std::stringstream ss(argv[4]);
+    ss >> std::boolalpha >> expected;
+    assert(!ss.fail());
+  }
+
+  std::cout << "Second test (EPICK):" << std::endl;
+  r += test_self_intersections<EPICK>(filename, expected);
+
+  std::cout << "Second test (EPECK):" << std::endl;
+  r += test_self_intersections<EPECK>(filename, expected);
+
+  // Third test ----------------------------------------------------------------
+  expected = true;
+  filename = (argc > 5) ? argv[5] : "data/overlapping_triangles.off";
+  if(argc > 5) {
+    assert(argc > 6);
+    std::stringstream ss(argv[6]);
+    ss >> std::boolalpha >> expected;
+    assert(!ss.fail());
+  }
+
+  std::cout << "Third test (EPICK):" << std::endl;
+  r += test_self_intersections<EPICK>(filename, expected);
+
+  std::cout << "Third test (EPECK):" << std::endl;
+  r += test_self_intersections<EPECK>(filename, expected);
+
+  // Fourth test ----------------------------------------------------------------
+  expected = true;
+  filename = (argc > 7) ? argv[7] : "data_degeneracies/degtri_single.off";
+  if(argc > 7) {
+    assert(argc > 8);
+    std::stringstream ss(argv[8]);
+    ss >> std::boolalpha >> expected;
+    assert(!ss.fail());
+  }
+
+  std::cout << "Fourth test (EPICK):" << std::endl;
+  r += test_self_intersections<EPICK>(filename, expected);
+
+  std::cout << "Fourth test (EPECK):" << std::endl;
+  r += test_self_intersections<EPECK>(filename, expected);
+
+  filename = (argc > 9) ? argv[9] : CGAL::data_file_path("meshes/mannequin-devil.off");
+
+  std::cout << "Test with maximum_number (EPICK):" << std::endl;
+  r += test_limited_self_intersections<EPICK>(filename);
+
+  std::cout << "Test with maximum_number (EPECK):" << std::endl;
+  r += test_limited_self_intersections<EPECK>(filename);
+
+  // extra hand written tests
+  {
+    // shared edge with same point: no self-intersection
+    typedef EPICK::Point_3 Point_3;
+    std::vector<EPICK::Point_3> points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,-1,0)};
+    std::vector< std::array<std::size_t, 3> > triangles={{0,1,2}, {0,1,3}};
+    assert(!PMP::does_triangle_soup_self_intersect(points, triangles));
+  }
+  {
+    // shared edge with duplicated points: self-intersection
+    typedef EPICK::Point_3 Point_3;
+    std::vector<EPICK::Point_3> points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0),Point_3(0,0,0), Point_3(1,0,0), Point_3(0,-1,0)};
+    std::vector< std::array<std::size_t, 3> > triangles={{0,1,2}, {3,4,5}};
+    assert(PMP::does_triangle_soup_self_intersect(points, triangles));
+  }
+  {
+    // shared vertex with same point: no self-intersection
+    typedef EPICK::Point_3 Point_3;
+    std::vector<EPICK::Point_3> points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,2,0), Point_3(1,2,0)};
+    std::vector< std::array<std::size_t, 3> > triangles={{0,1,2}, {3,4,2}};
+    assert(!PMP::does_triangle_soup_self_intersect(points, triangles));
+  }
+  {
+    // shared vertex with duplicated points: self-intersection
+    typedef EPICK::Point_3 Point_3;
+    std::vector<EPICK::Point_3> points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,2,0), Point_3(1,2,0), Point_3(0,1,0)};
+    std::vector< std::array<std::size_t, 3> > triangles={{0,1,2}, {3,4,5}};
+    assert(PMP::does_triangle_soup_self_intersect(points, triangles));
+  }
+  {
+    // 4 triangles around a shared edge: no self-intersection
+    typedef EPICK::Point_3 Point_3;
+    std::vector<EPICK::Point_3> points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,-1,0), Point_3(0,0,1), Point_3(0,0,-1)};
+    std::vector< std::array<std::size_t, 3> > triangles={{0,1,2},{0,1,3},{0,1,4},{0,1,5}};
+    assert(!PMP::does_triangle_soup_self_intersect(points, triangles));
+  }
+  {
+    // 4 triangles around a shared edge but two triangles intersecting: self-intersection
+    typedef EPICK::Point_3 Point_3;
+    std::vector<EPICK::Point_3> points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,0.5,0), Point_3(0,0,1), Point_3(0,0,-1)};
+    std::vector< std::array<std::size_t, 3> > triangles={{0,1,2},{0,1,3},{0,1,4},{0,1,5}};
+    assert(PMP::does_triangle_soup_self_intersect(points, triangles));
+  }
+
+  return r;
+}

Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_faces_test.cpp

@@ -1,13 +1,14 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-
 #include <CGAL/Surface_mesh.h>
 
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+
 #include <CGAL/boost/graph/copy_face_graph.h>
 #include <CGAL/boost/graph/Dual.h>
 #include <CGAL/boost/graph/named_params_helper.h>
 #include <CGAL/centroid.h>
-#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/Polygon_mesh_processing/polygon_mesh_to_polygon_soup.h>
 
 #include <boost/graph/filtered_graph.hpp>
 
@@ -20,6 +21,8 @@ template <typename K>
 bool
 test_triangulate_faces()
 {
+  std::cout << "\n--- test_triangulate_faces(" << typeid(K).name() << ") ---" << std::endl;
+
   typedef typename K::Point_3                    Point;
   typedef CGAL::Surface_mesh<Point>          Surface_mesh;
 
@@ -42,6 +45,8 @@ template <typename K>
 bool
 test_triangulate_faces_with_named_parameters()
 {
+  std::cout << "\n--- test_triangulate_faces_with_named_parameters(" << typeid(K).name() << ") ---" << std::endl;
+
   typedef typename K::Point_3                    Point;
   typedef CGAL::Surface_mesh<Epic::Point_3>      Surface_mesh;
 
@@ -77,13 +82,15 @@ test_triangulate_faces_with_named_parameters()
 
 template <typename K>
 bool
-test_triangulate_face_range()
+test_triangulate_face_range(const std::string& filename)
 {
+  std::cout << "\n--- test_triangulate_face_range(" << typeid(K).name() << ") ---" << std::endl;
+
   typedef typename K::Point_3                    Point;
   typedef CGAL::Surface_mesh<Point>          Surface_mesh;
 
   Surface_mesh mesh;
-  std::ifstream input(CGAL::data_file_path("meshes/cube_quad.off"));
+  std::ifstream input(filename);
 
   if (!input || !(input >> mesh) || mesh.is_empty())
   {
@@ -92,6 +99,18 @@ test_triangulate_face_range()
   }
 
   bool success = CGAL::Polygon_mesh_processing::triangulate_faces(faces(mesh), mesh);
+
+  for(auto f : faces(mesh))
+  {
+    if(!is_triangle(halfedge(f, mesh), mesh))
+    {
+      std::cout << "non triangular face:" << std::endl;
+      for(auto h : halfedges_around_face(halfedge(f, mesh), mesh))
+        std::cout << "  " << mesh.point(target(h, mesh)) << std::endl;
+      assert(false);
+    }
+  }
+
   assert(CGAL::is_triangle_mesh(mesh));
 
   // For compilation
@@ -105,6 +124,8 @@ template <typename K>
 bool
 test_triangulate_face()
 {
+  std::cout << "\n--- test_triangulate_face(" << typeid(K).name() << ") ---" << std::endl;
+
   typedef typename K::Point_3                    Point;
   typedef CGAL::Surface_mesh<Point>          Surface_mesh;
 
@@ -117,21 +138,17 @@ test_triangulate_face()
     return false;
   }
 
-  unsigned int nb = 0;
   for(typename boost::graph_traits<Surface_mesh>::face_descriptor fit : faces(mesh))
   {
-    if (nb > 4)
+    if (next(next(halfedge(fit, mesh), mesh), mesh) !=   prev(halfedge(fit, mesh), mesh))
-      break;
-    else if (next(next(halfedge(fit, mesh), mesh), mesh)
-             !=   prev(halfedge(fit, mesh), mesh))
     {
-      if(CGAL::Polygon_mesh_processing::triangulate_face(fit, mesh))
+      if(!CGAL::Polygon_mesh_processing::triangulate_face(fit, mesh))
-        ++nb;
-      else
         assert(false);
     }
   }
 
+  assert(CGAL::is_triangle_mesh(mesh));
+
   return true;
 }
 
@@ -139,6 +156,8 @@ template <typename K>
 bool
 test_triangulate_triangle_face()
 {
+  std::cout << "\n--- test_triangulate_triangle_face(" << typeid(K).name() << ") ---" << std::endl;
+
   typedef typename K::Point_3                    Point;
   typedef CGAL::Surface_mesh<Point>          Surface_mesh;
 
@@ -156,6 +175,9 @@ test_triangulate_triangle_face()
     if(!CGAL::Polygon_mesh_processing::triangulate_face(fit, mesh, CGAL::parameters::geom_traits(K())))
       assert(false);
   }
+
+  assert(CGAL::is_triangle_mesh(mesh));
+
   return true;
 }
 
@@ -211,6 +233,8 @@ template <typename K>
 bool
 test_dual_with_various_faces()
 {
+  std::cout << "\n--- test_dual_with_various_faces(" << typeid(K).name() << ") ---" << std::endl;
+
   typedef typename K::Point_3                    Point;
   typedef CGAL::Surface_mesh<Point>          Surface_mesh;
 
@@ -235,27 +259,85 @@ test_dual_with_various_faces()
 
   for(typename boost::graph_traits<Surface_mesh>::face_descriptor fit : faces(sm_dual))
   {
-    if(!CGAL::Polygon_mesh_processing::triangulate_face(fit, sm_dual))
+    if(!CGAL::Polygon_mesh_processing::triangulate_face(fit, sm_dual,
+                                                        CGAL::parameters::use_2d_constrained_delaunay_triangulation(true)))
       assert(false);
   }
+
+  assert(CGAL::is_triangle_mesh(sm_dual));
+
   return true;
 }
 
-int main()
+template <typename K>
+bool
+test_triangulate_soup()
 {
+  std::cout << "\n--- test_triangulate_soup(" << typeid(K).name() << ") ---" << std::endl;
+
+  typedef typename K::Point_3                    Point;
+  typedef CGAL::Surface_mesh<Point>          Surface_mesh;
+
+  Surface_mesh mesh;
+  std::ifstream input(CGAL::data_file_path("meshes/elephant.off"));
+
+  if (!input || !(input >> mesh) || mesh.is_empty())
+  {
+    std::cerr << "Not a valid off file." << std::endl;
+    return false;
+  }
+
+  typedef typename boost::property_map<Surface_mesh, boost::vertex_point_t>::type Pmap;
+  Pmap vpmap = get_property_map(boost::vertex_point, mesh);
+
+  CGAL::Dual<Surface_mesh> dual(mesh);
+  // copy dual to a sm
+  Surface_mesh sm_dual;
+  CGAL::copy_face_graph(dual, sm_dual,
+                        CGAL::parameters::vertex_point_map(
+                          Dual_vpm<Surface_mesh, Point, Pmap>(mesh, vpmap)));
+
+  std::vector<Point> points;
+  std::vector<std::vector<std::size_t> > polygons;
+  CGAL::Polygon_mesh_processing::polygon_mesh_to_polygon_soup(sm_dual, points, polygons);
+
+  bool success = CGAL::Polygon_mesh_processing::triangulate_polygons(points, polygons,
+                                                                     CGAL::parameters::geom_traits(K())
+                                                                                      .use_2d_constrained_delaunay_triangulation(false));
+  for(std::size_t i = 0; i < polygons.size(); ++i)
+  {
+    assert(polygons[i].size() == 3);
+  }
+
+  // For compilation
+  success = CGAL::Polygon_mesh_processing::triangulate_polygons(points, polygons);
+
+  return success;
+}
+
+int main(int argc, char** argv)
+{
+  if(argc > 1)
+  {
+    assert(test_triangulate_face_range<Epic>(argv[1]));
+  }
+
   assert(test_triangulate_faces<Epic>());
   assert(test_triangulate_faces_with_named_parameters<Epic>());
-  assert(test_triangulate_face_range<Epic>());
+  assert(test_triangulate_face_range<Epic>(CGAL::data_file_path("meshes/cube_quad.off")));
   assert(test_triangulate_face<Epic>());
   assert(test_triangulate_triangle_face<Epic>());
   assert(test_dual_with_various_faces<Epic>());
+  assert(test_triangulate_soup<Epic>());
 
   assert(test_triangulate_faces<Epec>());
   assert(test_triangulate_faces_with_named_parameters<Epec>());
-  assert(test_triangulate_face_range<Epec>());
+  assert(test_triangulate_face_range<Epec>(CGAL::data_file_path("meshes/cube_quad.off")));
   assert(test_triangulate_face<Epec>());
   assert(test_triangulate_triangle_face<Epec>());
   assert(test_dual_with_various_faces<Epec>());
+  assert(test_triangulate_soup<Epec>());
 
+  std::cout << "Done" << std::endl;
   return EXIT_SUCCESS;
 }

Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt

@@ -105,7 +105,7 @@ add_custom_target(self_intersection_plugin)
 add_dependencies(self_intersection_plugin selection_plugin)
 
 polyhedron_demo_plugin(triangulate_facets_plugin Triangulate_facets_plugin KEYWORDS PMP)
-target_link_libraries(triangulate_facets_plugin PUBLIC scene_surface_mesh_item scene_selection_item)
+target_link_libraries(triangulate_facets_plugin PUBLIC scene_surface_mesh_item scene_selection_item scene_polygon_soup_item)
 
 polyhedron_demo_plugin(corefinement_plugin Corefinement_plugin KEYWORDS PMP)
 target_link_libraries(corefinement_plugin PUBLIC scene_surface_mesh_item)

Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp

@@ -178,12 +178,13 @@ class Polyhedron_demo_isotropic_remeshing_plugin :
   typedef std::unordered_set<edge_descriptor>    Edge_set;
   typedef Scene_polyhedron_selection_item::Is_constrained_map<Edge_set> Edge_constrained_pmap;
 
-  struct Visitor
+  struct Selection_updater_visitor
+    : public CGAL::Polygon_mesh_processing::Hole_filling::Default_visitor
   {
     typedef typename Scene_polyhedron_selection_item::Selection_set_facet Container;
     Container& faces;
 
-    Visitor(Container& container)
+    Selection_updater_visitor(Container& container)
       : faces(container)
     {}
 
@@ -493,7 +494,7 @@ public Q_SLOTS:
                        (QMessageBox::Ok | QMessageBox::Cancel),
                        QMessageBox::Ok))
                   {
-                    Visitor visitor(selection_item->selected_facets);
+                    Selection_updater_visitor visitor(selection_item->selected_facets);
                     CGAL::Polygon_mesh_processing::triangulate_faces(selection_item->selected_facets,
                       pmesh,
                       CGAL::parameters::visitor(visitor));

Polyhedron/demo/Polyhedron/Plugins/PMP/Triangulate_facets_plugin.cpp

@@ -4,10 +4,13 @@
 #include "Messages_interface.h"
 #include <CGAL/Three/Polyhedron_demo_plugin_helper.h>
 #include <CGAL/Three/Three.h>
+
 #include "Scene_surface_mesh_item.h"
 #include "Scene_polyhedron_selection_item.h"
+#include "Scene_polygon_soup_item.h"
 
 #include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+
 using namespace CGAL::Three;
 class Polyhedron_demo_triangulate_facets_plugin :
   public QObject,
@@ -20,12 +23,13 @@ class Polyhedron_demo_triangulate_facets_plugin :
   typedef Scene_surface_mesh_item::Face_graph FaceGraph;
   typedef boost::graph_traits<FaceGraph>::face_descriptor face_descriptor;
 
-  struct Visitor
+  struct Selection_updater_visitor
+    : public CGAL::Polygon_mesh_processing::Hole_filling::Default_visitor
   {
     typedef typename Scene_polyhedron_selection_item::Selection_set_facet Container;
     Container& faces;
 
-    Visitor(Container& container)
+    Selection_updater_visitor(Container& container)
       : faces(container)
     {}
     void before_subface_creations(face_descriptor fd)
@@ -66,6 +70,8 @@ public:
         return true;
       if ( qobject_cast<Scene_polyhedron_selection_item*>(scene->item(index)))
         return true;
+      if(qobject_cast<Scene_polygon_soup_item*>(scene->item(index)))
+        return true;
     }
     return false;
   }
@@ -82,39 +88,54 @@ public Q_SLOTS:
       Scene_polyhedron_selection_item* selection_item =
         qobject_cast<Scene_polyhedron_selection_item*>(scene->item(index));
 
-      SMesh* pMesh = (sm_item != nullptr)
+      Scene_polygon_soup_item* soup_item =
-                    ? sm_item->polyhedron()
+        qobject_cast<Scene_polygon_soup_item*>(scene->item(index));
-                    : selection_item->polyhedron();
 
-      if(!pMesh) continue;
+      if (soup_item)
-      if(is_triangle_mesh(*pMesh)) {
-        CGAL::Three::Three::warning(tr("The polyhedron  \"%1\"  is already triangulated.")
-                          .arg(sm_item->name()) );
-        continue;
-      }
-      if (sm_item)
       {
-        if (!CGAL::Polygon_mesh_processing::triangulate_faces(*pMesh))
+        soup_item->triangulate();
-          CGAL::Three::Three::warning(tr("Some facets could not be triangulated."));
       }
-      else if (selection_item)
+      else
       {
-        Visitor visitor(selection_item->selected_facets);
+        SMesh* pMesh = (sm_item != nullptr) ? sm_item->polyhedron()
-        if (!CGAL::Polygon_mesh_processing::triangulate_faces(
+                                            : selection_item->polyhedron();
-                 selection_item->selected_facets,
-                 *pMesh,
-                 CGAL::parameters::visitor(visitor)))
-          CGAL::Three::Three::warning(tr("Some facets could not be triangulated."));
 
-        sm_item = selection_item->polyhedron_item();
+        if(!pMesh)
-        selection_item->set_num_faces(num_faces(*sm_item->face_graph()));
+          continue;
 
-        selection_item->invalidateOpenGLBuffers();
+        if(is_triangle_mesh(*pMesh))
-        selection_item->itemChanged();
+        {
+          CGAL::Three::Three::warning(tr("The polyhedron  \"%1\"  is already triangulated.")
+                            .arg(sm_item->name()) );
+          continue;
+        }
+
+        if (sm_item)
+        {
+          if (!CGAL::Polygon_mesh_processing::triangulate_faces(*pMesh))
+            CGAL::Three::Three::warning(tr("Some facets could not be triangulated."));
+
+          sm_item->invalidateOpenGLBuffers();
+        }
+        else if (selection_item)
+        {
+          Selection_updater_visitor visitor(selection_item->selected_facets);
+          if (!CGAL::Polygon_mesh_processing::triangulate_faces(
+                  selection_item->selected_facets,
+                  *pMesh,
+                  CGAL::parameters::visitor(visitor)))
+            CGAL::Three::Three::warning(tr("Some facets could not be triangulated."));
+
+          sm_item = selection_item->polyhedron_item();
+          selection_item->set_num_faces(num_faces(*sm_item->face_graph()));
+
+          selection_item->invalidateOpenGLBuffers();
+          selection_item->itemChanged();
+        }
+
+        sm_item->resetColors(); // @todo should have a visitor to give the color of the parent face
       }
 
-      sm_item->resetColors();
-      sm_item->invalidateOpenGLBuffers();
       scene->itemChanged(sm_item);
     } // end of the loop on the selected items
 

Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.cpp

@@ -22,6 +22,7 @@
 #include <CGAL/Polygon_mesh_processing/orientation.h>
 #include <CGAL/Polygon_mesh_processing/repair.h>
 #include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
 #include <CGAL/Polygon_2.h>
 #include <CGAL/version.h>
 
@@ -419,6 +420,32 @@ void Scene_polygon_soup_item::inside_out()
   invalidateOpenGLBuffers();
 }
 
+void Scene_polygon_soup_item::repair(bool erase_dup, bool req_same_orientation)
+{
+  QApplication::setOverrideCursor(Qt::BusyCursor);
+  CGAL::Polygon_mesh_processing::repair_polygon_soup(
+        d->soup->points,
+        d->soup->polygons,
+        CGAL::parameters::erase_all_duplicates(erase_dup)
+                         .require_same_orientation(req_same_orientation));
+  QApplication::restoreOverrideCursor();
+  invalidateOpenGLBuffers();
+}
+
+bool Scene_polygon_soup_item::triangulate()
+{
+  QApplication::setOverrideCursor(Qt::BusyCursor);
+
+  bool success = true;
+
+  CGAL::Polygon_mesh_processing::triangulate_polygons(d->soup->points, d->soup->polygons);
+
+  QApplication::restoreOverrideCursor();
+  invalidateOpenGLBuffers();
+
+  return success;
+}
+
 bool
 Scene_polygon_soup_item::orient(std::vector<std::size_t>& non_manifold_vertices)
 {
@@ -894,20 +921,6 @@ void Scene_polygon_soup_item::computeElements() const
   QApplication::restoreOverrideCursor();
 }
 
-void Scene_polygon_soup_item::repair(bool erase_dup, bool req_same_orientation)
-{
-  QApplication::setOverrideCursor(Qt::BusyCursor);
-  CGAL::Polygon_mesh_processing::repair_polygon_soup(
-        d->soup->points,
-        d->soup->polygons,
-        CGAL::parameters::
-        erase_all_duplicates(erase_dup)
-        .require_same_orientation(req_same_orientation));
-  QApplication::restoreOverrideCursor();
-
- // CGAL::Three::Three::information(
-}
-
 CGAL::Three::Scene_item::Header_data Scene_polygon_soup_item::header() const
 {
   CGAL::Three::Scene_item::Header_data data;

Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.h

@@ -187,6 +187,7 @@ public Q_SLOTS:
     bool exportAsSurfaceMesh(SMesh*);
     void inside_out();
     void repair(bool erase_dup, bool req_same_orientation);
+    bool triangulate();
 
     void setDisplayNonManifoldEdges(const bool);
     bool displayNonManifoldEdges() const;

Principal_component_analysis/include/CGAL/PCA_util.h

@@ -391,7 +391,7 @@ assemble_covariance_matrix_3(InputIterator first,
                    0.0, radius, 0.0,
                    0.0, 0.0, radius};
     Matrix transformation = init_matrix<FT>(3,delta);
-    FT volume = (FT)(4.0/3.0) * radius * t.squared_radius();
+    FT volume = radius * t.squared_radius();
 
                 // skip zero measure primitives
     if(volume == (FT)0.0)
@@ -400,8 +400,9 @@ assemble_covariance_matrix_3(InputIterator first,
     // Find the 2nd order moment for the sphere wrt to the origin by an affine transformation.
 
     // Transform the standard 2nd order moment using the transformation matrix
-    transformation = (3.0/4.0) * volume * transformation * moment * LA::transpose(transformation);
+    transformation = volume * transformation * moment * LA::transpose(transformation);
 
+    volume *= FT(4.0/3.0);
     // Translate the 2nd order moment to the center of the sphere.
     FT x0 = t.center().x();
     FT y0 = t.center().y();
@@ -476,7 +477,7 @@ assemble_covariance_matrix_3(InputIterator first,
                    0.0,    radius, 0.0,
                    0.0,    0.0,    radius};
     Matrix transformation = init_matrix<FT>(3,delta);
-    FT area = (FT)4.0 * t.squared_radius();
+    FT area = t.squared_radius();
 
                 // skip zero measure primitives
     if(area == (FT)0.0)
@@ -485,8 +486,9 @@ assemble_covariance_matrix_3(InputIterator first,
     // Find the 2nd order moment for the sphere wrt to the origin by an affine transformation.
 
     // Transform the standard 2nd order moment using the transformation matrix
-    transformation = (1.0/4.0) * area * transformation * moment * LA::transpose(transformation);
+    transformation = area * transformation * moment * LA::transpose(transformation);
 
+    area *= FT(4.0);
     // Translate the 2nd order moment to the center of the sphere.
     FT x0 = t.center().x();
     FT y0 = t.center().y();

Principal_component_analysis/include/CGAL/PCA_util_Eigen.h

@@ -368,7 +368,7 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation << radius, 0.0, 0.0,
                       0.0, radius, 0.0,
                       0.0, 0.0, radius;
-    FT volume = (FT)(4.0/3.0) * radius * t.squared_radius();
+    FT volume =  radius * t.squared_radius();
 
     // skip zero measure primitives
     if(volume == (FT)0.0)
@@ -377,8 +377,9 @@ assemble_covariance_matrix_3(InputIterator first,
     // Find the 2nd order moment for the sphere wrt to the origin by an affine transformation.
 
     // Transform the standard 2nd order moment using the transformation matrix
-    transformation = (3.0/4.0) * volume * transformation * moment * transformation.transpose();
+    transformation = volume * transformation * moment * transformation.transpose();
 
+    volume *= FT(4.0 / 3.0);
     // Translate the 2nd order moment to the center of the sphere.
     FT x0 = t.center().x();
     FT y0 = t.center().y();
@@ -453,7 +454,7 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation << radius, 0.0,    0.0,
                       0.0,    radius, 0.0,
                       0.0,    0.0,    radius;
-    FT area = (FT)4.0 * t.squared_radius();
+    FT area = t.squared_radius();
 
     // skip zero measure primitives
     if(area == (FT)0.0)
@@ -462,8 +463,9 @@ assemble_covariance_matrix_3(InputIterator first,
     // Find the 2nd order moment for the sphere wrt to the origin by an affine transformation.
 
     // Transform the standard 2nd order moment using the transformation matrix
-    transformation = (1.0/4.0) * area * transformation * moment * transformation.transpose();
+    transformation = area * transformation * moment * transformation.transpose();
 
+    area *= FT(4.0);
     // Translate the 2nd order moment to the center of the sphere.
     FT x0 = t.center().x();
     FT y0 = t.center().y();