diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/snap_polygon_soup.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/snap_polygon_soup.cpp
index d0265ee9fab..e23fcdbaa9f 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/snap_polygon_soup.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/snap_polygon_soup.cpp
@@ -1,12 +1,13 @@
-#define PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+// #define PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
 // #define CGAL_PMP_AUTOREFINE_USE_DEFAULT_VERBOSE
 
-#include <CGAL/Cartesian.h>
+#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
 #include <CGAL/Polygon_mesh_processing/autorefinement.h>
 #include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
 #include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
 #include <CGAL/IO/polygon_soup_io.h>
 #include <CGAL/Real_timer.h>
@@ -15,9 +16,9 @@
 
 #include <iostream>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-typedef CGAL::Cartesian<double> Cartesian;
-typedef Kernel::Point_3                                     Point;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK;
+typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK;
+typedef CGAL::Simple_cartesian<double> Cartesian;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -26,31 +27,76 @@ int main(int argc, char** argv)
   const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
                                          : std::string(argv[1]);
 
+  // const std::string out_file = argc <= 2 ? "rounded_soup.off"
+  //                                        : std::string(argv[2]);
+  const std::string out_file = "rounded_soup.off";
+
   const int grid_size = argc <= 2 ? 23
                                   : std::stoi(std::string(argv[2]));
 
-  std::vector<Point> input_points;
-  std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
-  if (!CGAL::IO::read_polygon_soup(filename, input_points, input_triangles))
+  const bool epeck = argc <= 3 ? false
+                               : std::string(argv[3])=="EPECK";
+
+  if(epeck)
   {
-    std::cerr << "Cannot read " << filename << "\n";
-    return 1;
+    std::vector<EPECK::Point_3> input_points;
+    std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
+
+    std::cout << "Snap rounding apply on " << filename << " with EPECK\n";
+    if (!CGAL::IO::read_polygon_soup(filename, input_points, input_triangles))
+    {
+      std::cerr << "Cannot read " << filename << "\n";
+      return 1;
+    }
+    std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << std::endl;
+    std::cout << "Is 2-manifold: " << PMP::is_polygon_soup_a_polygon_mesh(input_triangles) << std::endl;
+
+    PMP::repair_polygon_soup(input_points, input_triangles);
+    PMP::triangulate_polygons(input_points, input_triangles);
+
+    CGAL::Real_timer t;
+    t.start();
+    PMP::autorefine_triangle_soup(input_points, input_triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size));
+    t.stop();
+    std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << " in " << t.time() << " sec." << std::endl;
+    std::cout << "Does self-intersect: " << PMP::does_triangle_soup_self_intersect(input_points, input_triangles) << std::endl;
+    std::cout << "Is 2-manifold: " << PMP::orient_polygon_soup(input_points, input_triangles) << "\n\n"  << std::endl;
+
+    std::vector<Cartesian::Point_3> output_points;
+    for(auto &p: input_points)
+      output_points.emplace_back(CGAL::to_double(p.x()),CGAL::to_double(p.y()),CGAL::to_double(p.z()));
+
+    CGAL::IO::write_polygon_soup(out_file, output_points, input_triangles, CGAL::parameters::stream_precision(17));
   }
-  PMP::repair_polygon_soup(input_points, input_triangles);
-  PMP::triangulate_polygons(input_points, input_triangles);
+  else
+  {
+    std::vector<EPICK::Point_3> input_points;
+    std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
+    std::cout << "Snap rounding on " << filename << " with EPICK\n";
+    if (!CGAL::IO::read_polygon_soup(filename, input_points, input_triangles))
+    {
+      std::cerr << "Cannot read " << filename << "\n";
+      return 1;
+    }
+    std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << std::endl;
+    std::cout << "Is 2-manifold: " << PMP::is_polygon_soup_a_polygon_mesh(input_triangles) << std::endl;
 
-  CGAL::Real_timer t;
-  t.start();
-  PMP::autorefine_triangle_soup(input_points, input_triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size));
-  t.stop();
-  std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << " in " << t.time() << " sec." << std::endl;
+    std::vector<std::pair<std::size_t, std::size_t>> pairs_of_intersecting_triangles;
+    PMP::triangle_soup_self_intersections(input_points, input_triangles, std::back_inserter(pairs_of_intersecting_triangles));
+    std::cout << "Nb of pairs of intersecting triangles: " << pairs_of_intersecting_triangles.size() << std::endl;
 
-  std::vector<Cartesian::Point_3> output_points;
-  std::cout << "Does self-intersect: " << PMP::does_triangle_soup_self_intersect(input_points, input_triangles) << std::endl;
-  for(auto &p: input_points)
-    output_points.emplace_back(CGAL::to_double(p.x()),CGAL::to_double(p.y()),CGAL::to_double(p.z()));
-
-  CGAL::IO::write_polygon_soup("rounded_soup.off", output_points, input_triangles, CGAL::parameters::stream_precision(17));
+    PMP::repair_polygon_soup(input_points, input_triangles);
+    PMP::triangulate_polygons(input_points, input_triangles);
 
+    CGAL::Real_timer t;
+    t.start();
+    PMP::autorefine_triangle_soup(input_points, input_triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size));
+    t.stop();
+    std::cout << "\nOutput:" << std::endl;
+    std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << " in " << t.time() << " sec." << std::endl;
+    std::cout << "Does self-intersect: " << PMP::does_triangle_soup_self_intersect(input_points, input_triangles) << std::endl;
+    std::cout << "Is 2-manifold: " << PMP::orient_polygon_soup(input_points, input_triangles) << "\n\n"  << std::endl;
+    CGAL::IO::write_polygon_soup(out_file, input_points, input_triangles, CGAL::parameters::stream_precision(17));
+  }
   return 0;
 }