diff --git a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp
index 09383c8f7bc..f92c252014a 100644
--- a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp
+++ b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp
@@ -45,7 +45,7 @@ int main(int argc, char** argv)
   SMS::Count_ratio_stop_predicate<Surface_mesh> stop(ratio);
 
   // Garland&Heckbert simplification policies
-  typedef typename SMS::GarlandHeckbert_policies<Surface_mesh, Kernel>          GH_policies;
+  typedef typename SMS::GarlandHeckbert_plane_policies<Surface_mesh, Kernel>    GH_policies;
   typedef typename GH_policies::Get_cost                                        GH_cost;
   typedef typename GH_policies::Get_placement                                   GH_placement;
   typedef SMS::Bounded_normal_change_placement<GH_placement>                    Bounded_GH_placement;
diff --git a/Surface_mesh_simplification/examples/Surface_mesh_simplification/gh_statistics.cpp b/Surface_mesh_simplification/examples/Surface_mesh_simplification/gh_statistics.cpp
index e95c76c1f31..6e11e5dce20 100644
--- a/Surface_mesh_simplification/examples/Surface_mesh_simplification/gh_statistics.cpp
+++ b/Surface_mesh_simplification/examples/Surface_mesh_simplification/gh_statistics.cpp
@@ -4,28 +4,27 @@
 #include <CGAL/Surface_mesh_simplification/edge_collapse.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_ratio_stop_predicate.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_placement.h>
-#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_policies.h>
-#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_tri_policies.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_plane_policies.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_triangle_policies.h>
-#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h>
-
-#include <iostream>
-#include <fstream>
-
-#include <chrono>
-#include <vector>
-#include <array>
-#include <functional>
-
-#include <boost/format.hpp> 
-#include <boost/filesystem.hpp>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_plane_policies.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_triangle_policies.h>
 
 #include <CGAL/Polygon_mesh_processing/measure.h>
 #include <CGAL/Polygon_mesh_processing/distance.h>
 
+#include <boost/format.hpp>
+#include <boost/filesystem.hpp>
+
+#include <array>
+#include <chrono>
+#include <iostream>
+#include <fstream>
+#include <functional>
+#include <vector>
+
 /**
  * this is a largely undocumented file that was used to gather most statistics
- * for my GSoC 2021 project on mesh decimation, maybe it is useful in the future  
+ * for my GSoC 2021 project on mesh decimation, maybe it is useful in the future
  */
 
 typedef CGAL::Simple_cartesian<double>                          Kernel;
@@ -40,10 +39,10 @@ namespace SMS = CGAL::Surface_mesh_simplification;
 namespace PMP = CGAL::Polygon_mesh_processing;
 namespace fs = boost::filesystem;
 
-typedef SMS::GarlandHeckbert_policies<Surface_mesh, Kernel> Classic_plane;
-typedef SMS::GarlandHeckbert_probabilistic_policies<Surface_mesh, Kernel> Prob_plane;
+typedef SMS::GarlandHeckbert_plane_policies<Surface_mesh, Kernel> Classic_plane;
 typedef SMS::GarlandHeckbert_triangle_policies<Surface_mesh, Kernel> Classic_tri;
-typedef SMS::GarlandHeckbert_probabilistic_tri_policies<Surface_mesh, Kernel> Prob_tri;
+typedef SMS::GarlandHeckbert_probabilistic_plane_policies<Surface_mesh, Kernel> Prob_plane;
+typedef SMS::GarlandHeckbert_probabilistic_triangle_policies<Surface_mesh, Kernel> Prob_tri;
 
 // settings for benchmarking - throw away the first n_burns results and keep the n_samples
 // samples
@@ -55,96 +54,90 @@ constexpr std::size_t prob_plane_index = 1;
 constexpr std::size_t classic_tri_index = 2;
 constexpr std::size_t prob_tri_index = 3;
 
