diff --git a/Classification/examples/Classification/gis_tutorial_example.cpp b/Classification/examples/Classification/gis_tutorial_example.cpp
index e02d128a307..ba3c15d3825 100644
--- a/Classification/examples/Classification/gis_tutorial_example.cpp
+++ b/Classification/examples/Classification/gis_tutorial_example.cpp
@@ -473,8 +473,7 @@ int main (int argc, char** argv)
   // Fill all holes except the bigest (which is the outer hull of the mesh)
   for (Mesh::Halfedge_index hi : holes)
     if (hi != outer_hull)
-      CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
-        (dtm_mesh, hi, CGAL::Emptyset_iterator(), CGAL::Emptyset_iterator());
+      CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole (dtm_mesh, hi);
 
   // Save DTM with holes filled
   std::ofstream dtm_filled_ofile ("dtm_filled.ply", std::ios_base::binary);
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 e04ebf64113..f574822c926 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
@@ -390,7 +390,7 @@ namespace Polygon_mesh_processing {
     Vertex_output_iterator vertex_out = choose_parameter<Emptyset_iterator>(get_parameter(np, internal_np::vertex_output_iterator));
 
     std::vector<typename boost::graph_traits<PolygonMesh>::face_descriptor> patch;
-    triangulate_hole(pmesh, border_halfedge, std::back_inserter(patch), np);
+    triangulate_hole(pmesh, border_halfedge, np.face_output_iterator(std::back_inserter(patch)));
     face_out = std::copy(patch.begin(), patch.end(), face_out);
 
     Hole_filling::Default_visitor default_visitor;
@@ -584,7 +584,7 @@ namespace Polygon_mesh_processing {
 
     std::vector<typename boost::graph_traits<PolygonMesh>::vertex_descriptor> patch;
     face_out = triangulate_and_refine_hole
-      (pmesh, border_halfedge, face_out, std::back_inserter(patch), np).first;
+      (pmesh, border_halfedge, np.face_output_iterator(face_out).vertex_output_iterator(std::back_inserter(patch))).first;
 
     CGAL_postcondition(CGAL::is_triangle_mesh(pmesh));
 
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp
index af6459bc13f..a11648436a5 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp
@@ -112,7 +112,7 @@ void test_triangulate_hole_weight(const std::string file_name, std::size_t nb_re
   for(typename std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it) {
     std::vector<Facet_handle> patch;
     CGAL::Polygon_mesh_processing::triangulate_hole(
-      poly, *it, back_inserter(patch),CGAL::parameters::use_delaunay_triangulation(true));
+      poly, *it, CGAL::parameters::use_delaunay_triangulation(true).face_output_iterator(back_inserter(patch)));
     if(patch.empty()) { continue; }
   }
 
@@ -161,15 +161,15 @@ void test_triangulate_hole_should_be_no_output(const std::string file_name) {
 
   for(typename std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it) {
     std::vector<Facet_handle> patch;
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch),
-      CGAL::parameters::use_delaunay_triangulation(false));
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it,
+      CGAL::parameters::use_delaunay_triangulation(false).face_output_iterator(back_inserter(patch)));
     if(!patch.empty()) {
       std::cerr << "  Error: patch should be empty" << std::endl;
       assert(false);
     }
 
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch),
-      CGAL::parameters::use_delaunay_triangulation(true));
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it,
+      CGAL::parameters::use_delaunay_triangulation(true).face_output_iterator(back_inserter(patch)));
     if(!patch.empty()) {
       std::cerr << "  Error: patch should be empty" << std::endl;
       assert(false);
@@ -258,11 +258,11 @@ void test_ouput_iterators_triangulate_hole(const std::string file_name) {
   typename std::vector<Halfedge_handle>::iterator it_2 = border_reps_2.begin();
   for(typename std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it, ++it_2) {
     std::vector<Facet_handle> patch;
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch));
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, CGAL::parameters::face_output_iterator(back_inserter(patch)));
 
     std::vector<Facet_handle> patch_2 = patch;
     Facet_handle* output_it =
