diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md
index 5d3bbdd0cf8..a8cd52195eb 100644
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@@ -40,6 +40,8 @@ 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.
 
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
index c347d4a7ae3..b17c41399d0 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
+++ b/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
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp
index f46d09e0f2b..9b246f7678d 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp
+++ b/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;
 
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polygon_mesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polygon_mesh.h
index ab837b05c88..7b36d3295ec 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polygon_mesh.h
+++ b/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;
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polyline.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polyline.h
index 797a6912f26..5e84ba9e694 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polyline.h
+++ b/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) {
-    const std::size_t ip = (i + 1) % size;
-    if (vertices[i] != vertices[ip]) {
-      cdt.insert_constraint(vertices[i], vertices[ip]);
+  try
+  {
+    for (std::size_t i = 0; i < size; ++i) {
+      const std::size_t ip = (i + 1) % size;
+      if (vertices[i] != vertices[ip]) {
+        cdt.insert_constraint(vertices[i], vertices[ip]);
+      }
     }
   }
+  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(),
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
index cb24ce8ae7c..bca94855450 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
+++ b/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"
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_faces.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_faces.h
index d7b88034d60..1e9135658db 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_faces.h
+++ b/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
-          , typename VertexPointMap
-          , typename Kernel
-          , typename Visitor>
-class Triangulate_modifier
+template <typename PolygonMesh>
+class Triangulate_polygon_mesh_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;
-  typedef typename boost::graph_traits<PM>::halfedge_descriptor halfedge_descriptor;
-  typedef typename boost::graph_traits<PM>::face_descriptor face_descriptor;
-  typedef typename boost::graph_traits<PM>::edge_descriptor edge_descriptor;
-  typedef typename Kernel::Point_3 Point;
-
-  struct Face_info {
-    typename boost::graph_traits<PM>::halfedge_descriptor e[3];
-    bool is_external;
-  };
-
-  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)
+private:
+  template <typename VPM,
+            typename Visitor,
+            typename NamedParameters>
+  bool triangulate_face_with_hole_filling(face_descriptor f,
+                                          PolygonMesh& pmesh,
+                                          const VPM vpm,
+                                          Visitor visitor,
+                                          const NamedParameters& np)
   {
     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),
-                                   parameters::geom_traits(_traits));
+    PMP::triangulate_hole_polyline(hole_points, std::back_inserter(patch), np);
 
     if(patch.empty())
       return false;
 
     // triangulate the hole
-    std::map< std::pair<int, int> , halfedge_descriptor > halfedge_map;
-    int i=0;
+    std::map<std::pair<int, int>, halfedge_descriptor > halfedge_map;
+    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)
-        first=false;
+      if(first)
+        first = false;
       else
-        f=add_face(pmesh);
+        f = add_face(pmesh);
+
       visitor.after_subface_created(f);
 