-template<typename Policy>
-Surface_mesh edge_collapse(Surface_mesh& mesh, double ratio = 0.2)
+template <typename Policy>
+Surface_mesh edge_collapse(Surface_mesh& mesh,
+                           const double ratio = 0.2)
 {
-  typedef typename Policy::Get_cost Cost; 
+  typedef typename Policy::Get_cost Cost;
   typedef typename Policy::Get_placement Placement;
-  
+
   typedef SMS::Bounded_normal_change_placement<Placement> Bounded_placement;
-  
-  const Policy p {mesh, 100};
-  
+
+  const Policy p { mesh, 100 };
+
   const Cost& cost = p.get_cost();
   const Placement& unbounded_placement = p.get_placement();
 
   Bounded_placement placement(unbounded_placement);
 
   SMS::Count_ratio_stop_predicate<Surface_mesh> stop(ratio);
-  
+
   // collapse edges ignoring result code
-  SMS::edge_collapse(mesh, stop,
-      CGAL::parameters::get_cost(cost).get_placement(placement));
-  
+  SMS::edge_collapse(mesh, stop, CGAL::parameters::get_cost(cost).get_placement(placement));
+
   return mesh;
 }
 
-template<typename Policy>
-void time_all_vector(const std::vector<Surface_mesh>& meshes, const fs::path& output_file)
+template <typename Policy>
+void time_all_vector(const fs::path& output_file,
+                     const std::vector<Surface_mesh>& meshes)
 {
-  namespace time = std::chrono;
+  fs::ofstream out_stream { output_file };
 
-  fs::ofstream out_stream {output_file};
-  
   // always decimate meshes in this vector to avoid timing the
   // copying of the meshes
   std::vector<Surface_mesh> temp_meshes = meshes;
-  
-  for(int i = 0; i < n_burns; ++i)
+
+  for(int i=0; i<n_burns; ++i)
   {
     for(Surface_mesh& mesh : temp_meshes)
-    {
       edge_collapse<Policy>(mesh);
-    }
   }
+
   temp_meshes = meshes;
 
-  // measure each run 
-  for(int i = 0; i < n_samples; ++i)
+  // measure each run
+  for(int i=0; i<n_samples; ++i)
   {
     // measure time taken by the edge_collapse function over each mesh individually
-    time::time_point start_time = time::steady_clock::now();
+    std::chrono::time_point<std::chrono::steady_clock> start_time = std::chrono::steady_clock::now();
 
     for(Surface_mesh& mesh : temp_meshes)
-    {
       edge_collapse<Policy>(mesh);
-    }
 
-    time::time_point end_time = time::steady_clock::now();
+    std::chrono::time_point<std::chrono::steady_clock> end_time = std::chrono::steady_clock::now();
 
     // get elapsed time in nanoseconds
-    unsigned long elapsed_ns = time::duration_cast<time::nanoseconds>
-      (end_time - start_time).count();
+    unsigned long elapsed_ns = std::chrono::duration_cast<std::chrono::nanoseconds>(end_time - start_time).count();
 
     // add the elapsed time as data point
-    out_stream << std::to_string(elapsed_ns) << '\n';   
+    out_stream << std::to_string(elapsed_ns) << '\n';
 
     temp_meshes = meshes;
   }
 }
 