-      CGAL::Polygon_mesh_processing::triangulate_hole(poly_2, *it_2, &*patch_2.begin());
+      CGAL::Polygon_mesh_processing::triangulate_hole(poly_2, *it_2, CGAL::parameters::face_output_iterator(& *patch_2.begin()));
 
     if(patch.size() != (std::size_t)(output_it - &*patch_2.begin())) {
       std::cerr << "  Error: returned facet output iterator is not valid!" << std::endl;
@@ -291,8 +291,8 @@ void test_ouput_iterators_triangulate_and_refine_hole(const std::string file_nam
   for(typename std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it, ++it_2) {
     std::vector<Facet_handle> patch_facets;
     std::vector<Vertex_handle> patch_vertices;
-    CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly,
-      *it, back_inserter(patch_facets), back_inserter(patch_vertices));
+    CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly, *it,
+      CGAL::parameters::face_output_iterator(back_inserter(patch_facets)).vertex_output_iterator(back_inserter(patch_vertices)));
     // create enough space to hold outputs
     std::vector<Facet_handle> patch_facets_2 = patch_facets;
     std::vector<Vertex_handle> patch_vertices_2 = patch_vertices;
@@ -300,7 +300,7 @@ void test_ouput_iterators_triangulate_and_refine_hole(const std::string file_nam
 
     std::pair<Facet_handle*, Vertex_handle*> output_its =
       CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly_2,
-        *it_2, &*patch_facets_2.begin(), &*patch_vertices_2.begin());
+        *it_2, CGAL::parameters::face_output_iterator(& *patch_facets_2.begin()).vertex_output_iterator(& *patch_vertices_2.begin()));
 
     if(patch_facets.size() != (std::size_t) (output_its.first - &*patch_facets_2.begin())) {
       std::cout << "  Error: returned facet output iterator is not valid!" << std::endl;
@@ -337,22 +337,29 @@ void test_triangulate_refine_and_fair_hole_compile() {
   // use all param
   read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
   CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
-  (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices),
+  (poly, border_reps[0],
   CGAL::parameters::
+    face_output_iterator(back_inserter(patch_facets)).
+    vertex_output_iterator(back_inserter(patch_vertices)).
     weight_calculator(CGAL::Weights::Uniform_weight<Polyhedron>()).
     sparse_linear_solver(Default_solver()));
 
   // default solver
   read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
   CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
-    (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices),
+    (poly, border_reps[0],
     CGAL::parameters::
+      face_output_iterator(back_inserter(patch_facets)).
+      vertex_output_iterator(back_inserter(patch_vertices)).
       weight_calculator(CGAL::Weights::Uniform_weight<Polyhedron>()));
 
   // default solver and weight
   read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
   CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
-    (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices));
+    (poly, border_reps[0],
+        CGAL::parameters::
+        face_output_iterator(back_inserter(patch_facets)).
+        vertex_output_iterator(back_inserter(patch_vertices)));
 }
 
 template <class Polyhedron>
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp
index 9933d3ab217..146eed2fd87 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp
@@ -138,8 +138,10 @@ void test_triangulate_hole(const std::string file_name, bool use_cdt) {
 
   for(std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it) {
     std::vector<Facet_handle> patch;
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch),
-                                                    CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it,
+                                                    CGAL::parameters::
+                                                      face_output_iterator(std::back_inserter(patch)).
+                                                      use_2d_constrained_delaunay_triangulation(use_cdt));
     if(patch.empty()) {
       std::cerr << "  Error: empty patch created." << std::endl;
       assert(false);
@@ -163,16 +165,17 @@ void test_triangulate_hole_should_be_no_output(const std::string file_name, bool
 
   for(std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it) {
     std::vector<Facet_handle> patch;
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch),
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it,
       CGAL::parameters::use_delaunay_triangulation(false)
+        .face_output_iterator(back_inserter(patch))
         .use_2d_constrained_delaunay_triangulation(use_cdt));
     if(!patch.empty()) {
       std::cerr << "  Error: patch should be empty" << std::endl;
       assert(false);
     }
 
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch),
-      CGAL::parameters::use_delaunay_triangulation(true));
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it,
+      CGAL::parameters::use_delaunay_triangulation(true).face_output_iterator(back_inserter(patch)));
     if(!patch.empty()) {
       std::cerr << "  Error: patch should be empty" << std::endl;
       assert(false);
@@ -193,8 +196,10 @@ void test_triangulate_and_refine_hole(const std::string file_name, bool use_cdt)
     std::vector<Facet_handle> patch_facets;
     std::vector<Vertex_handle> patch_vertices;
     CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly, *it,
-      back_inserter(patch_facets), back_inserter(patch_vertices),
-                                                               CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
+      CGAL::parameters::
+        face_output_iterator(std::back_inserter(patch_facets)).
+        vertex_output_iterator(std::back_inserter(patch_vertices)).
+        use_2d_constrained_delaunay_triangulation(use_cdt));
 
     if(patch_facets.empty()) {
       std::cerr << "  Error: empty patch created." << std::endl;
@@ -220,9 +225,11 @@ void test_triangulate_refine_and_fair_hole(const std::string file_name, bool use
   for(std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it) {
     std::vector<Facet_handle> patch_facets;
     std::vector<Vertex_handle> patch_vertices;
-    CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole(poly,
-      *it, back_inserter(patch_facets), back_inserter(patch_vertices),
-      CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
+    CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole(poly, *it,
+      CGAL::parameters::
+        face_output_iterator(std::back_inserter(patch_facets)).
+        vertex_output_iterator(std::back_inserter(patch_vertices)).
+        use_2d_constrained_delaunay_triangulation(use_cdt));
 
     if(patch_facets.empty()) {
       std::cerr << "  Error: empty patch created." << std::endl;
@@ -251,12 +258,14 @@ void test_ouput_iterators_triangulate_hole(const std::string file_name, bool use
   std::vector<Halfedge_handle>::iterator it_2 = border_reps_2.begin();
   for(std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it, ++it_2) {
     std::vector<Facet_handle> patch;
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it, back_inserter(patch),
-              CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, *it,
+      CGAL::parameters::
+        face_output_iterator(std::back_inserter(patch)).
+        use_2d_constrained_delaunay_triangulation(use_cdt));
 
     std::vector<Facet_handle> patch_2 = patch;
     Facet_handle* output_it =
-      CGAL::Polygon_mesh_processing::triangulate_hole(poly_2, *it_2, &*patch_2.begin());
+      CGAL::Polygon_mesh_processing::triangulate_hole(poly_2, *it_2, CGAL::parameters::face_output_iterator(& *patch_2.begin()));
 
     if(patch.size() != (std::size_t)(output_it - &*patch_2.begin())) {
       std::cerr << "  Error: returned facet output iterator is not valid!" << std::endl;
@@ -282,18 +291,22 @@ void test_ouput_iterators_triangulate_and_refine_hole(const std::string file_nam
   for(std::vector<Halfedge_handle>::iterator it = border_reps.begin(); it != border_reps.end(); ++it, ++it_2) {
     std::vector<Facet_handle> patch_facets;
     std::vector<Vertex_handle> patch_vertices;
-    CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly,
-      *it, back_inserter(patch_facets), back_inserter(patch_vertices),
-      CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
+    CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly, *it,
+      CGAL::parameters::
+        face_output_iterator(std::back_inserter(patch_facets)).
+        vertex_output_iterator(std::back_inserter(patch_vertices)).
+        use_2d_constrained_delaunay_triangulation(use_cdt));
     // create enough space to hold outputs
     std::vector<Facet_handle> patch_facets_2 = patch_facets;
     std::vector<Vertex_handle> patch_vertices_2 = patch_vertices;
     if(patch_vertices_2.empty()) { patch_vertices_2.push_back(Vertex_handle()); } //just allocate space for dereferencing
 
     std::pair<Facet_handle*, Vertex_handle*> output_its =
-      CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly_2,
-        *it_2, &*patch_facets_2.begin(), &*patch_vertices_2.begin(),
-        CGAL::parameters::use_2d_constrained_delaunay_triangulation(use_cdt));
+      CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly_2, *it_2,
+        CGAL::parameters::
+          face_output_iterator(&*patch_facets_2.begin()).
+          vertex_output_iterator(&*patch_vertices_2.begin()).
+          use_2d_constrained_delaunay_triangulation(use_cdt));
 
     if(patch_facets.size() != (std::size_t) (output_its.first - &*patch_facets_2.begin())) {
       std::cout << "  Error: returned facet output iterator is not valid!" << std::endl;
@@ -327,23 +340,30 @@ void test_triangulate_refine_and_fair_hole_compile() {
   // use all param
   read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
   CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
-  (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices),
-  CGAL::parameters::
-    weight_calculator(CGAL::Weights::Uniform_weight<Polyhedron>()).
-    sparse_linear_solver(Default_solver()).
-    use_2d_constrained_delaunay_triangulation(false));
+  (poly, border_reps[0],
+   CGAL::parameters::
+     face_output_iterator(back_inserter(patch_facets)).
+     vertex_output_iterator(back_inserter(patch_vertices)).
+     weight_calculator(CGAL::Weights::Uniform_weight<Polyhedron>()).
+     sparse_linear_solver(Default_solver()).
+     use_2d_constrained_delaunay_triangulation(false));
 
   // default solver
   read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
   CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
-    (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices),
-    CGAL::parameters::
-      weight_calculator(CGAL::Weights::Uniform_weight<Polyhedron>()));
+    (poly, border_reps[0],
+     CGAL::parameters::
+       face_output_iterator(back_inserter(patch_facets)).
+       vertex_output_iterator(back_inserter(patch_vertices)).
+       weight_calculator(CGAL::Weights::Uniform_weight<Polyhedron>()));
 
   // default solver and weight
   read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
   CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
-    (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices));
+    (poly, border_reps[0],
+     CGAL::parameters::
+       face_output_iterator(back_inserter(patch_facets)).
+       vertex_output_iterator(back_inserter(patch_vertices)));
 }
 
 void generate_elephant_with_hole()
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_with_cdt_2_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_with_cdt_2_test.cpp
index 8ab2e2afd71..8f1ede58b93 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_with_cdt_2_test.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_with_cdt_2_test.cpp
@@ -86,9 +86,9 @@ void test_triangulate_hole_with_cdt_2(
     CGAL::Polygon_mesh_processing::triangulate_hole(
       pmesh,
       h,
-      std::back_inserter(patch_faces),
-      CGAL::parameters::vertex_point_map(
-        get(CGAL::vertex_point, pmesh)).
+      CGAL::parameters::
+        face_output_iterator(std::back_inserter(patch_faces)).
+        vertex_point_map(get(CGAL::vertex_point, pmesh)).
         use_2d_constrained_delaunay_triangulation(true).
         geom_traits(GeomTraits()));
 
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Hole_filling_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Hole_filling_plugin.cpp
index fbc5fe8af1b..da3c06dbf72 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Hole_filling_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Hole_filling_plugin.cpp
@@ -701,9 +701,10 @@ bool Polyhedron_demo_hole_filling_plugin::fill
   CGAL::Timer timer; timer.start();
   std::vector<fg_face_descriptor> patch;
   if(action_index == 0) {
-    CGAL::Polygon_mesh_processing::triangulate_hole(poly,
-             it, std::back_inserter(patch),
-             CGAL::parameters::use_delaunay_triangulation(use_DT));
+    CGAL::Polygon_mesh_processing::triangulate_hole(poly, it,
+             CGAL::parameters::
+             face_output_iterator(std::back_inserter(patch)).
+             use_delaunay_triangulation(use_DT));
   }
   else if(action_index == 1) {
     CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(poly, it,