-      std::array<int, 4> indices =
-        make_array( triangle.first,
-                    triangle.second,
-                    triangle.third,
-                    triangle.first );
+      std::array<int, 4> indices = make_array(triangle.first,
+                                              triangle.second,
+                                              triangle.third,
+                                              triangle.first);
       for (int i=0; i<3; ++i)
       {
         typename std::map< std::pair<int, int> , halfedge_descriptor >::iterator insert_res =
-          halfedge_map.insert(
-            std::make_pair( std::make_pair(indices[i], indices[i+1]),
-                            boost::graph_traits<PM>::null_halfedge() ) ).first;
-        if (insert_res->second == boost::graph_traits<PM>::null_halfedge())
+          halfedge_map.emplace(std::make_pair(indices[i], indices[i+1]),
+                               boost::graph_traits<PolygonMesh>::null_halfedge()).first;
+        if(insert_res->second == boost::graph_traits<PolygonMesh>::null_halfedge())
         {
           halfedge_descriptor nh = halfedge(add_edge(pmesh), pmesh);
-          insert_res->second=nh;
-          halfedge_map[std::make_pair(indices[i+1], indices[i])]=opposite(nh, pmesh);
+          insert_res->second = nh;
+          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>
-  bool operator()(FaceRange face_range, PM& pmesh, bool use_cdt, Visitor visitor)
+public:
+  template <typename NamedParameters>
+  bool operator()(face_descriptor f,
+                  PolygonMesh& pmesh,
+                  const NamedParameters& np)
   {
-   bool result = true;
-    // One need to store facet handles into a vector, because the list of
-    // 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)));
+    using Traits = typename GetGeomTraits<PolygonMesh, NamedParameters>::type;
+    using VPM = typename GetVertexPointMap<PolygonMesh, NamedParameters>::type;
 
-    //only consider non-triangular faces
-    for(face_descriptor fit : face_range)
-      if ( next( next( halfedge(fit, pmesh), pmesh), pmesh)
-        !=       prev( halfedge(fit, pmesh), pmesh) )
-        facets.push_back(fit);
+    using FT = typename Traits::FT;
+    using Point_ref = typename boost::property_traits<VPM>::reference;
 
-    // Iterates on the vector of face descriptors
-    for(face_descriptor f : facets)
+    using Visitor = typename internal_np::Lookup_named_param_def<
+                                            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))
-       result = false;
+      halfedge_descriptor v0, v1, v2, v3;
+      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)
-  {
-    //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
+} // 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;
-  using parameters::get_parameter;
-
-  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);
+  internal::Triangulate_polygon_mesh_modifier<PolygonMesh> modifier;
+  return modifier(f, pmesh, np);
 }
 
 /**
@@ -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 parameters::get_parameter;
+  using face_descriptor = typename boost::graph_traits<PolygonMesh>::face_descriptor;
 
-  //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));
+  bool result = true;
 
-  //Kernel
-  typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type Kernel;
-  Kernel traits = choose_parameter<Kernel>(get_parameter(np, internal_np::geom_traits));
+  // One needs to store the facets into a vector, because the list of
+  // facets of the polyhedron will be modified during the loop, and
+  // that invalidates the range [facets_begin(), facets_end()[.
+  std::vector<face_descriptor> facets(std::begin(face_range), std::end(face_range));
 
-  //Option
-  bool use_cdt = choose_parameter(get_parameter(np, internal_np::use_delaunay_triangulation), true);
+  internal::Triangulate_polygon_mesh_modifier<PolygonMesh> modifier;
+  for(face_descriptor f : facets)
+  {
+    if(!modifier(f, pmesh, np))
+      result = false;
+  }
 
-  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(face_range, pmesh, use_cdt, visitor);
+  return result;
 }
 
 /**
@@ -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
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h
index 46c3a6904be..ba55ceca89e 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h
+++ b/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 =
-#ifdef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
-      false;
-#else
-      choose_parameter(get_parameter(np, internal_np::use_2d_constrained_delaunay_triangulation), false);
+#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2
+    bool use_cdt = 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
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_faces_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_faces_test.cpp
index 5f239c7dfaa..66085b83c74 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_faces_test.cpp
+++ b/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)
-      break;
-    else if (next(next(halfedge(fit, mesh), mesh), mesh)
-             !=   prev(halfedge(fit, mesh), mesh))
+    if (next(next(halfedge(fit, mesh), mesh), mesh) !=   prev(halfedge(fit, mesh), mesh))
     {
-      if(CGAL::Polygon_mesh_processing::triangulate_face(fit, mesh))
-        ++nb;
-      else
+      if(!CGAL::Polygon_mesh_processing::triangulate_face(fit, mesh))
         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;
 }
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt b/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt
index 6ff2b7c356d..09cf3baaa08 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt
+++ b/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)
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp
index a15e1f14dfd..37dcedc8adb 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp
+++ b/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));
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Triangulate_facets_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Triangulate_facets_plugin.cpp
index 3f86836f4bf..6b85c6fb45b 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Triangulate_facets_plugin.cpp
+++ b/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)
-                    ? sm_item->polyhedron()
-                    : selection_item->polyhedron();
+      Scene_polygon_soup_item* soup_item =
+        qobject_cast<Scene_polygon_soup_item*>(scene->item(index));
 
-      if(!pMesh) continue;
-      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 (soup_item)
       {
-        if (!CGAL::Polygon_mesh_processing::triangulate_faces(*pMesh))
-          CGAL::Three::Three::warning(tr("Some facets could not be triangulated."));
+        soup_item->triangulate();
       }
-      else if (selection_item)
+      else
       {
-        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."));
+        SMesh* pMesh = (sm_item != nullptr) ? sm_item->polyhedron()
+                                            : selection_item->polyhedron();
 
-        sm_item = selection_item->polyhedron_item();
-        selection_item->set_num_faces(num_faces(*sm_item->face_graph()));
+        if(!pMesh)
+          continue;
 
-        selection_item->invalidateOpenGLBuffers();
-        selection_item->itemChanged();
+        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))
+            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
 
diff --git a/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.cpp b/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.cpp
index b258e0c5714..8ba8813f349 100644
--- a/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.cpp
+++ b/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;
diff --git a/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.h b/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.h
index 5dea8c11f02..9fd1aed124f 100644
--- a/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.h
+++ b/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;