-bool read_from_file(Surface_mesh& mesh, const fs::path& p)
+bool read_from_file(const fs::path& p,
+                    Surface_mesh& mesh)
 {
-  fs::ifstream in {p};
+  fs::ifstream in { p };
 
   // make sure the mesh is empty and ready for the next set of data
   mesh.clear();
 
-  if (!in || !(in >> mesh))
-  {
+  if(!in || !(in >> mesh))
     return false;
-  }
 
   return true;
 }
 
 std::vector<std::pair<Surface_mesh, std::string>> get_all_meshes(const fs::path& dir)
 {
-  // vector for storing results 
+  // vector for storing results
   std::vector<std::pair<Surface_mesh, std::string>> meshes { };
 
   // read all meshes in the given directory into the vector
-  if (fs::is_directory(dir))
+  if(fs::is_directory(dir))
   {
     fs::directory_iterator end_iter;
     Surface_mesh mesh;
@@ -153,25 +146,24 @@ std::vector<std::pair<Surface_mesh, std::string>> get_all_meshes(const fs::path&
     for(fs::directory_iterator dir_iter(dir); dir_iter != end_iter; ++dir_iter)
     {
       // look for .off files
-      if (fs::is_regular_file(dir_iter->status())
-          && dir_iter->path().extension() == ".off")
+      if(fs::is_regular_file(dir_iter->status()) &&
+         dir_iter->path().extension() == ".off")
       {
         const fs::path& path = dir_iter->path();
 
-        // try to read the file into the surface mesh, upon failure, we continue looping 
+        // try to read the file into the surface mesh, upon failure, we continue looping
         // iterating through the directory
-        if (!read_from_file(mesh, path))
+        if(!read_from_file(path, mesh))
         {
-          std::cerr << "Failed to read input mesh " << path 
-            << " in directory " << dir << "." << std::endl;
+          std::cerr << "Failed to read input mesh " << path << " in directory " << dir << "." << std::endl;
         }
-        else if (!CGAL::is_triangle_mesh(mesh))
+        else if(!CGAL::is_triangle_mesh(mesh))
         {
           std::cerr << "Input geometry is not triangulated." << std::endl;
         }
         else
         {
-          // we can move this mesh since we don't need it anymore - the next one will be read into 
+          // we can move this mesh since we don't need it anymore - the next one will be read into
           // the same variable
           meshes.emplace_back(std::move(mesh), path.filename().string());
         }
@@ -182,43 +174,39 @@ std::vector<std::pair<Surface_mesh, std::string>> get_all_meshes(const fs::path&
   return meshes;
 }
 
-  template<typename Policy>
+template <typename Policy>
 void time_policy(const fs::path& dir, const fs::path& output_file)
 {
   auto vec = get_all_meshes(dir);
-  std::vector<Surface_mesh> meshes { };
+  std::vector<Surface_mesh> meshes;
 
-  for(const auto& p : vec) {
+  for(const auto& p : vec)
     std::cout << p.second << '\n';
-  }
 
   // keep the first values from the pair (we don't need the names)(we don't need the names)
   std::transform(vec.begin(), vec.end(), std::back_inserter(meshes),
-      [] (const std::pair<Surface_mesh, std::string>& p) { return p.first; }
-      );
+                 [] (const std::pair<Surface_mesh, std::string>& p) { return p.first; } );
 
-  return time_all_vector<Policy>(meshes, output_file);
+  return time_all_vector<Policy>(output_file, meshes);
 }
 
-  template<typename Policy>
+template <typename Policy>
 double hausdorff_error(const Surface_mesh& mesh, double ratio = 0.2)
 {
-  // an arbitrarily chosen small value 
+  // an arbitrarily chosen small value
   constexpr double error_bound = 0.00001;
 
   // make a copy of the mesh so we can compare later
   Surface_mesh temp = mesh;
 
   edge_collapse<Policy>(temp, ratio);
-  return PMP::bounded_error_symmetric_Hausdorff_distance<CGAL::Sequential_tag>
-    (mesh, temp, error_bound);
+  return PMP::bounded_error_symmetric_Hausdorff_distance<CGAL::Parallel_if_available_tag>(mesh, temp, error_bound);
 }
 
-// calculate approximate Hausdorff errors for all different policies at 
-// the same decimation ratio
+// calculate approximate Hausdorff errors for all different policies at the same decimation ratio
 std::array<FT, 4> hausdorff_errors(const Surface_mesh& mesh, double ratio)
 {
-  std::array<FT, 4> ret { {0, 0, 0, 0} }; 
+  std::array<FT, 4> ret { {0, 0, 0, 0} };
 
   ret[classic_plane_index] = hausdorff_error<Classic_plane>(mesh, ratio);
   ret[prob_plane_index] = hausdorff_error<Prob_plane>(mesh, ratio);
@@ -228,10 +216,10 @@ std::array<FT, 4> hausdorff_errors(const Surface_mesh& mesh, double ratio)
   return ret;
 }
 
-  template<typename InputIt>
-void hausdorff_errors_mesh(const Surface_mesh& mesh, const fs::path& out, InputIt begin, InputIt end)
+template <typename InputIt>
+void hausdorff_errors_mesh(const fs::path& out, const Surface_mesh& mesh, InputIt begin, InputIt end)
 {
-  fs::ofstream out_stream {out};
+  fs::ofstream out_stream { out };
 
   for(InputIt it = begin; it != end; ++it)
   {
@@ -240,34 +228,28 @@ void hausdorff_errors_mesh(const Surface_mesh& mesh, const fs::path& out, InputI
     out_stream << "ratio: " << *it << '\n';
 
     for(double e : errs)
-    {
       out_stream << std::to_string(e) << '\n';
-    }
   }
 }
 
-  template<typename InputIt>
-void hausdorff_errors_mesh(const fs::path& in, const fs::path& out, InputIt begin, InputIt end)
+template <typename InputIt>
+void hausdorff_errors_mesh(const fs::path& out, const fs::path& in, InputIt begin, InputIt end)
 {
   fs::ifstream is {in};
   Surface_mesh mesh;
 
-  if (!read_from_file(mesh, in))
-  {
-    std::cerr << "Failed to read input mesh " << in << '\n'; 
-  }
+  if(!read_from_file(in, mesh))
+    std::cerr << "Failed to read input mesh " << in << '\n';
   else
-  {
-    hausdorff_errors_mesh(mesh, out, begin, end);
-  }
-} 
+    hausdorff_errors_mesh(out, mesh, begin, end);
+}
 
 void hausdorff_errors_dir(const fs::path& dir, const fs::path& out, double ratio = 0.2)
 {
   // get a vector of pairs of a mesh and its name by loading all .off files in a directory
-  const auto meshes = get_all_meshes(dir); 
+  const auto meshes = get_all_meshes(dir);
 
-  fs::ofstream out_stream {out};
+  fs::ofstream out_stream { out };
 
   for(auto& p : meshes)
   {
@@ -276,31 +258,27 @@ void hausdorff_errors_dir(const fs::path& dir, const fs::path& out, double ratio
 
     out_stream << p.second << '\n';
     for(double err : errors)
-    {
       out_stream << std::to_string(err) << '\n';
-    }
   }
 }
 
-void write_aspect_ratios(fs::ofstream& out, const Surface_mesh& mesh) 
+void write_aspect_ratios(fs::ofstream& out, const Surface_mesh& mesh)
 {
   for(auto face : faces(mesh))
-  {
-    out << std::to_string(PMP::face_aspect_ratio(face, mesh)) << '\n'; 
-  }
+    out << std::to_string(PMP::face_aspect_ratio(face, mesh)) << '\n';
 }
 
 void gather_face_aspect_ratio(const fs::path& file, const fs::path& out)
 {
-  Surface_mesh cp { }; 
-  Surface_mesh ct { }; 
-  Surface_mesh pp { }; 
-  Surface_mesh pt { }; 
+  Surface_mesh cp { };
+  Surface_mesh ct { };
+  Surface_mesh pp { };
+  Surface_mesh pt { };
 
-  read_from_file(cp, file);
-  read_from_file(ct, file);
-  read_from_file(pp, file);
-  read_from_file(pt, file);
+  read_from_file(file, cp);
+  read_from_file(file, ct);
+  read_from_file(file, pp);
+  read_from_file(file, pt);
 
   edge_collapse<Classic_plane>(cp);
   edge_collapse<Prob_plane>(pp);
@@ -321,40 +299,32 @@ void gather_face_aspect_ratio(const fs::path& file, const fs::path& out)
 
 enum class Mode { time, hausdorff_ratio, hausdorff_mesh, run, apect_ratio };
 
-void time_policy_string(const fs::path& in, const fs::path& out, const std::string& policy) {
-  if (policy == "classic_plane")
-  {
+void time_policy_string(const fs::path& in, const fs::path& out, const std::string& policy)
+{
+  if(policy == "classic_plane")
     time_policy<Classic_plane>(in, out);
-  }
-  else if (policy == "prob_plane")
-  {
+  else if(policy == "prob_plane")
     time_policy<Prob_plane>(in, out);
-  }
-  else if (policy == "classic_tri")
-  {
+  else if(policy == "classic_tri")
     time_policy<Classic_tri>(in, out);
-  }
-  else if (policy == "prob_tri")
-  {
+  else if(policy == "prob_tri")
     time_policy<Prob_tri>(in, out);
-  }
 }
 
 int main(int argc, char* argv[])
 {
-
   // default is time
   Mode m = Mode::time;
 
-  if (argc > 1)
+  if(argc > 1)
   {
     std::string input { argv[1] };
 
-    if (input == "time")
+    if(input == "time")
     {
       m = Mode::time;
 
-      if (argc != 5)
+      if(argc != 5)
       {
         std::cerr << "Expected 3 arugments to time: [input_dir] [output_file] [policy].\n";
         return EXIT_FAILURE;
@@ -362,36 +332,36 @@ int main(int argc, char* argv[])
 
       time_policy_string(argv[2], argv[3], argv[4]);
     }
-    else if (input == "ratio")
+    else if(input == "ratio")
     {
       m = Mode::hausdorff_ratio;
     }
-    else if (input == "mesh")
+    else if(input == "mesh")
     {
       m = Mode::hausdorff_mesh;
     }
-    else if (input == "run")
+    else if(input == "run")
     {
       m = Mode::run;
-      if (argc != 3 && argc != 4)
+      if(argc != 3 && argc != 4)
       {
         std::cerr << "Expected 1 (2) argument(s) to run: [input_file] [ratio].\n";
         return EXIT_FAILURE;
       }
       // default init empty mesh
-      Surface_mesh m { };
+      Surface_mesh sm;
 
       FT ratio = 0.2;
-      if (argc == 4) {
+      if(argc == 4) {
         ratio = std::stod(argv[3]);
       }
 
-      read_from_file(m, argv[2]);
+      read_from_file(argv[2], sm);
 
-      Surface_mesh cp = m;
-      Surface_mesh pp = m;
-      Surface_mesh ct = m;
-      Surface_mesh pt = m;
+      Surface_mesh cp = sm;
+      Surface_mesh pp = sm;
+      Surface_mesh ct = sm;
+      Surface_mesh pt = sm;
 
       edge_collapse<Classic_plane>(cp, ratio);
       edge_collapse<Prob_plane>(pp, ratio);
@@ -408,27 +378,27 @@ int main(int argc, char* argv[])
       classic_tri_stream << ct;
       prob_tri_stream << pt;
     }
-    else if (input == "aspect")
+    else if(input == "aspect")
     {
       gather_face_aspect_ratio(argv[2], argv[3]);
-    } 
-    else 
+    }
+    else
     {
-      std::cerr << "Didn't recognise command line argument " << input << 
-        "(options are time, ratio or mesh)\n"; 
+      std::cerr << "Didn't recognise command line argument " << input <<
+        "(options are time, ratio or mesh)\n";
 
       return EXIT_FAILURE;
     }
   }
 
-  if (m == Mode::hausdorff_mesh)
+  if(m == Mode::hausdorff_mesh)
   {
-    if (argc != 3)
+    if(argc != 3)
     {
       std::cerr << "Expected 1 argument to mesh: [input_file].\n";
       return EXIT_FAILURE;
     }
-    else 
+    else
     {
       std::array<double, 10> range {{ 0.7, 0.6, 0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15 }};
       hausdorff_errors_mesh(argv[2], "mesh_hausdorff_err", range.begin(), range.end());