Merge pull request #8744 from LeoValque/PMP_triangle_soup_rounding-GF

Add do_snap parameter to PMP::autorefine_triangle_soup

Documentation/doc/biblio/geom.bib

@@ -152087,6 +152087,18 @@ keywords = {polygonal surface mesh, Surface reconstruction, kinetic framework, s
   publisher={Elsevier}
 }
 
+@unpublished{lazard:hal-04907149,
+  TITLE = {{Removing self-intersections in 3D meshes while preserving floating-point coordinates}},
+  AUTHOR = {Lazard, Sylvain and Valque, Leo},
+  URL = {https://inria.hal.science/hal-04907149},
+  NOTE = {working paper or preprint},
+  YEAR = {2025},
+  MONTH = Jan,
+  KEYWORDS = {Snap rounding ; mesh intersection ; robustness},
+  PDF = {https://inria.hal.science/hal-04907149v1/file/Snap-HAL.pdf},
+  HAL_ID = {hal-04907149},
+  HAL_VERSION = {v1},
+}
 
 @inproceedings{si2005meshing,
   title={Meshing piecewise linear complexes by constrained {Delaunay} tetrahedralizations},

Installation/CHANGES.md

@@ -39,6 +39,7 @@
 - New implementation of `CGAL::Polygon_mesh_processing::clip()` with a plane as clipper that is much faster and is now able to handle non-triangulated surface meshes.
 - New implementation of `CGAL::Polygon_mesh_processing::split()` with a plane as clipper that is much faster and is now able to handle non-triangulated surface meshes.
 - Added the function `CGAL::Polygon_mesh_processing::refine_with_plane()`, which enables users to refine a mesh with their intersection with a plane.
+- Added the parameter `apply_iterative_snap_rounding` to the function `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()`. When set to `true`, the coordinates are rounded to fit in double and may perform additional subdivisions to ensure the output is free of self-intersections.
 
 ### [Point Set Processing](https://doc.cgal.org/6.1/Manual/packages.html#PkgPointSetProcessing3)
 - Added `poisson_eliminate()` to downsample a point cloud to a target size while providing Poisson disk property, i.e., a larger minimal distance between points.

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/CMakeLists.txt

@@ -37,7 +37,21 @@ else()
 endif()
 
 create_single_source_cgal_program("fast.cpp")
+create_single_source_cgal_program("rotated_cubes_autorefinement.cpp")
+create_single_source_cgal_program("coplanar_cubes_autorefinement.cpp")
+
+create_single_source_cgal_program("Performance/performance_snap_polygon_soup.cpp")
+create_single_source_cgal_program("Robustness/robustness_snap_polygon_soup.cpp")
+create_single_source_cgal_program("Quality/quality_snap_polygon_soup.cpp")
 
 create_single_source_cgal_program("polygon_mesh_slicer.cpp")
 target_link_libraries(polygon_mesh_slicer PRIVATE CGAL::Eigen3_support)
 
+find_package(TBB QUIET)
+include(CGAL_TBB_support)
+if(TARGET CGAL::TBB_support)
+  target_link_libraries(rotated_cubes_autorefinement PRIVATE CGAL::TBB_support)
+  target_link_libraries(coplanar_cubes_autorefinement PRIVATE CGAL::TBB_support)
+else()
+  message(STATUS "NOTICE: Intel TBB was not found. Sequential code will be used.")
+endif()

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Performance/performance_snap_polygon_soup.cpp

@@ -0,0 +1,50 @@
+#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/triangulate_faces.h>
+#include <CGAL/IO/polygon_soup_io.h>
+
+#include <boost/container/small_vector.hpp>
+
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef typename Kernel::Point_3 Point_3;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+enum EXIT_CODES { VALID_OUTPUT=0,
+                  INVALID_INPUT=1,
+                  ROUNDING_FAILED=2,
+                  SIGSEGV=10,
+                  SIGSABRT=11,
+                  SIGFPE=12,
+                  TIMEOUT=13
+                  };
+
+int main(int argc, char** argv)
+{
+  if(argc<4){
+    std::cout << "Invalid argument" << std::endl;
+    return 1;
+  }
+
+  const std::string filename = std::string(argv[1]);
+  const int grid_size = std::stoi(std::string(argv[2]));
+  const bool erase_duplicate = std::stoi(argv[3])==1;
+
+  std::vector<Point_3> points;
+  std::vector<boost::container::small_vector<std::size_t, 3>> triangles;
+
+  if (!CGAL::IO::read_polygon_soup(filename, points, triangles))
+  {
+    std::cerr << "Cannot read " << filename << "\n";
+    return 1;
+  }
+
+  PMP::repair_polygon_soup(points, triangles);
+  PMP::triangulate_polygons(points, triangles);
+
+  PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(erase_duplicate).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).number_of_iterations(15));
+
+  return 0;
+}

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Performance/run_performance.sh

@@ -0,0 +1,28 @@
+#!/bin/bash
+
+set -e
+
+if [ "$#" -lt 4 ]; then
+    echo "Usage: $0 <input_file> <timeout> [component_params...]"
+    exit 1
+fi
+
+INPUT_FILE=$1
+TIMEOUT=$2
+GRID_SIZE=$3
+ERASE_ALL_DUPLICATE=$4
+
+# Use /usr/bin/time for memory usage (maximum resident set size in KB)
+TMP_LOG=$(mktemp)
+
+# Run the benchmarked command
+/usr/bin/time -f "TIME:%e\nMEM:%M" timeout "$TIMEOUT"s performance_snap_polygon_soup "$INPUT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" 2> "$TMP_LOG"
+
+# Parse time and memory
+SECONDS=$(grep "TIME" "$TMP_LOG" | cut -d':' -f2)
+MEMORY=$(grep "MEM" "$TMP_LOG" | cut -d':' -f2)
+
+rm -f "$TMP_LOG"
+
+# Output JSON
+echo "{\"seconds\": \"$SECONDS\", \"memory_peaks\": \"$MEMORY\"}"

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Quality/quality_snap_polygon_soup.cpp

@@ -0,0 +1,60 @@
+#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/triangulate_faces.h>
+#include <CGAL/Polygon_mesh_processing/distance.h>
+#include <CGAL/IO/polygon_soup_io.h>
+
+#include <CGAL/Bbox_3.h>
+#include <CGAL/boost/graph/helpers.h>
+
+#include <CGAL/Surface_mesh.h>
+
+#include <boost/container/small_vector.hpp>
+#include <CGAL/Polygon_mesh_processing/orientation.h>
+
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+typedef typename Kernel::Point_3 Point_3;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char** argv)
+{
+  if(argc<4){
+    std::cout << "Invalid argument" << std::endl;
+    return 1;
+  }
+  const std::string filename = std::string(argv[1]);
+  const int grid_size = std::stoi(std::string(argv[2]));
+  const bool erase_duplicate = std::stoi(argv[3])==1;
+
+  std::vector<Point_3> points;
+  std::vector<boost::container::small_vector<std::size_t, 3>> triangles;
+
+  CGAL::Bbox_3 bb = CGAL::bbox_3(points.begin(), points.end());
+  double diag_length=std::sqrt((bb.xmax()-bb.xmin())*(bb.xmax()-bb.xmin()) + (bb.ymax()-bb.ymin())*(bb.ymax()-bb.ymin()) + (bb.zmax()-bb.zmin())*(bb.zmax()-bb.zmin()));
+  if (!CGAL::IO::read_polygon_soup(filename, points, triangles))
+  {
+    std::cerr << "Cannot read " << filename << "\n";
+    return 1;
+  }
+
+  std::vector<Point_3> input_points(points.begin(), points.end());
+
+  PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(erase_duplicate).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).number_of_iterations(15));
+
+
+  std::cout << "{" <<
+    "\"Nb_output_points\": \"" << points.size() << "\",\n" <<
+    "\"Nb_output_triangles\": \"" << triangles.size() << "\",\n" <<
+    "\"Is_2_manifold\": \"" << (PMP::orient_polygon_soup(points, triangles)?"True":"False") << "\",\n";
+  CGAL::Surface_mesh<Point_3> sm;
+  PMP::polygon_soup_to_polygon_mesh(points, triangles, sm);
+
+  std::cout << std::setprecision(17) <<
+   "\"Hausdorff_distance_output_to_input_(divide_by_bbox_diag)\": \"" << PMP::max_distance_to_triangle_mesh<CGAL::Parallel_if_available_tag>(input_points, sm) / diag_length << "\",\n" <<
+   "\"Closed_output\": \"" << (CGAL::is_closed(sm)?"True":"False") << "\",\n" <<
+   "\"Ouput_bound_a_volume\": \"" << (PMP::does_bound_a_volume(sm)?"True":"False") << "\"\n}"
+  << std::endl;
+
+  return 0;
+}

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Quality/run_quality.sh

@@ -0,0 +1,19 @@
+#!/bin/bash
+
+set -e
+
+if [ "$#" -lt 4 ]; then
+    echo "Usage: $0 <input_file> <timeout> [component_params...]"
+    exit 1
+fi
+
+INPUT_FILE=$1
+TIMEOUT=$2
+GRID_SIZE=$3
+ERASE_ALL_DUPLICATE=$4
+
+TMP_LOG=$(mktemp)
+timeout "$TIMEOUT"s quality_snap_polygon_soup "$INPUT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE"  > "$TMP_LOG"
+
+cat $TMP_LOG
+rm -f "$TMP_LOG"

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Robustness/robustness_snap_polygon_soup.cpp

@@ -0,0 +1,53 @@
+#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/triangulate_faces.h>
+#include <CGAL/IO/polygon_soup_io.h>
+
+#include <boost/container/small_vector.hpp>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef typename Kernel::Point_3 Point_3;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+enum EXIT_CODES { VALID_OUTPUT=0,
+                  INVALID_INPUT=1,
+                  ROUNDING_FAILED=2,
+                  SELF_INTERSECTING_OUTPUT=3,
+                  SIGSEGV=10,
+                  SIGSABRT=11,
+                  SIGFPE=12,
+                  TIMEOUT=13
+                  };
+
+int main(int argc, char** argv)
+{
+  if(argc<4){
+    std::cout << "Invalid argument" << std::endl;
+    return 1;
+  }
+
+  const std::string filename = std::string(argv[1]);
+  const int grid_size = std::stoi(std::string(argv[2]));
+  const bool erase_duplicate = std::stoi(argv[3])==1;
+
+  std::vector<Point_3> points;
+  std::vector<boost::container::small_vector<std::size_t, 3>> triangles;
+
+  if (!CGAL::IO::read_polygon_soup(filename, points, triangles) || points.size()==0 || triangles.size()==0)
+  {
+    return INVALID_INPUT;
+  }
+
+  PMP::repair_polygon_soup(points, triangles);
+  PMP::triangulate_polygons(points, triangles);
+
+  bool success=PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(erase_duplicate).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).number_of_iterations(15));
+
+  if(!success)
+    return ROUNDING_FAILED;
+  if( PMP::does_triangle_soup_self_intersect<CGAL::Parallel_if_available_tag>(points, triangles) )
+    return SELF_INTERSECTING_OUTPUT;
+
+  return VALID_OUTPUT;
+}

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Robustness/run_robustness.sh

@@ -0,0 +1,38 @@
+#!/bin/bash
+
+if [ "$#" -lt 4 ]; then
+    echo "Usage: $0 <input_file> <timeout> [component_params...]"
+    exit 1
+fi
+
+timeout_bis() {
+  timeout 5 sleep 10
+}
+
+INPUT_FILE=$1
+TIMEOUT=$2
+GRID_SIZE=$3
+ERASE_ALL_DUPLICATE=$4
+
+# Run with timeout, capture exit code
+timeout "--foreground" "$TIMEOUT"s robustness_snap_polygon_soup "$INPUT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE"
+EXIT_CODE=$?
+
+# Interpret exit codes
+declare -A TAGS
+TAGS[0]="VALID_OUTPUT"
+TAGS[1]="INPUT_IS_INVALID"
+TAGS[2]="ROUNDING_FAILED"
+TAGS[3]="SELF_INTERSECTING_OUTPUT"
+TAGS[139]="SIGSEGV"
+TAGS[11]="SIGSEGV"
+TAGS[6]="SIGABRT"
+TAGS[8]="SIGFPE"
+TAGS[132]="SIGILL"
+TAGS[124]="TIMEOUT"
+
+TAG_NAME=${TAGS[$EXIT_CODE]:-UNKNOWN}
+TAG_DESC=$([[ "$EXIT_CODE" -eq 0 ]] && echo "OK" || echo "Error")
+
+# Output JSON
+echo "{\"TAG_NAME\": \"$TAG_NAME\", \"TAG\": \"$TAG_DESC\"}"

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/benchmarking_snap_polygon_soup.sh

@@ -0,0 +1,116 @@
+#!/bin/bash
+
+# Temp directory for individual result JSONs
+TMP_RESULT_DIR=$(mktemp -d)
+
+# Job control
+JOBS=0
+MAX_JOBS=$NUM_THREADS
+
+# Function to process a single file
+process_file() {
+    INPUT_PATH="$1"
+    INPUT_ID=$(basename "$INPUT_PATH" | cut -d. -f1)
+    COMPONENT_NAME="$2"
+    PROJECT_DIR="$3"
+    TIMEOUT="$4"
+    OUTPUT_DIR="$5"
+    TMP_RESULT_FILE="$6"
+    GRID_SIZE="$7"
+    ERASE_ALL_DUPLICATE="$8"
+    {
+        echo "    \"$INPUT_ID\": {"
+        echo "      \"path\": \"$INPUT_PATH\","
+
+        PERF_OUTPUT=$(bash "$PROJECT_DIR/Performance/run_performance.sh" "$INPUT_PATH" "$TIMEOUT" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" 2>> "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Performance/$INPUT_ID.log")
+        echo "      \"Performance\": $PERF_OUTPUT,"
+
+        QUALITY_OUTPUT=$(bash "$PROJECT_DIR/Quality/run_quality.sh" "$INPUT_PATH" "$TIMEOUT" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" 2>> "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Quality/$INPUT_ID.log")
+        echo "      \"Quality\": $QUALITY_OUTPUT,"
+
+        ROBUST_OUTPUT=$(bash "$PROJECT_DIR/Robustness/run_robustness.sh" "$INPUT_PATH" "$TIMEOUT" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" 2>> "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Robustness/$INPUT_ID.log")
+        echo "      \"Robustness\": $ROBUST_OUTPUT"
+
+        echo "    }"
+    } > "$TMP_RESULT_FILE"
+}
+export -f process_file
+
+# Usage function
+usage() {
+    echo "Usage: $0 <project_dir> <input_data_dir> <output_results_dir> <timeout> <num_threads> [component_params...]"
+    exit 1
+}
+
+# Check parameters
+if [ "$#" -lt 5 ]; then
+    usage
+fi
+
+# Arguments
+PROJECT_DIR=$1
+INPUT_DIR=$2
+OUTPUT_DIR=$3
+TIMEOUT=$4
+NUM_THREADS=$5
+GRID_SIZE=$6
+ERASE_ALL_DUPLICATE=$7
+
+# Get component name from the project directory name
+COMPONENT_NAME=$(basename "$PROJECT_DIR")
+DATE_TAG=$(date +"%Y-%m-%d")
+TIMESTAMP=$(date +"%Y-%m-%d %H:%M:%S")
+RESULT_JSON="$OUTPUT_DIR/${COMPONENT_NAME}_results_${DATE_TAG}.json"
+
+# Compile
+# Do not forget to define CGAL_DIR
+cmake "$PROJECT_DIR" "-DCMAKE_BUILD_TYPE=Release" "-DCMAKE_CXX_FLAGS=-O3"
+make -j $NUM_THREADS
+
+# Prepare log directories
+mkdir -p "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Performance"
+mkdir -p "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Quality"
+mkdir -p "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Robustness"
+
+# Initialize JSON
+echo "{" > "$RESULT_JSON"
+echo "  \"$COMPONENT_NAME\": {" >> "$RESULT_JSON"
+echo "    \"Thingi10K\": {" >> "$RESULT_JSON"
+
+#process_file "$INPUT_DIR/100036.stl" "$COMPONENT_NAME" "$PROJECT_DIR" "$TIMEOUT" "$OUTPUT_DIR" "$TMP_RESULT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE"
+# Loop input files and spawn parallel jobs
+for INPUT_FILE in "$INPUT_DIR"/*; do
+    INPUT_ID=$(basename "$INPUT_FILE" | cut -d. -f1)
+    TMP_RESULT_FILE="$TMP_RESULT_DIR/$INPUT_ID.json"
+
+    process_file "$INPUT_FILE" "$COMPONENT_NAME" "$PROJECT_DIR" "$TIMEOUT" "$OUTPUT_DIR" "$TMP_RESULT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE"
+
+    ((JOBS+=1))
+    if [ "$JOBS" -ge "$NUM_THREADS" ]; then
+        wait
+        JOBS=0
+    fi
+done
+
+wait
+
+# Merge all partial JSONs
+echo "{" > "$RESULT_JSON"
+echo "  \"$COMPONENT_NAME\": {" >> "$RESULT_JSON"
+echo "    \"Thingi10K\": {" >> "$RESULT_JSON"
+
+FIRST_ENTRY=true
+for FILE in "$TMP_RESULT_DIR"/*.json; do
+    if [ "$FIRST_ENTRY" = true ]; then
+        FIRST_ENTRY=false
+    else
+        echo "," >> "$RESULT_JSON"
+    fi
+    cat "$FILE" >> "$RESULT_JSON"
+done
+
+echo "" >> "$RESULT_JSON"
+echo "    }," >> "$RESULT_JSON"
+echo "    \"finished_at\": \"$TIMESTAMP\"" >> "$RESULT_JSON"
+echo "  }" >> "$RESULT_JSON"
+echo "}" >> "$RESULT_JSON"

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/coplanar_cubes_autorefinement.cpp

@@ -0,0 +1,135 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/corefinement.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/Polygon_mesh_processing/transform.h>
+#include <CGAL/Real_timer.h>
+
+#include <boost/container/small_vector.hpp>
+
+#include <iostream>
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel K;
+typedef CGAL::Surface_mesh<K::Point_3>                      Mesh;
+typedef K::Point_3 Point_3;
+
+typedef CGAL::Simple_cartesian<double> Cartesian;
+typedef Cartesian::Point_3 Double_Point_3;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+namespace params = CGAL::parameters;
+
+struct Sphere_function {
+    double radius;
+    Sphere_function(double r) : radius(r) {}
+    Kernel::FT operator()(const Kernel::Point_3& p) const {
+        return p.x()*p.x() + p.y()*p.y() + p.z()*p.z() - radius*radius;
+    }
+};
+
+//Thanks Roberto!
+K::Aff_transformation_3
+random_rotation(CGAL::Random &gen)
+{
+  double a=gen.get_double(0,2*CGAL_PI);
+
+  double ca = cos(a);
+  double sa = sin(a);
+
+  K::Aff_transformation_3 aff(1, 0, 0,
+                              0, ca,-sa,
+                              0, sa, ca);
+  std::cout << "Rotation by " << a << std::endl;
+  return aff;
+}
+
+int main(int argc, char** argv)
+{
+  Mesh cube;
+  std::vector<Point_3> points;
+  std::vector<boost::container::small_vector<std::size_t, 3>> faces;
+
+  CGAL::make_hexahedron(
+      Point_3(-1,-1,-1),
+      Point_3(1,-1,-1),
+      Point_3(1,1,-1),
+      Point_3(-1,1,-1),
+      Point_3(-1,1,1),
+      Point_3(-1,-1,1),
+      Point_3(1,-1,1),
+      Point_3(1,1,1),
+      cube,
+      CGAL::parameters::do_not_triangulate_faces(false)
+      );
+
+  std::cout << "Iterative intersection of rotative cubes with snapping" << std::endl;
+
+  int i=0;
+  CGAL::Random random_gen = argc == 1 ? CGAL::get_default_random() : CGAL::Random(std::stoi(argv[1]));
+
+  Mesh inter=cube;
+  while(true)
+  {
+    std::cout << "Iteration " << i << std::endl;
+
+    CGAL::Real_timer t;
+    t.start();
+
+    std::cout << "Add a randomly rotated cube to the scene" << std::endl;
+    Mesh rotated_cube=cube;
+    PMP::transform(random_rotation(random_gen), rotated_cube);
+
+    std::cout << "compute_intersection" << std::endl;
+    bool OK = PMP::corefine_and_compute_intersection(inter, rotated_cube, inter);
+
+    if(!OK){
+      std::cout << "No manifold, stop experiment" << std::endl;
+      exit(0);
+    }
+
+    points.clear();
+    faces.clear();
+    PMP::polygon_mesh_to_polygon_soup(inter, points, faces);
+
+    std::cout << "Snapped the points on double" << std::endl;
+    bool success=PMP::autorefine_triangle_soup(points, faces, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag()));
+    t.stop();
+    if(!success){
+      std::cout << "Rounding failed" << std::endl;
+      exit(0);
+    }
+
+    //dump model every 100 iterations
+    if(i%100==0){
+      std::cout << "Dump model" << std::endl;
+      std::vector<Double_Point_3> double_points;
+      for(auto &p: points)
+        double_points.emplace_back(CGAL::to_double(p.x()),CGAL::to_double(p.y()),CGAL::to_double(p.z()));
+      std::ofstream outfile("cubes_"+std::to_string(i)+".off");
+      outfile.precision(17);
+      outfile << "OFF\n" << points.size() << " " << faces.size() << " 0\n";
+      for(auto p : points)
+      outfile << p.x() << " " << p.y() << " " << p.z() << std::endl;
+
+      for(auto &t : faces)
+          outfile << "3" << " " << t[0] << " " << t[1] << " " << t[2] << std::endl;
+      outfile.close();//
+    }
+
+    std::cout << "#points = " << points.size() << " and #triangles = " << faces.size() << " in " << t.time() << " sec.\n\n" << std::endl;
+
+    inter.clear();
+    PMP::polygon_soup_to_polygon_mesh(points, faces, inter);
+    ++i;
+  }
+
+  return 0;
+}
+

Polygon_mesh_processing/benchmark/Polygon_mesh_processing/rotated_cubes_autorefinement.cpp

@@ -0,0 +1,137 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/corefinement.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/Polygon_mesh_processing/transform.h>
+#include <CGAL/Real_timer.h>
+
+#include <boost/container/small_vector.hpp>
+
+#include <iostream>
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel K;
+typedef CGAL::Surface_mesh<K::Point_3>                      Mesh;
+typedef K::Point_3 Point_3;
+
+typedef CGAL::Simple_cartesian<double> Cartesian;
+typedef Cartesian::Point_3 Double_Point_3;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+namespace params = CGAL::parameters;
+
+struct Sphere_function {
+    double radius;
+    Sphere_function(double r) : radius(r) {}
+    Kernel::FT operator()(const Kernel::Point_3& p) const {
+        return p.x()*p.x() + p.y()*p.y() + p.z()*p.z() - radius*radius;
+    }
+};
+
+//Thanks Roberto!
+K::Aff_transformation_3
+random_rotation(CGAL::Random &gen)
+{
+  double a=gen.get_double(0,2*CGAL_PI);
+  double b=gen.get_double(0,2*CGAL_PI);
+  double c=gen.get_double(0,2*CGAL_PI);
+
+  double ca = cos(a), cb = cos(b), cc = cos(c);
+  double sa = sin(a), sb = sin(b), sc = sin(c);
+
+  K::Aff_transformation_3 aff(cb * cc, cc* sa* sb - ca * sc, ca* cc* sb + sa * sc,
+                              cb* sc, ca* cc + sa * sb * sc, ca* sb* sc - cc * sa,
+                              -sb, cb* sa, ca* cb);
+  std::cout << "Rotation by " << a << " " << b << " " << c << std::endl;
+  return aff;
+}
+
+int main(int argc, char** argv)
+{
+  Mesh cube;
+  std::vector<Point_3> points;
+  std::vector<boost::container::small_vector<std::size_t, 3>> faces;
+
+  CGAL::make_hexahedron(
+      Point_3(-1,-1,-1),
+      Point_3(1,-1,-1),
+      Point_3(1,1,-1),
+      Point_3(-1,1,-1),
+      Point_3(-1,1,1),
+      Point_3(-1,-1,1),
+      Point_3(1,-1,1),
+      Point_3(1,1,1),
+      cube,
+      CGAL::parameters::do_not_triangulate_faces(false)
+      );
+
+  std::cout << "Iterative intersection of rotative cubes with snapping" << std::endl;
+
+  int i=0;
+  CGAL::Random random_gen = argc == 1 ? CGAL::get_default_random() : CGAL::Random(std::stoi(argv[1]));
+
+  Mesh inter=cube;
+  while(true)
+  {
+    std::cout << "Iteration " << i << std::endl;
+
+    CGAL::Real_timer t;
+    t.start();
+
+    std::cout << "Add a randomly rotated cube to the scene" << std::endl;
+    Mesh rotated_cube=cube;
+    PMP::transform(random_rotation(random_gen), rotated_cube);
+
+    std::cout << "compute_intersection" << std::endl;
+    bool OK = PMP::corefine_and_compute_intersection(inter, rotated_cube, inter);
+
+    if(!OK){
+      std::cout << "No manifold, stop experiment" << std::endl;
+      exit(0);
+    }
+
+    points.clear();
+    faces.clear();
+    PMP::polygon_mesh_to_polygon_soup(inter, points, faces);
+
+    std::cout << "Snapped the points on double" << std::endl;
+    bool success=PMP::autorefine_triangle_soup(points, faces, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag()));
+    t.stop();
+    if(!success){
+      std::cout << "Rounding failed" << std::endl;
+      exit(0);
+    }
+
+    //dump model every 100 iterations
+    if(i%100==0){
+      std::cout << "Dump model" << std::endl;
+      std::vector<Double_Point_3> double_points;
+      for(auto &p: points)
+        double_points.emplace_back(CGAL::to_double(p.x()),CGAL::to_double(p.y()),CGAL::to_double(p.z()));
+      std::ofstream outfile("cubes_"+std::to_string(i)+".off");
+      outfile.precision(17);
+      outfile << "OFF\n" << points.size() << " " << faces.size() << " 0\n";
+      for(auto p : points)
+      outfile << p.x() << " " << p.y() << " " << p.z() << std::endl;
+
+      for(auto &t : faces)
+          outfile << "3" << " " << t[0] << " " << t[1] << " " << t[2] << std::endl;
+      outfile.close();//
+    }
+
+    std::cout << "#points = " << points.size() << " and #triangles = " << faces.size() << " in " << t.time() << " sec.\n\n" << std::endl;
+
+    inter.clear();
+    PMP::polygon_soup_to_polygon_mesh(points, faces, inter);
+    ++i;
+  }
+
+  return 0;
+}
+

Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h

@@ -23,5 +23,7 @@ public:
   /// called for each subtriangle created from a triangle with intersection, `tgt_id` is the position in the triangle container after calling
   /// `autorefine_triangle_soup()` of the subtriangle, while `src_id` was the position of the original support triangle before calling the function.
   void new_subtriangle(std::size_t tgt_id, std::size_t src_id);
+  /// called for each input triangle absent in the output because it was degenerate. `src_id` is the position of that triangle in the input range. Additionally, if `apply_iterative_snap_rounding()` is set to `true`, some extra triangles might become degenerate and will be absent in the output. Those triangles are also reported by this function.
+  void delete_triangle(std::size_t src_id);
 /// @}
 };

Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt

@@ -4,7 +4,7 @@ namespace CGAL {
 \anchor Chapter_PolygonMeshProcessing
 
 \cgalAutoToc
-\authors David Coeurjolly, Jaques-Olivier Lachaud,  Konstantinos Katrioplas, Sébastien Loriot, Ivan Pađen, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, Sébastien Valette, and Ilker %O. Yaz
+\authors David Coeurjolly, Jaques-Olivier Lachaud,  Konstantinos Katrioplas, Sébastien Loriot, Ivan Pađen, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, Sébastien Valette, Léo Valque, and Ilker %O. Yaz
 
 \image html neptun_head.jpg
 \image latex neptun_head.jpg
@@ -911,8 +911,10 @@ would then also includes overlaps of duplicated points.
 The function `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` provides a way to refine a triangle soup
 using the intersections of the triangles from the soup. In particular, if some points are duplicated they will be
 merged. Note that if a kernel with exact predicates but inexact constructions is used, some new self-intersections
-might be introduced due to rounding issues of points coordinates.
+might be introduced due to the rounding of the coordinates of intersection points. The `apply_iterative_snap_rounding` option can be used to resolve this issue.
-To guarantee that the triangle soup is free from self-intersections, a kernel with exact constructions must be used.
+When set to `true`, it ensures that the coordinates are rounded to fit in `double` with potential additional subdivisions,
+preventing any self-intersections from occurring.
+
 
 \subsection PMPRemoveCapsNeedles Removal of Almost Degenerate Triangle Faces
 Triangle faces of a mesh made up of almost collinear points are badly shaped elements that
@@ -1477,6 +1479,7 @@ supervision of Sébastien Valette and Sébastien Loriot, who later finalized the
 The implementation is based on \cgalCite{cgal:vcp-grtmmdvd-08}. and
 preceding work. <a href="https://www.creatis.insa-lyon.fr/~valette/public/project/acvd/">ACVD's implementation</a> was also
 used as a reference during the project.
+The `apply_iterative_snap_rounding` option for autorefinement was implemented in 2025, by Léo Valque, based on his work with Sylvain Lazard \cgalCite{lazard:hal-04907149}.
 
 */
 } /* namespace CGAL */

Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt

@@ -53,4 +53,5 @@
 \example Polygon_mesh_processing/acvd_remeshing_example.cpp
 \example Polygon_mesh_processing/sample_example.cpp
 \example Polygon_mesh_processing/soup_autorefinement.cpp
+\example Polygon_mesh_processing/snap_polygon_soup.cpp
 */

Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

@@ -58,6 +58,7 @@ create_single_source_cgal_program("isotropic_remeshing_with_custom_sizing_exampl
 create_single_source_cgal_program("isotropic_remeshing_with_allow_move.cpp")
 create_single_source_cgal_program("triangle_mesh_autorefinement.cpp")
 create_single_source_cgal_program("soup_autorefinement.cpp")
+create_single_source_cgal_program("snap_polygon_soup.cpp")
 
 find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater)
 include(CGAL_Eigen3_support)
@@ -140,6 +141,7 @@ if(TARGET CGAL::TBB_support)
   target_link_libraries(hausdorff_distance_remeshing_example PRIVATE CGAL::TBB_support)
   target_link_libraries(hausdorff_bounded_error_distance_example PRIVATE CGAL::TBB_support)
   target_link_libraries(soup_autorefinement PRIVATE CGAL::TBB_support)
+  target_link_libraries(snap_polygon_soup PRIVATE CGAL::TBB_support)
 
   create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp")
   target_link_libraries(corefinement_parallel_union_meshes PRIVATE CGAL::TBB_support)

Polygon_mesh_processing/examples/Polygon_mesh_processing/snap_polygon_soup.cpp

@@ -0,0 +1,62 @@
+#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>
+
+#include <boost/container/small_vector.hpp>
+
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char** argv)
+{
+  const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
+                                         : std::string(argv[1]);
+
+  const int grid_size = argc <= 2 ? 23
+                                  : std::stoi(std::string(argv[2]));
+
+  const std::string out_file = "rounded_soup.off";
+
+  std::vector<typename Kernel::Point_3> points;
+  std::vector<boost::container::small_vector<std::size_t, 3>> triangles;
+
+  std::cout << "Snap rounding on " << filename << "\n";
+  if (!CGAL::IO::read_polygon_soup(filename, points, triangles))
+  {
+    std::cerr << "Cannot read " << filename << "\n";
+    return 1;
+  }
+
+  PMP::repair_polygon_soup(points, triangles);
+  PMP::triangulate_polygons(points, triangles);
+
+  std::cout << "#points = " << points.size() << " and #triangles = " << triangles.size() << std::endl;
+  std::cout << "Is 2-manifold: " << PMP::is_polygon_soup_a_polygon_mesh(triangles) << std::endl;
+
+  std::vector<std::pair<std::size_t, std::size_t>> pairs_of_intersecting_triangles;
+  PMP::triangle_soup_self_intersections(points, triangles, std::back_inserter(pairs_of_intersecting_triangles));
+  std::cout << "Nb of pairs of intersecting triangles: " << pairs_of_intersecting_triangles.size() << std::endl;
+
+  CGAL::Real_timer t;
+  t.start();
+  bool success=PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(false).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).number_of_iterations(15));
+  t.stop();
+
+  std::cout << "\nOutput:" << std::endl;
+  std::cout << "#points = " << points.size() << " and #triangles = " << triangles.size() << " in " << t.time() << " sec." << std::endl;
+  if(success)
+    std::cout << "Does self-intersect: " << PMP::does_triangle_soup_self_intersect<CGAL::Parallel_if_available_tag>(points, triangles) << std::endl;
+  else
+    std::cout << "ROUNDING FAILED" << std::endl;
+
+  CGAL::IO::write_polygon_soup(out_file, points, triangles, CGAL::parameters::stream_precision(17));
+  std::cout << "Is 2-manifold: " << PMP::orient_polygon_soup(points, triangles) << "\n\n"  << std::endl;
+  return 0;
+}

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h

@@ -33,7 +33,6 @@
 #include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
 #include <CGAL/Polygon_mesh_processing/polygon_mesh_to_polygon_soup.h>
 
-
 #ifdef CGAL_PMP_AUTOREFINE_USE_DEFAULT_VERBOSE
 #define CGAL_PMP_AUTOREFINE_VERBOSE(X) std::cout << X << "\n";
 #endif
@@ -90,6 +89,7 @@ struct Default_visitor
   inline void number_of_output_triangles(std::size_t /*nbt*/) {}
   inline void verbatim_triangle_copy(std::size_t /*tgt_id*/, std::size_t /*src_id*/) {}
   inline void new_subtriangle(std::size_t /*tgt_id*/, std::size_t /*src_id*/) {}
+  inline void delete_triangle(std::size_t /*src_id*/) {}
 };
 
 } // end of Autorefinement visitor
@@ -975,59 +975,19 @@ void generate_subtriangles(std::size_t ti,
 } // end of autorefine_impl
 #endif
 
-/**
+namespace autorefine_impl{
-* \ingroup PMP_corefinement_grp
+// Forward declaration
-*
+struct Wrap_snap_visitor;
-* refines a soup of triangles so that no pair of triangles intersects.
+
-* Output triangles may share a common edge or a common vertex (but with the same indexed position in `points`).
+template <typename PointRange, typename PolygonRange, class NamedParameters = parameters::Default_named_parameters>
-* Note that points in `soup_points` can only be added (intersection points) at the end of the container, with the initial order preserved.
+bool polygon_soup_snap_rounding(PointRange &points,
-* Note that if `soup_points` contains two or more identical points then only the first copy (following the order in the `soup_points`)
+                                PolygonRange &triangles,
-* will be used in `soup_triangles`.
+                                const NamedParameters& np = parameters::default_values());
-* `soup_triangles` will be updated to contain both the input triangles and the new subdivided triangles. Degenerate triangles will be removed.
+
-* Also triangles in `soup_triangles` will be triangles without intersection first, followed by triangles coming from a subdivision induced
-* by an intersection. The named parameter `visitor()` can be used to track
-*
-* @tparam PointRange a model of the concept `RandomAccessContainer`
-* whose value type is the point type
-* @tparam TriangleRange a model of the concepts `RandomAccessContainer`, `BackInsertionSequence` and `Swappable`, whose
-* value type is a model of the concept `RandomAccessContainer` whose value type is convertible to `std::size_t` and that
-* is constructible from an `std::initializer_list<std::size_t>` of size 3.
-* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
-*
-* @param soup_points points of the soup of polygons
-* @param soup_triangles each element in the range describes a triangle using the indexed position of the points in `soup_points`
-* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
-*
-* \cgalNamedParamsBegin
-*   \cgalParamNBegin{concurrency_tag}
-*     \cgalParamDescription{a tag indicating if the task should be done using one or several threads.}
-*     \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`}
-*     \cgalParamDefault{`CGAL::Sequential_tag`}
-*   \cgalParamNEnd
-*   \cgalParamNBegin{point_map}
-*     \cgalParamDescription{a property map associating points to the elements of the range `soup_points`}
-*     \cgalParamType{a model of `ReadWritePropertyMap` whose value type is a point type}
-*     \cgalParamDefault{`CGAL::Identity_property_map`}
-*   \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 point type.}
-*   \cgalParamNEnd
-*   \cgalParamNBegin{visitor}
-*     \cgalParamDescription{a visitor used to track the creation of new faces}
-*     \cgalParamType{a class model of `PMPAutorefinementVisitor`}
-*     \cgalParamDefault{`Autorefinement::Default_visitor`}
-*     \cgalParamExtra{The visitor will be copied.}
-*   \cgalParamNEnd
-* \cgalNamedParamsEnd
-*
-*/
 template <class PointRange, class TriangleRange, class NamedParameters = parameters::Default_named_parameters>
-void autorefine_triangle_soup(PointRange& soup_points,
+bool autorefine_triangle_soup(PointRange& soup_points,
-                              TriangleRange& soup_triangles,
+                               TriangleRange& soup_triangles,
-                              const NamedParameters& np = parameters::default_values())
+                               const NamedParameters& np = parameters::default_values())
 {
   using parameters::choose_parameter;
   using parameters::get_parameter;
@@ -1050,7 +1010,6 @@ void autorefine_triangle_soup(PointRange& soup_points,
   > ::type Visitor;
   Visitor visitor(choose_parameter<Visitor>(get_parameter(np, internal_np::visitor)));
 
-
   constexpr bool parallel_execution = std::is_same_v<Parallel_tag, Concurrency_tag>;
 
 #ifndef CGAL_LINKED_WITH_TBB
@@ -1077,7 +1036,7 @@ void autorefine_triangle_soup(PointRange& soup_points,
       for(std::size_t i=0; i<soup_triangles.size(); ++i)
         visitor.verbatim_triangle_copy(i, i);
     }
-    return;
+    return true;
   }
 
   // mark degenerate faces so that we can ignore them
@@ -1455,15 +1414,17 @@ void autorefine_triangle_soup(PointRange& soup_points,
   std::vector<std::size_t> tri_inter_ids_inverse(triangles.size());
   for (Input_TID f=0; f<soup_triangles.size(); ++f)
   {
-    if (is_degen[f]) continue; //skip degenerate faces
+    if (is_degen[f])
+    {
+      visitor.delete_triangle(f);
+      continue; //skip degenerate faces
+    }
 
     int tiid = tri_inter_ids[f];
     if (tiid == -1)
     {
       visitor.verbatim_triangle_copy(soup_triangles_out.size(), f);
-      soup_triangles_out.push_back(
+      soup_triangles_out.push_back(soup_triangles[f]);
-        {soup_triangles[f][0], soup_triangles[f][1], soup_triangles[f][2]}
-      );
     }
     else
     {
@@ -1529,22 +1490,41 @@ void autorefine_triangle_soup(PointRange& soup_points,
       [&](const tbb::blocked_range<size_t>& r) {
         for (size_t ti = r.begin(); ti != r.end(); ++ti) {
           const std::array<EK::Point_3, 3>& t = new_triangles[ti].first;
-          visitor.new_subtriangle(offset+ti, tri_inter_ids_inverse[new_triangles[ti].second]);
           triangle_buffer[ti] = CGAL::make_array(concurrent_get_point_id(t[0]), concurrent_get_point_id(t[1]), concurrent_get_point_id(t[2]));
         }
       }
     );
-    tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
+
-      [&](const tbb::blocked_range<size_t>& r) {
+    // The constexpr  was initially inside the lambda, but that did not compile  with VC 2017
-          for (size_t ti = r.begin(); ti != r.end(); ++ti)
+    if constexpr(std::is_same_v<Visitor, Wrap_snap_visitor>){
-          {
+      tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
-            soup_triangles_out[offset + ti] =
+        [&](const tbb::blocked_range<size_t>& r) {
-              { triangle_buffer[ti][0]->second,
+            for (size_t ti = r.begin(); ti != r.end(); ++ti)
-                triangle_buffer[ti][1]->second,
+            {
-                triangle_buffer[ti][2]->second };
+              soup_triangles_out[offset + ti] =
+                { triangle_buffer[ti][0]->second,
+                  triangle_buffer[ti][1]->second,
+                  triangle_buffer[ti][2]->second };
+              visitor.new_subdivision(soup_triangles_out[offset + ti], soup_triangles[tri_inter_ids_inverse[new_triangles[ti].second]]);
+            }
           }
-        }
+      );
-    );
+    }
+    else
+    {
+      tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
+        [&](const tbb::blocked_range<size_t>& r) {
+            for (size_t ti = r.begin(); ti != r.end(); ++ti)
+            {
+              soup_triangles_out[offset + ti] =
+                { triangle_buffer[ti][0]->second,
+                  triangle_buffer[ti][1]->second,
+                  triangle_buffer[ti][2]->second };
+              visitor.new_subtriangle(offset+ti, tri_inter_ids_inverse[new_triangles[ti].second]);
+            }
+          }
+      );
+    }
 #else
     //option 2 (without mutex)
     /// Lambda concurrent_get_point_id()
@@ -1591,17 +1571,33 @@ void autorefine_triangle_soup(PointRange& soup_points,
       }
     );
 
-    tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
+    // The constexpr  was initially inside the lambda, but that did not compile  with VC 2017
-      [&](const tbb::blocked_range<size_t>& r) {
+    if constexpr(std::is_same_v<Visitor, Wrap_snap_visitor>){
-          for (size_t ti = r.begin(); ti != r.end(); ++ti)
+      tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
-          {
+        [&](const tbb::blocked_range<size_t>& r) {
-            soup_triangles_out[offset + ti] =
+            for (size_t ti = r.begin(); ti != r.end(); ++ti)
-              { triangle_buffer[ti][0]->second,
+            {
-                triangle_buffer[ti][1]->second,
+              soup_triangles_out[offset + ti] =
-                triangle_buffer[ti][2]->second };
+                { triangle_buffer[ti][0]->second,
+                  triangle_buffer[ti][1]->second,
+                  triangle_buffer[ti][2]->second };
+                visitor.new_subdivision(soup_triangles_out[offset + ti], soup_triangles[tri_inter_ids_inverse[new_triangles[ti].second]]);
+            }
           }
-        }
+      );
-    );
+    }else{
+      tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
+        [&](const tbb::blocked_range<size_t>& r) {
+            for (size_t ti = r.begin(); ti != r.end(); ++ti)
+            {
+              soup_triangles_out[offset + ti] =
+                { triangle_buffer[ti][0]->second,
+                  triangle_buffer[ti][1]->second,
+                  triangle_buffer[ti][2]->second };
+            }
+          }
+      );
+    }
 #endif
   }
   else
@@ -1613,10 +1609,13 @@ void autorefine_triangle_soup(PointRange& soup_points,
     soup_triangles_out.reserve(offset + new_triangles.size());
     for (const std::pair<std::array<EK::Point_3,3>, std::size_t>& t_and_id : new_triangles)
     {
-      visitor.new_subtriangle(soup_triangles_out.size(), tri_inter_ids_inverse[t_and_id.second]);
       soup_triangles_out.push_back({ get_point_id(t_and_id.first[0]),
                                      get_point_id(t_and_id.first[1]),
                                      get_point_id(t_and_id.first[2]) });
+      if constexpr(std::is_same_v<Visitor, Wrap_snap_visitor>)
+        visitor.new_subdivision(soup_triangles_out[soup_triangles_out.size()-1], soup_triangles[tri_inter_ids_inverse[t_and_id.second]]);
+      else
+        visitor.new_subtriangle(soup_triangles_out.size()-1, tri_inter_ids_inverse[t_and_id.second]);
     }
   }
 
@@ -1642,8 +1641,108 @@ void autorefine_triangle_soup(PointRange& soup_points,
   swap(soup_triangles, soup_triangles_out);
 
   CGAL_PMP_AUTOREFINE_VERBOSE("done");
+  return true;
 }
 
+} // end of autorefine_impl
+
+/**
+* \ingroup PMP_corefinement_grp
+*
+* refines a soup of triangles so that no pair of triangles intersects.
+* Output triangles may share a common edge or a common vertex (but with the same indexed position in `points`).
+* Note that if `apply_iterative_snap_rounding` option is set to `false`, points in `soup_points` can only be added (intersection points) at the end of the container, with the initial order preserved.
+* Note that if `soup_points` contains two or more identical points then only the first copy (following the order in the `soup_points`)
+* will be used in `soup_triangles`. if `apply_iterative_snap_rounding` is set to `true`, all duplicates points are removed.
+* `soup_triangles` will be updated to contain both the input triangles and the new subdivided triangles. Degenerate triangles will be removed.
+* Also if `apply_iterative_snap_rounding` option is set to `false`, triangles in `soup_triangles` will be triangles without intersection first, followed by triangles coming from a subdivision induced
+* by an intersection. The named parameter `visitor()` can be used to track the creation and removal of triangles independantly of
+* the `apply_iterative_snap_rounding` option.
+* If the `apply_iterative_snap_rounding` parameter is set to `true`, the coordinates of the vertices are rounded to fit within the precision of a double-precision floating point,
+* while trying to make the triangle soup free of intersections. The `snap_grid_size()` parameter limits the drift of the snapped vertices.
+* A smaller value is more likely to output an intersection free output and perform more vertex collapses, but it may increase the Hausdorff distance from the input.
+*
+* @tparam PointRange a model of the concept `RandomAccessContainer`
+* whose value type is the point type
+* @tparam TriangleRange a model of the concepts `RandomAccessContainer`, `BackInsertionSequence` and `Swappable`, whose
+* value type is a model of the concept `RandomAccessContainer` whose value type is convertible to `std::size_t` and that
+* is constructible from an `std::initializer_list<std::size_t>` of size 3.
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* @param soup_points points of the soup of polygons
+* @param soup_triangles each element in the range describes a triangle using the indexed position of the points in `soup_points`
+* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{concurrency_tag}
+*     \cgalParamDescription{a tag indicating if the task should be done using one or several threads.}
+*     \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`}
+*     \cgalParamDefault{`CGAL::Sequential_tag`}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{point_map}
+*     \cgalParamDescription{a property map associating points to the elements of the range `soup_points`}
+*     \cgalParamType{a model of `ReadWritePropertyMap` whose value type is a point type}
+*     \cgalParamDefault{`CGAL::Identity_property_map`}
+*   \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 point type.}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{visitor}
+*     \cgalParamDescription{a visitor used to track the creation of new faces}
+*     \cgalParamType{a class model of `PMPAutorefinementVisitor`}
+*     \cgalParamDefault{`Autorefinement::Default_visitor`}
+*     \cgalParamExtra{The visitor will be copied.}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{apply_iterative_snap_rounding}
+*     \cgalParamDescription{Enable the rounding of the coordinates so that they fit in doubles.}
+*     \cgalParamType{Boolean}
+*     \cgalParamDefault{false}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{snap_grid_size}
+*     \cgalParamDescription{A value `gs` used to scale the points to `[-2^gs, 2^gs]` before rounding them on integers. Used only if `apply_iterative_snap_rounding()` is set to `true`}
+*     \cgalParamType{`unsigned int`}
+*     \cgalParamDefault{23}
+*     \cgalParamExtra{Must be lower than 52.}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{number_of_iterations}
+*     \cgalParamDescription{Maximum number of iterations performed by the snap rounding algorithm. Used only if `apply_iterative_snap_rounding` is true.}
+*     \cgalParamType{unsigned int}
+*     \cgalParamDefault{5}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+* \return `true` if `apply_iterative_snap_rounding` is set to `false`.
+* Otherwise, returns `true` if the modified triangle soup is free of self-intersections,
+* or `false` if the algorithm was unable to produce such a result within the allowed number of iterations.
+* In the latter case, the output triangle soup represents a partial result from the final iteration,
+* with no guarantee of its validity.
+*
+*/
+template <class PointRange, class TriangleRange, class NamedParameters = parameters::Default_named_parameters>
+bool autorefine_triangle_soup(PointRange& soup_points,
+                              TriangleRange& soup_triangles,
+                              const NamedParameters& np = parameters::default_values())
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  //Dispatch the execution according the apply_iterative_snap_rounding parameter
+  const bool do_snap = choose_parameter(get_parameter(np, internal_np::apply_iterative_snap_rounding), false);
+  if(do_snap)
+  {
+    CGAL_PMP_AUTOREFINE_VERBOSE("Snap polygon soup");
+    return autorefine_impl::polygon_soup_snap_rounding(soup_points, soup_triangles, np);
+  }
+  else
+  {
+    return autorefine_impl::autorefine_triangle_soup(soup_points, soup_triangles, np);
+  }
+}
+
+
 /**
  * \ingroup PMP_corefinement_grp
  * refines a triangle mesh so that no triangles intersects in their interior.
@@ -1717,4 +1816,6 @@ autorefine(      TriangleMesh& tm,
 #endif
 #endif
 
+#include <CGAL/Polygon_mesh_processing/internal/triangle_soup_snap_rounding.h>
+
 #endif // CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/triangle_soup_snap_rounding.h

@@ -0,0 +1,585 @@
+// Copyright (c) 2025 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Léo Valque
+
+#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_POLYGON_SOUP_SNAP_ROUNDING_H
+#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_POLYGON_SOUP_SNAP_ROUNDING_H
+
+#include <CGAL/license/Polygon_mesh_processing/geometric_repair.h>
+
+#include <CGAL/Algebraic_structure_traits.h>
+#include <CGAL/number_utils.h>
+
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <CGAL/Fraction_traits.h>
+#include <CGAL/Lazy_exact_nt.h>
+
+namespace CGAL
+{
+
+namespace Polygon_mesh_processing
+{
+
+namespace autorefine_impl
+{
+
+// Certified ceil function for exact number types
+template <class NT> double double_ceil(const Lazy_exact_nt< NT > &x);
+template <class NT> double double_ceil(const NT &x);
+
+template <class NT>
+double double_ceil(const Lazy_exact_nt< NT > &x){
+  // If both sides are in the same ceil, return this ceil
+  double ceil_left=std::ceil(to_interval(x).first);
+  if(ceil_left==std::ceil(to_interval(x).second))
+    return ceil_left;
+  // If not refine the interval by contracting the DAG and try again
+  x.exact();
+  ceil_left=std::ceil(to_interval(x).first);
+  if(ceil_left==std::ceil(to_interval(x).second))
+    return ceil_left;
+  // If not return the ceil of the exact value
+  return double_ceil( x.exact());
+};
+
+template <class NT>
+double double_ceil(const NT &x){
+  using FT = Fraction_traits<NT>;
+  if constexpr(FT::Is_fraction::value){
+    // If NT is a fraction, the ceil value is the result of the euclidian division of the numerator and the denominator.
+    typename FT::Numerator_type num, r, e;
+    typename FT::Denominator_type denom;
+    typename FT::Decompose()(x,num,denom);
+    div_mod(num, denom, r, e);
+    if((r>=0) && e!=0) //If the result is positive, the ceil value is one above
+      return to_double(r+1);
+    return to_double(r);
+  } else {
+    // Return the ceil of the approximation
+    return std::ceil(to_double(x));
+  }
+};
+
+// Provides an index to the range, which is used by the visitor to track the correspondence between the input and the output
+template <typename Range>
+class Indexes_range : public Range{
+  typedef std::remove_cv_t<typename std::iterator_traits<typename Range::iterator>::value_type> Value_type;
+public:
+  typedef typename Range::const_iterator const_iterator;
+  typedef typename Range::iterator iterator;
+
+  Indexes_range() = default;
+  Indexes_range(const std::initializer_list<std::size_t>  &l): Range(l), m_id(0), modified(true){}
+  Indexes_range(Range &p): Range(p), modified(true){}
+  Indexes_range(Range &&p): Range(std::move(p)), modified(true){}
+  Indexes_range(Range &p, std::size_t id): Range(p), m_id(id),modified(false){}
+  Indexes_range(Range &&p, std::size_t id): Range(std::move(p)), m_id(id),modified(false){}
+
+  inline std::size_t id() const { return m_id; }
+  inline void set_id(std::size_t id){ m_id=id; }
+  inline bool was_modified() const { return modified; }
+
+private:
+  std::size_t m_id;
+  bool modified;
+};
+
+// Repair_polygon_soup without remove_pinched_polygons since our polygon are triangles
+template <typename PointRange, typename PolygonRange,
+          typename Polygon = typename ::CGAL::internal::Polygon_types<PointRange, PolygonRange>::Polygon_3>
+struct Triangle_soup_fixer
+{
+  template <typename NamedParameters>
+  void operator()(PointRange& points,
+                  PolygonRange& polygons,
+                  const NamedParameters& np) const
+  {
+    using parameters::get_parameter;
+    using parameters::choose_parameter;
+
+    using Traits = typename GetPolygonGeomTraits<PointRange, PolygonRange, NamedParameters>::type;
+    Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
+
+    CGAL::Polygon_mesh_processing::merge_duplicate_points_in_polygon_soup(points, polygons, np);
+    CGAL::Polygon_mesh_processing::internal::simplify_polygons_in_polygon_soup(points, polygons, traits);
+    CGAL::Polygon_mesh_processing::internal::remove_invalid_polygons_in_polygon_soup(points, polygons);
+    CGAL::Polygon_mesh_processing::merge_duplicate_polygons_in_polygon_soup(points, polygons, np);
+    CGAL::Polygon_mesh_processing::remove_isolated_points_in_polygon_soup(points, polygons);
+  }
+};
+
+// Specialization for array and Indexes_range of array
+template <typename PointRange, typename PolygonRange, typename PID, std::size_t N>
+struct Triangle_soup_fixer<PointRange, PolygonRange, std::array<PID, N> >
+{
+  template <typename NamedParameters>
+  void operator()(PointRange& points,
+                  PolygonRange& polygons,
+                  const NamedParameters& np) const
+  {
+    repair_polygon_soup(points, polygons, np);
+  }
+};
+
+template <typename PointRange, typename PolygonRange, typename PID, std::size_t N>
+struct Triangle_soup_fixer<PointRange, PolygonRange, Indexes_range< std::array<PID, N> > >
+{
+  template <typename NamedParameters>
+  void operator()(PointRange& points,
+                  PolygonRange& polygons,
+                  const NamedParameters& np) const
+  {
+    using parameters::get_parameter;
+    using parameters::choose_parameter;
+
+    using Traits = typename GetPolygonGeomTraits<PointRange, PolygonRange, NamedParameters>::type;
+    Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
+
+    CGAL::Polygon_mesh_processing::merge_duplicate_points_in_polygon_soup(points, polygons, np);
+    CGAL::Polygon_mesh_processing::internal::remove_invalid_polygons_in_array_polygon_soup(points, polygons, traits);
+    CGAL::Polygon_mesh_processing::merge_duplicate_polygons_in_polygon_soup(points, polygons, np);
+    CGAL::Polygon_mesh_processing::remove_isolated_points_in_polygon_soup(points, polygons);
+  }
+};
+
+template <typename PointRange, typename PolygonRange, typename NamedParameters>
+void repair_triangle_soup(PointRange& points,
+                          PolygonRange& polygons,
+                          const NamedParameters& np)
+{
+  Triangle_soup_fixer<PointRange, PolygonRange> fixer;
+  fixer(points, polygons, np);
+}
+
+// A visitor of Autorefinement to track the correspondance between input and output triangles
+struct Wrap_snap_visitor : public Autorefinement::Default_visitor
+{
+  template< typename Triangle>
+  void new_subdivision(Triangle& new_t, const Triangle& old_t) {
+    new_t.set_id(old_t.id());
+  }
+};
+
+/**
+*
+* rounds the coordinates of the points so that they fit in doubles while making and keeping the model intersection free by potentially subdividing the triangles.
+* The input can be any triangle soup and the output is an intersection-free triangle soup with Hausdorff distance
+* between the input and the output bounded by `M*2^-gs*k` where `M`is the maximum absolute coordinate in the model, `gs` the `snap_grid_size` (see below) and `k` the number of iterations
+* performed by the algorithm.
+*
+* @tparam PointRange a model of the concept `RandomAccessContainer`
+* whose value type is the point type
+* @tparam TriangleRange a model of the concepts `RandomAccessContainer`, `BackInsertionSequence` and `Swappable`, whose
+* value type is a model of the concept `RandomAccessContainer` whose value type is convertible to `std::size_t` and that
+* is constructible from an `std::initializer_list<std::size_t>` of size 3.
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* @param soup_points points of the soup of polygons
+* @param soup_triangles each element in the range describes a triangle using the indexed position of the points in `soup_points`
+* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{concurrency_tag}
+*     \cgalParamDescription{a tag indicating if the task should be done using one or several threads.}
+*     \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`}
+*     \cgalParamDefault{`CGAL::Sequential_tag`}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{point_map}
+*     \cgalParamDescription{a property map associating points to the elements of the range `soup_points`}
+*     \cgalParamType{a model of `ReadWritePropertyMap` whose value type is a point type}
+*     \cgalParamDefault{`CGAL::Identity_property_map`}
+*   \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 point type.}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{snap_grid_size}
+*     \cgalParamDescription{Scale the points to [-2^gs, 2^gs] where gs is the snap_grid_size before to round them on integer.}
+*     \cgalParamType{unsigned int}
+*     \cgalParamDefault{23}
+*     \cgalParamExtra{Must be lower than 52.}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{number_of_iterations}
+*     \cgalParamDescription{Maximum number of iteration performed by the algorithm.}
+*     \cgalParamType{unsigned int}
+*     \cgalParamDefault{5}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+* \return `true` if the modified triangle soup is free from self-intersection, and `false` if the algorithm was not
+* able to provide such a triangle soup within the number of iterations.
+*/
+template <typename PointRange, typename PolygonRange, class NamedParameters>
+bool polygon_soup_snap_rounding_impl(PointRange &points,
+                                     PolygonRange &triangles,
+                                     const NamedParameters& np)
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  // typedef typename GetPolygonSoupGeomTraits<PointRange, NamedParameters>::type GT;
+  typedef typename GetPointMap<PointRange, NamedParameters>::const_type    Point_map;
+  Point_map pm = choose_parameter<Point_map>(get_parameter(np, internal_np::point_map));
+
+  typedef typename internal_np::Lookup_named_param_def <
+    internal_np::concurrency_tag_t,
+    NamedParameters,
+    Sequential_tag
+  > ::type Concurrency_tag;
+
+  // visitor
+  typedef typename internal_np::Lookup_named_param_def <
+    internal_np::visitor_t,
+    NamedParameters,
+    Autorefinement::Default_visitor
+  >::type Visitor;
+  Visitor visitor(choose_parameter<Visitor>(get_parameter(np, internal_np::visitor)));
+
+  constexpr bool has_visitor = !std::is_same_v<Autorefinement::Default_visitor, Visitor>;
+  constexpr bool parallel_execution = std::is_same_v<Parallel_tag, Concurrency_tag>;
+  const size_t number_of_input_triangles = triangles.size();
+
+#ifndef CGAL_LINKED_WITH_TBB
+  static_assert (!parallel_execution,
+                 "Parallel_tag is enabled but TBB is unavailable.");
+#endif
+
+  using Point_3 = std::remove_cv_t<typename std::iterator_traits<typename PointRange::const_iterator>::value_type>;
+  using Kernel = typename Kernel_traits<Point_3>::Kernel;
+
+  // Get the grid size from the named parameter, the grid size could not be greater than 52
+  const unsigned int grid_size = (std::min)(52,choose_parameter(get_parameter(np, internal_np::snap_grid_size), 23));
+  const unsigned int max_nb_of_iteration = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 5);
+
+#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+  std::cout << "Compute the scaling of the coordinates" << std::endl;
+  std::cout << grid_size << std::endl;
+#endif
+
+  auto exp = [](const double v)
+  {
+    int n;
+    frexp(v, &n);
+    return n;
+  };
+  auto pow_2 = [](const int k)
+  {
+      return (k>=0)?std::pow(2,k):1./std::pow(2,-k);
+  };
+
+  // Get max absolute value of each dimension
+  Bbox_3 bb = bbox_3(points.begin(), points.end());
+  std::array<double, 3> max_abs{(std::max)(-bb.xmin(), bb.xmax()),
+                                (std::max)(-bb.ymin(), bb.ymax()),
+                                (std::max)(-bb.zmin(), bb.zmax())};
+  // Compute scale so that the exponent of max absolute value are 52-1.
+  std::array<double, 3> scale{pow_2(grid_size - exp(max_abs[0]) - 1),
+                              pow_2(grid_size - exp(max_abs[1]) - 1),
+                              pow_2(grid_size - exp(max_abs[2]) - 1)};
+
+#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+  std::cout << "Scaling: " << scale[0] << " " << scale[1] << " " << scale[2] << std::endl;
+#endif
+
+  auto snap = [](typename Kernel::FT x, double scale)
+  {
+    // Scale the coordinate, round to nearest integer and scale back
+    return double_ceil((x * scale) - 0.5) / scale;
+  };
+  auto snap_p = [scale, snap](const Point_3 &p)
+  {
+    return Point_3( snap(p.x(),scale[0]),
+                    snap(p.y(),scale[1]),
+                    snap(p.z(),scale[2]) );
+  };
+
+  for (std::size_t i = 0; i < max_nb_of_iteration; ++i)
+  {
+#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+    std::cout << "Start Iteration " << i << std::endl;
+    std::cout << "Model size: " << points.size() << " " << triangles.size() << std::endl;
+#endif
+
+#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+    std::cout << "Round all coordinates on doubles" << std::endl;
+#endif
+
+#ifdef CGAL_LINKED_WITH_TBB
+    if constexpr (parallel_execution)
+    {
+      tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()),
+                        [&](const tbb::blocked_range<size_t>& r){
+                          for(size_t pi = r.begin(); pi != r.end(); ++pi)
+                            points[pi] = Point_3(to_double(points[pi].x()), to_double(points[pi].y()), to_double(points[pi].z()));
+                        }
+                        );
+    } else
+#endif
+    for (Point_3 &p : points)
+      p = Point_3(to_double(p.x()), to_double(p.y()), to_double(p.z()));
+    repair_triangle_soup(points, triangles, np);
+
+    // Get all intersecting triangles
+    std::vector<std::pair<std::size_t, std::size_t>> pairs_of_intersecting_triangles;
+    triangle_soup_self_intersections<Concurrency_tag>(points, triangles, std::back_inserter(pairs_of_intersecting_triangles), np);
+
+    if (pairs_of_intersecting_triangles.empty())
+    {
+#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+      std::cout << "End of the snapping" << std::endl;
+#endif
+      if constexpr(has_visitor)
+      {
+        std::vector<std::vector<size_t> > map_io(number_of_input_triangles);
+        size_t id=0;
+        for(auto &t: triangles)
+          map_io[t.id()].push_back(id++);
+
+        visitor.number_of_output_triangles(triangles.size());
+        for(size_t src_id=0; src_id!=map_io.size(); ++src_id){
+          if(map_io[src_id].size()==0)
+            visitor.delete_triangle(src_id);
+          else if(map_io[src_id].size()==1 && !triangles[map_io[src_id][0]].was_modified())
+            visitor.verbatim_triangle_copy(map_io[src_id][0],src_id);
+          else
+            for(size_t new_id: map_io[src_id])
+              visitor.new_subtriangle(new_id,src_id);
+        }
+      }
+      return true;
+    }
+
+#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+    std::cout << "Number of pairs of intersecting triangles: " << pairs_of_intersecting_triangles.size() << std::endl;
+#endif
+
+#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+    std::cout << "Snap the coordinates of the vertices of the intersecting triangles" << std::endl;
+#endif
+
+// Variants of the rounding process depending of the macros defined
+#if defined(PMP_ROUNDING_VERTICES_NAIVE_SNAP)
+    // Nothing
+
+#elif defined(PMP_ROUNDING_VERTICES_CLOSEST_POINT_SNAP)
+    // Version where points in a voxel are rounded to the closest point
+
+    // Group the points of the vertices of the intersecting triangles by their voxel
+    std::map<Point_3, size_t> snap_points;
+    std::size_t index=0;
+    for (auto &pair : pairs_of_intersecting_triangles)
+    {
+      for (size_t pi : triangles[pair.first])
+      {
+        auto res=snap_points.emplace(snap_p(points[pi]), index);
+        if(res.second)
+          ++index;
+      }
+      for (size_t pi : triangles[pair.second])
+      {
+        auto res=snap_points.emplace(snap_p(points[pi]), index);
+        if(res.second)
+          ++index;
+      }
+    }
+
+    std::vector<std::vector<size_t>> identical_points(index);
+    for (size_t i=0; i!=points.size(); ++i)
+    {
+      Point_3 p_snap = snap_p(points[i]);
+      auto it=snap_points.find(p_snap);
+      if (it!=snap_points.end()){
+        identical_points[it->second].push_back(i);
+      }
+    }
+
+    // Replace all points in a voxel by the closest point
+    for(const auto &v: identical_points){
+      if(v.size()>1){
+        Point_3 center = snap_p(points[v[0]]);
+        size_t argmin(0);
+        for(size_t i=1; i!=v.size(); ++i){
+          if(Kernel().compare_distance_3_object()(center, points[v[i]], points[v[argmin]])==SMALLER)
+          {
+            argmin=i;
+          }
+        }
+        for(auto i: v){
+          points[i]=points[v[argmin]];
+        }
+      }
+    }
+
+#elif defined(PMP_ROUNDING_VERTICES_BARYCENTER_SNAP)
+    // Version where points in a voxel are rounded to their barycenter.
+
+    // Group the points of the vertices of the intersecting triangles by their voxel
+    std::map<Point_3, size_t> snap_points;
+    std::size_t index=0;
+    for (const auto &pair : pairs_of_intersecting_triangles)
+    {
+      for (size_t pi : triangles[pair.first])
+      {
+        auto res=snap_points.emplace(snap_p(points[pi]), index);
+        if(res.second)
+          ++index;
+      }
+      for (size_t pi : triangles[pair.second])
+      {
+        auto res=snap_points.emplace(snap_p(points[pi]), index);
+        if(res.second)
+          ++index;
+      }
+    }
+
+    std::vector<std::vector<size_t>> identical_points(index);
+    for (size_t i=0; i!=points.size(); ++i)
+    {
+      Point_3 p_snap = snap_p(points[i]);
+      auto it=snap_points.find(p_snap);
+      if (it!=snap_points.end()){
+        identical_points[it->second].push_back(i);
+      }
+    }
+
+    // Replace all points in a voxel by their barycenter
+    for(const auto &v: identical_points){
+      if(v.size()>1){
+        std::array<double, 3> a{0,0,0};
+        for(auto i: v){
+          a[0]+=to_double(points[i].x());
+          a[1]+=to_double(points[i].y());
+          a[2]+=to_double(points[i].z());
+        }
+        a[0]/=v.size();
+        a[1]/=v.size();
+        a[2]/=v.size();
+        for(auto i: v){
+          points[i]=Point_3(a[0],a[1],a[2]);
+        }
+      }
+    }
+#elif defined(PMP_ROUNDING_VERTICES_FLOAT_SNAP)
+    // Version where points are rounded to the closest simple precision float.
+
+    for (const auto &pair : pairs_of_intersecting_triangles)
+    {
+      for (size_t pi : triangles[pair.first])
+        points[pi]=Point_3( (float) to_double(points[pi].x()), (float) to_double(points[pi].y()), (float) to_double(points[pi].z()));
+      for (size_t pi : triangles[pair.second])
+        points[pi]=Point_3( (float) to_double(points[pi].x()), (float) to_double(points[pi].y()), (float) to_double(points[pi].z()));
+    }
+#else // Default Version
+    // Version where points are rounded to the center of their voxel.
+
+    // Get all the snap version of the points of the vertices of the intersecting triangles
+    // Note: points will not be modified here, they will be modified in the next for loop
+
+    std::vector<Point_3> snap_points;
+    snap_points.reserve(pairs_of_intersecting_triangles.size() * 3);
+
+    for (const auto &pair : pairs_of_intersecting_triangles)
+    {
+      for (size_t pi : triangles[pair.first])
+        snap_points.emplace_back(snap_p(points[pi]));
+      for (size_t pi : triangles[pair.second])
+        snap_points.emplace_back(snap_p(points[pi]));
+    }
+
+#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+    std::cout << "Snap the coordinates of the vertices close-by the previous ones" << std::endl;
+#endif
+
+    std::sort(snap_points.begin(), snap_points.end());
+    snap_points.erase(std::unique(snap_points.begin(), snap_points.end()), snap_points.end());
+
+    // If the snapped version of a point correspond to one of the previous point, we snap it
+#ifdef CGAL_LINKED_WITH_TBB
+    if constexpr(parallel_execution)
+    {
+      tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()),
+                        [&](const tbb::blocked_range<size_t>& r){
+                          for(size_t pi = r.begin(); pi != r.end(); ++pi){
+                            Point_3 p_snap=snap_p(points[pi]);
+                            if (std::binary_search(snap_points.begin(), snap_points.end(), p_snap))
+                              points[pi] = p_snap;
+                          }
+                        }
+                        );
+    } else
+#endif
+    for (Point_3 &p : points)
+    {
+      Point_3 p_snap = snap_p(p);
+      if (std::binary_search(snap_points.begin(), snap_points.end(), p_snap))
+        p = p_snap;
+    }
+#endif
+
+    repair_triangle_soup(points, triangles, np);
+#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+    std::cout << "Model size: " << points.size() << " " << triangles.size() << std::endl;
+    std::cout << "Autorefine the soup" << std::endl;
+#endif
+    if constexpr(has_visitor)
+    {
+      Wrap_snap_visitor visitor;
+      autorefine_triangle_soup(points, triangles, parameters::point_map(pm).concurrency_tag(Concurrency_tag()).visitor(visitor));
+    }
+    else
+    {
+      autorefine_triangle_soup(points, triangles, parameters::point_map(pm).concurrency_tag(Concurrency_tag()));
+    }
+  }
+  return false;
+}
+
+template <typename PointRange, typename PolygonRange, class NamedParameters>
+bool polygon_soup_snap_rounding(PointRange &soup_points,
+                                PolygonRange &soup_triangles,
+                                const NamedParameters& np)
+{
+  typedef typename internal_np::Lookup_named_param_def <
+    internal_np::visitor_t,
+    NamedParameters,
+    Autorefinement::Default_visitor//default
+  > ::type Visitor;
+
+  constexpr bool has_visitor = !std::is_same_v<Autorefinement::Default_visitor, Visitor>;
+  if constexpr(has_visitor)
+  {
+    //If a visitor is provided, we color the triangles with the index of their source to correctly track the modification
+    using Triangle = std::remove_cv_t<typename std::iterator_traits<typename PolygonRange::iterator>::value_type>;
+    std::vector<Indexes_range<Triangle> > indexes_soup_triangles;
+    size_t id=0;
+    for(typename PolygonRange::iterator it=soup_triangles.begin(); it!=soup_triangles.end(); ++it)
+      indexes_soup_triangles.emplace_back((*it), id++);
+
+    bool res=polygon_soup_snap_rounding_impl(soup_points, indexes_soup_triangles, np);
+
+    soup_triangles.clear();
+    soup_triangles.reserve(indexes_soup_triangles.size());
+    for(const Indexes_range<Triangle> &t: indexes_soup_triangles)
+      soup_triangles.push_back({t[0],t[1],t[2]});
+    return res;
+  }
+  else
+  {
+    return polygon_soup_snap_rounding_impl(soup_points, soup_triangles, np);
+  }
+}
+
+} } } //end of CGAL::Polygon_mesh_processing::autorefine_impl namespace
+
+#endif //CGAL_POLYGON_MESH_PROCESSING_POLYGON_SOUP_SNAP_ROUNDING_H

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -43,6 +43,7 @@ create_single_source_cgal_program("test_coref_epic_points_identity.cpp")
 create_single_source_cgal_program("test_does_bound_a_volume.cpp")
 create_single_source_cgal_program("test_pmp_clip.cpp")
 create_single_source_cgal_program("test_autorefinement.cpp")
+create_single_source_cgal_program("test_snap_rounding.cpp")
 create_single_source_cgal_program("autorefinement_sm.cpp")
 create_single_source_cgal_program( "corefine_non_manifold.cpp" )
 create_single_source_cgal_program("triangulate_hole_polyline_test.cpp")
@@ -115,6 +116,7 @@ if(TARGET CGAL::TBB_support)
   target_link_libraries(orient_polygon_soup_test PRIVATE CGAL::TBB_support)
   target_link_libraries(self_intersection_surface_mesh_test PRIVATE CGAL::TBB_support)
   target_link_libraries(test_autorefinement PRIVATE CGAL::TBB_support)
+  target_link_libraries(test_snap_rounding PRIVATE CGAL::TBB_support)
 else()
   message(STATUS "NOTICE: Intel TBB was not found. Tests will use sequential code.")
 endif()

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_1.off

@@ -0,0 +1,16 @@
+OFF
+8 4 0
+
+0 0 0
+0 1 0
+1 1 0
+1 0 0
+0 0 1
+0 1 1
+1 1 -1e-10
+1 0 -1e-10
+3  0 1 2
+3  0 2 3
+3  4 5 7
+3  5 6 7
+

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_2.off

@@ -0,0 +1,19 @@
+OFF
+12 4 0
+
+-1 -1 -1
+-1 -1 1
+1 -1e-15 0
+-1 1 -1
+-1 1 1
+1 1e-15 0
+0 -1e-15 -1
+0 -1e-15 1
+0 -2 0
+0 1e-15 -1
+0 1e-15 1
+0 2 0
+3 0 1 2
+3 3 4 5
+3 6 7 8
+3 9 10 11

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_3.off

@@ -0,0 +1,19 @@
+OFF
+12 4 0
+
+-1 -1 -1
+-1 -1 1
+1 -1e-15 0
+1 1 -1
+1 1 1
+1 1e-15 0
+0 -1e-15 -1
+0 -1e-15 1
+0 -2 0
+0 1e-15 -1
+0 1e-15 1
+0 2 0
+3 0 1 2
+3 3 4 5
+3 6 7 8
+3 9 10 11

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/coplanar_cubes.off

@@ -0,0 +1,202 @@
+OFF
+80 120 0
+-1 -5.627423701743102e-06 1.4142135623618981
+1 -5.627423701743102e-06 1.4142135623618981
+1 -1.4142135623618981 -5.627423701743102e-06
+-1 -1.4142135623618981 -5.627423701743102e-06
+-1 5.627423701743102e-06 -1.4142135623618981
+-1 1.4142135623618981 5.627423701743102e-06
+1 1.4142135623618981 5.627423701743102e-06
+1 5.627423701743102e-06 -1.4142135623618981
+-1 -4.9339086289293424e-06 1.4142135623644878
+1 -4.9339086289293424e-06 1.4142135623644878
+1 -1.4142135623644878 -4.9339086289293424e-06
+-1 -1.4142135623644878 -4.9339086289293424e-06
+-1 4.9339086289293424e-06 -1.4142135623644878
+-1 1.4142135623644878 4.9339086289293424e-06
+1 1.4142135623644878 4.9339086289293424e-06
+1 4.9339086289293424e-06 -1.4142135623644878
+-1 -4.3072127542901049e-06 1.4142135623665355
+1 -4.3072127542901049e-06 1.4142135623665355
+1 -1.4142135623665355 -4.3072127542901049e-06
+-1 -1.4142135623665355 -4.3072127542901049e-06
+-1 4.3072127542901049e-06 -1.4142135623665355
+-1 1.4142135623665355 4.3072127542901049e-06
+1 1.4142135623665355 4.3072127542901049e-06
+1 4.3072127542901049e-06 -1.4142135623665355
+-1 -3.7949909043547426e-06 1.4142135623680028
+1 -3.7949909043547426e-06 1.4142135623680028
+1 -1.4142135623680028 -3.7949909043547426e-06
+-1 -1.4142135623680028 -3.7949909043547426e-06
+-1 3.7949909043547426e-06 -1.4142135623680028
+-1 1.4142135623680028 3.7949909043547426e-06
+1 1.4142135623680028 3.7949909043547426e-06
+1 3.7949909043547426e-06 -1.4142135623680028
+-1 -3.6070768958898969e-06 1.4142135623684944
+1 -3.6070768958898969e-06 1.4142135623684944
+1 -1.4142135623684944 -3.6070768958898969e-06
+-1 -1.4142135623684944 -3.6070768958898969e-06
+-1 3.6070768958898969e-06 -1.4142135623684944
+-1 1.4142135623684944 3.6070768958898969e-06
+1 1.4142135623684944 3.6070768958898969e-06
+1 3.6070768958898969e-06 -1.4142135623684944
+-1 -2.4904507319126939e-06 1.4142135623709018
+1 -2.4904507319126939e-06 1.4142135623709018
+1 -1.4142135623709018 -2.4904507319126939e-06
+-1 -1.4142135623709018 -2.4904507319126939e-06
+-1 2.4904507319126939e-06 -1.4142135623709018
+-1 1.4142135623709018 2.4904507319126939e-06
+1 1.4142135623709018 2.4904507319126939e-06
+1 2.4904507319126939e-06 -1.4142135623709018
+-1 0.99999891029865429 1.0000010897001574
+1 0.99999891029865429 1.0000010897001574
+1 -1.0000010897001574 0.99999891029865429
+-1 -1.0000010897001574 0.99999891029865429
+-1 -0.99999891029865429 -1.0000010897001574
+-1 1.0000010897001574 -0.99999891029865429
+1 1.0000010897001574 -0.99999891029865429
+1 -0.99999891029865429 -1.0000010897001574
+-1 1.4142135623724554 1.3448878106118537e-06
+1 1.4142135623724554 1.3448878106118537e-06
+1 -1.3448878106118537e-06 1.4142135623724554
+-1 -1.3448878106118537e-06 1.4142135623724554
+-1 -1.4142135623724554 -1.3448878106118537e-06
+-1 1.3448878106118537e-06 -1.4142135623724554
+1 1.3448878106118537e-06 -1.4142135623724554
+1 -1.4142135623724554 -1.3448878106118537e-06
+-1 1.0000000341206472 -0.99999996587935136
+1 1.0000000341206472 -0.99999996587935136
+1 0.99999996587935136 1.0000000341206472
+-1 0.99999996587935136 1.0000000341206472
+-1 -1.0000000341206472 0.99999996587935136
+-1 -0.99999996587935136 -1.0000000341206472
+1 -0.99999996587935136 -1.0000000341206472
+1 -1.0000000341206472 0.99999996587935136
+-1 1.5942993642314858e-08 -1.4142135623730947
+1 1.5942993642314858e-08 -1.4142135623730947
+1 1.4142135623730947 1.5942993642314858e-08
+-1 1.4142135623730947 1.5942993642314858e-08
+-1 -1.5942993642314858e-08 1.4142135623730947
+-1 -1.4142135623730947 -1.5942993642314858e-08
+1 -1.4142135623730947 -1.5942993642314858e-08
+1 -1.5942993642314858e-08 1.4142135623730947
+3 6 7 4
+3 3 2 1
+3 6 1 2
+3 5 0 1
+3 4 3 0
+3 7 2 3
+3 1 0 3
+3 4 5 6
+3 2 7 6
+3 1 6 5
+3 0 5 4
+3 3 4 7
+3 14 15 12
+3 11 10 9
+3 14 9 10
+3 13 8 9
+3 12 11 8
+3 15 10 11
+3 9 8 11
+3 12 13 14
+3 10 15 14
+3 9 14 13
+3 8 13 12
+3 11 12 15
+3 22 23 20
+3 19 18 17
+3 22 17 18
+3 21 16 17
+3 20 19 16
+3 23 18 19
+3 17 16 19
+3 20 21 22
+3 18 23 22
+3 17 22 21
+3 16 21 20
+3 19 20 23
+3 30 31 28
+3 27 26 25
+3 30 25 26
+3 29 24 25
+3 28 27 24
+3 31 26 27
+3 25 24 27
+3 28 29 30
+3 26 31 30
+3 25 30 29
+3 24 29 28
+3 27 28 31
+3 38 39 36
+3 35 34 33
+3 38 33 34
+3 37 32 33
+3 36 35 32
+3 39 34 35
+3 33 32 35
+3 36 37 38
+3 34 39 38
+3 33 38 37
+3 32 37 36
+3 35 36 39
+3 46 47 44
+3 43 42 41
+3 46 41 42
+3 45 40 41
+3 44 43 40
+3 47 42 43
+3 41 40 43
+3 44 45 46
+3 42 47 46
+3 41 46 45
+3 40 45 44
+3 43 44 47
+3 54 55 52
+3 51 50 49
+3 54 49 50
+3 53 48 49
+3 52 51 48
+3 55 50 51
+3 49 48 51
+3 52 53 54
+3 50 55 54
+3 49 54 53
+3 48 53 52
+3 51 52 55
+3 62 63 60
+3 59 58 57
+3 62 57 58
+3 61 56 57
+3 60 59 56
+3 63 58 59
+3 57 56 59
+3 60 61 62
+3 58 63 62
+3 57 62 61
+3 56 61 60
+3 59 60 63
+3 70 71 68
+3 67 66 65
+3 70 65 66
+3 69 64 65
+3 68 67 64
+3 71 66 67
+3 65 64 67
+3 68 69 70
+3 66 71 70
+3 65 70 69
+3 64 69 68
+3 67 68 71
+3 78 79 76
+3 75 74 73
+3 78 73 74
+3 77 72 73
+3 76 75 72
+3 79 74 75
+3 73 72 75
+3 76 77 78
+3 74 79 78
+3 73 78 77
+3 72 77 76
+3 75 76 79

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/cubes.off

@@ -0,0 +1,204 @@
+OFF
+80 120 0
+
+-0.61594615543034026 -1.0701741879286757 -1.2146347356723264
+1.0452131429007676 -1.3737111650653515 -0.14299342950741956
+1.5883443594490732 0.52704573278436739 -0.44652546552251943
+-0.072814938882034802 0.83058270992104299 -1.5181667716874265
+-1.0452131429007676 1.3737111650653515 0.14299342950741956
+-1.5883443594490732 -0.52704573278436739 0.44652546552251943
+0.072814938882034802 -0.83058270992104299 1.5181667716874265
+0.61594615543034026 1.0701741879286757 1.2146347356723264
+-0.61594639934604423 -1.0701743934138028 -1.2146344309355146
+1.0452122820086462 -1.3737119347288513 -0.14299232819208521
+1.5883448836888012 0.52704444686667884 -0.44652511853937799
+-0.072813797665889982 0.83058198818172735 -1.518167221282807
+-1.0452122820086462 1.3737119347288513 0.14299232819208521
+-1.5883448836888012 -0.52704444686667884 0.44652511853937799
+0.072813797665889982 -0.83058198818172735 1.518167221282807
+0.61594639934604423 1.0701743934138028 1.2146344309355146
+-0.6159463762173395 -1.07017459735539 -1.2146342629779363
+1.0452117929893034 -1.3737123992443807 -0.14299144015755649
+1.5883452549553072 0.52704365666978259 -0.44652473058622827
+-0.072812914251335584 0.83058145855877352 -1.518167553406609
+-1.0452117929893034 1.3737123992443807 0.14299144015755649
+-1.5883452549553072 -0.52704365666978259 0.44652473058622827
+0.072812914251335584 -0.83058145855877352 1.518167553406609
+0.6159463762173395 1.07017459735539 1.2146342629779363
+-0.61594668912789197 -1.0701744000985252 -1.2146342780961472
+1.0452104492672072 -1.3737135467380566 -0.14299023830665158
+1.5883461040242048 0.52704182836364122 -0.44652386832759344
+-0.072811034370893465 0.83058097500317263 -1.5181679081170896
+-1.0452104492672072 1.3737135467380566 0.14299023830665158
+-1.5883461040242048 -0.52704182836364122 0.44652386832759344
+0.072811034370893465 -0.83058097500317263 1.5181679081170896
+0.61594668912789197 1.0701744000985252 1.2146342780961472
+-0.6159468868925585 -1.070174149353629 -1.2146343987318162
+1.0452090174071582 -1.3737147790282049 -0.14298886603597585
+1.5883470704382798 0.52703987181399814 -0.44652274001266495
+-0.072808833861436925 0.83058050148857432 -1.5181682727085062
+-1.0452090174071582 1.3737147790282049 0.14298886603597585
+-1.5883470704382798 -0.52703987181399814 0.44652274001266495
+0.072808833861436925 -0.83058050148857432 1.5181682727085062
+0.6159468868925585 1.070174149353629 1.2146343987318162
+-0.61594620872112704 -1.0701743005704609 -1.2146346094034179
+1.0452091168725113 -1.3737147792450597 -0.14298813688653067
+1.5883473052656241 0.52703984477504562 -0.44652193661109801
+-0.072808020328014339 0.83058032344964505 -1.5181684091279861
+-1.0452091168725113 1.3737147792450597 0.14298813688653067
+-1.5883473052656241 -0.52703984477504562 0.44652193661109801
+0.072808020328014339 -0.83058032344964505 1.5181684091279861
+0.61594620872112704 1.0701743005704609 1.2146346094034179
+0.015812704102963069 -1.4305657319567158 -0.97633582590218304
+-0.063147971057831695 -1.3970994075162435 1.0218246323460018
+1.3734324159237437 -0.0060147457399582988 1.0552948506097966
+1.4523930910845384 -0.03948107018043074 -0.94286560763838845
+0.063147971057831695 1.3970994075162435 -1.0218246323460018
+-1.3734324159237437 0.0060147457399582988 -1.0552948506097966
+-1.4523930910845384 0.03948107018043074 0.94286560763838845
+-0.015812704102963069 1.4305657319567158 0.97633582590218304
+-0.017005997503037738 -1.7115480160220284 -0.26516822000379986
+-1.4526640017493855 -0.672335677974221 0.66161320584190753
+-0.062496589720092359 0.32138030895690106 1.7008259385628288
+1.3731614145262554 -0.71783202909090671 0.77404451271712071
+1.4526640017493855 0.672335677974221 -0.66161320584190753
+0.062496589720092359 -0.32138030895690106 -1.7008259385628288
+-1.3731614145262554 0.71783202909090671 -0.77404451271712071
+0.017005997503037738 1.7115480160220284 0.26516822000379986
+-0.67677662809399264 -1.5909903432533568 0.10355251489119292
+-1.5303314706074931 -0.030328946305495721 -0.81065760039877899
+-1.0732219475484153 1.1338825933415291 0.75000327720449067
+-0.21966710503491477 -0.42677880360633264 1.6642133924944624
+1.5303314706074931 0.030328946305495721 0.81065760039877899
+1.0732219475484153 -1.1338825933415291 -0.75000327720449067
+0.21966710503491477 0.42677880360633264 -1.6642133924944624
+0.67677662809399264 1.5909903432533568 -0.10355251489119292
+-1.207106772440474 -1.2071067303442111 -0.29289346439627145
+-0.20710789493838763 -0.20710532012834212 -1.7071068233208859
+-0.50000073181889593 1.5000005345367655 -0.70710512978621964
+-1.4999996093209822 0.4999991243208966 0.70710822913839488
+0.20710789493838763 0.20710532012834212 1.7071068233208859
+0.50000073181889593 -1.5000005345367655 0.70710512978621964
+1.4999996093209822 -0.4999991243208966 -0.70710822913839488
+1.207106772440474 1.2071067303442111 0.29289346439627145
+3  6 7 4
+3  3 2 1
+3  6 1 2
+3  5 0 1
+3  4 3 0
+3  7 2 3
+3  1 0 3
+3  4 5 6
+3  2 7 6
+3  1 6 5
+3  0 5 4
+3  3 4 7
+3  14 15 12
+3  11 10 9
+3  14 9 10
+3  13 8 9
+3  12 11 8
+3  15 10 11
+3  9 8 11
+3  12 13 14
+3  10 15 14
+3  9 14 13
+3  8 13 12
+3  11 12 15
+3  22 23 20
+3  19 18 17
+3  22 17 18
+3  21 16 17
+3  20 19 16
+3  23 18 19
+3  17 16 19
+3  20 21 22
+3  18 23 22
+3  17 22 21
+3  16 21 20
+3  19 20 23
+3  30 31 28
+3  27 26 25
+3  30 25 26
+3  29 24 25
+3  28 27 24
+3  31 26 27
+3  25 24 27
+3  28 29 30
+3  26 31 30
+3  25 30 29
+3  24 29 28
+3  27 28 31
+3  38 39 36
+3  35 34 33
+3  38 33 34
+3  37 32 33
+3  36 35 32
+3  39 34 35
+3  33 32 35
+3  36 37 38
+3  34 39 38
+3  33 38 37
+3  32 37 36
+3  35 36 39
+3  46 47 44
+3  43 42 41
+3  46 41 42
+3  45 40 41
+3  44 43 40
+3  47 42 43
+3  41 40 43
+3  44 45 46
+3  42 47 46
+3  41 46 45
+3  40 45 44
+3  43 44 47
+3  54 55 52
+3  51 50 49
+3  54 49 50
+3  53 48 49
+3  52 51 48
+3  55 50 51
+3  49 48 51
+3  52 53 54
+3  50 55 54
+3  49 54 53
+3  48 53 52
+3  51 52 55
+3  62 63 60
+3  59 58 57
+3  62 57 58
+3  61 56 57
+3  60 59 56
+3  63 58 59
+3  57 56 59
+3  60 61 62
+3  58 63 62
+3  57 62 61
+3  56 61 60
+3  59 60 63
+3  70 71 68
+3  67 66 65
+3  70 65 66
+3  69 64 65
+3  68 67 64
+3  71 66 67
+3  65 64 67
+3  68 69 70
+3  66 71 70
+3  65 70 69
+3  64 69 68
+3  67 68 71
+3  78 79 76
+3  75 74 73
+3  78 73 74
+3  77 72 73
+3  76 75 72
+3  79 74 75
+3  73 72 75
+3  76 77 78
+3  74 79 78
+3  73 78 77
+3  72 77 76
+3  75 76 79
+

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_1.off

@@ -0,0 +1,14 @@
+OFF
+7 3 0
+
+0 0 0
+1 2 0
+1 0 0
+0 1 0
+1e-30 0 -8
+1 2 16
+2 2 2
+3  0 1 3
+3  0 2 3
+3  4 5 6
+

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_2.off

@@ -0,0 +1,15 @@
+OFF
+9 3 0
+
+0 1 0
+0 0 0
+2 1 0
+1 1 0
+2 0 0
+1.00000000000002 0 0
+1.00000000000001 1e-30 0
+2 -1 0
+1.00000000000002 -1 0
+3 0 1 2
+3 3 4 5
+3 6 7 8

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_3.off

@@ -0,0 +1,19 @@
+OFF
+12 4 0
+
+0 0 0
+0 8.000000000001 0
+8.000000000001 0 0
+4.000000000002 4.000000000002 0
+4 8.000000000001 0
+8.000000000001 4 0
+6.000000000002 6.000000000002 0
+6 8.000000000001 0
+8.000000000001 6 0
+6.5 6.5 0
+6.5 8.000000000001 0
+8.000000000001 6.5 0
+3 0 1 2
+3 3 4 5
+3 6 7 8
+3 9 10 11

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/no_collapse_1.off

@@ -0,0 +1,16 @@
+OFF
+8 4 0
+
+0 0 0
+0 1 0
+1 1 0
+1 0 0
+0 0 1
+0 1 1
+1 1 1e-10
+1 0 1e-10
+3  0 1 2
+3  0 2 3
+3  4 5 7
+3  5 6 7
+

Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/no_collapse_2.off

@@ -0,0 +1,19 @@
+OFF
+12 4 0
+
+1 -1 -1
+1 -1 1
+3 -1e-15 0
+1 1 -1
+1 1 1
+3 1e-15 0
+0 -1e-15 -1
+0 -1e-15 1
+0 -2 0
+0 1e-15 -1
+0 1e-15 1
+0 2 0
+3 0 1 2
+3 3 4 5
+3 6 7 8
+3 9 10 11

Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cpp

@@ -30,7 +30,7 @@ struct My_exp_visitor :
   std::shared_ptr<int> i;
 };
 
-struct My_visitor
+struct My_visitor : public PMP::Autorefinement::Default_visitor
 {
   My_visitor(std::size_t nb_input, std::size_t expected_nb_output)
     : nb_input(nb_input)
@@ -66,6 +66,11 @@ struct My_visitor
     tgt_check[tgt_id]=1;
   }
 
+  void delete_triangle(std::size_t src_id)
+  {
+    assert(src_id<nb_input);
+  }
+
   std::size_t nb_input;
   std::size_t expected_nb_output;
   std::vector<int> tgt_check;

Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding.cmd

@@ -0,0 +1,37 @@
+data-autoref/test_01.off 4 2
+data-autoref/test_02.off 10 17
+data-autoref/test_03.off 13 20
+data-autoref/test_04.off 12 20
+data-autoref/test_05.off 12 20
+data-autoref/test_06.off 89 248
+data-autoref/test_07.off 8 12
+data-autoref/test_08.off 57 148
+data-autoref/test_09.off 4 3
+data-autoref/test_10.off 10 13
+data-autoref/test_11.off 9 12
+data-autoref/test_12.off 9 6
+data-autoref/test_13.off 12 22
+data-autoref/test_14.off 12 22
+data-autoref/test_15.off 12 22
+data-autoref/test_16.off 4 2
+data-autoref/test_17.off 4 2
+data-autoref/triple_inter_exception/triple.off 15 18
+data-autoref/triple_inter_exception/cubes_cpln_1.off 66 206
+data-autoref/triple_inter_exception/cubes_cpln_2.off 54 158
+data-autoref/triple_inter_exception/cubes_cpln_3.off 61 188
+data-autoref/open_01.off 1313 2622
+data-autoref/open_02.off 562 1056
+data-autoref/cpln_01.off 30 78
+data-autoref/cpln_02.off 24 56
+data-autoref/cpln_03.off 27 68
+data-autoref/four_cubes.off 106 352
+data-autoref/spiral.off 19 31
+data-snap/intersection_when_rounded_1.off 9 6
+data-snap/intersection_when_rounded_2.off 10 7
+data-snap/intersection_when_rounded_3.off 16 14
+data-snap/collapse_1.off 6 4
+data-snap/collapse_2.off 13 16
+data-snap/collapse_3.off 11 10
+data-snap/no_collapse_1.off 8 4
+data-snap/no_collapse_2.off 12 4
+data-snap/coplanar_cubes.off 1514 5544

Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding.cpp

@@ -0,0 +1,100 @@
+
+#include <CGAL/Polygon_mesh_processing/intersection.h>
+#include <CGAL/Polygon_mesh_processing/corefinement.h>
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+
+#include <CGAL/IO/polygon_soup_io.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <fstream>
+#include <sstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+struct My_visitor : public PMP::Autorefinement::Default_visitor
+{
+  My_visitor(std::size_t nb_input, std::size_t expected_nb_output)
+    : nb_input(nb_input)
+    , expected_nb_output(expected_nb_output)
+  {}
+
+  ~My_visitor()
+  {
+    for(std::size_t i=0; i<tgt_check.size(); ++i)
+    {
+      assert( tgt_check[i]==1 );
+    }
+  }
+  void number_of_output_triangles(std::size_t nbt)
+  {
+    tgt_check.assign(expected_nb_output, 0);
+    assert(nbt==expected_nb_output);
+  }
+
+  void verbatim_triangle_copy(std::size_t tgt_id, std::size_t src_id)
+  {
+    assert(src_id<nb_input);
+    assert(tgt_id<expected_nb_output);
+    assert(tgt_check.size()==expected_nb_output && tgt_check[tgt_id]==0);
+    tgt_check[tgt_id]=1;
+  }
+
+  void new_subtriangle(std::size_t tgt_id, std::size_t src_id)
+  {
+    assert(src_id<nb_input);
+    assert(tgt_id<expected_nb_output);
+    assert(tgt_check.size()==expected_nb_output && tgt_check[tgt_id]==0);
+    tgt_check[tgt_id]=1;
+  }
+
+  void delete_triangle(std::size_t src_id)
+  {
+    assert(src_id<nb_input);
+  }
+
+  std::size_t nb_input;
+  std::size_t expected_nb_output;
+  std::vector<int> tgt_check;
+};
+
+template <class Tag>
+void test(const char* fname, std::size_t nb_vertices_after_autorefine, std::size_t expected_nb_output, Tag tag)
+{
+  std::cout << "Running tests on " << fname;
+  if (std::is_same_v<Tag, CGAL::Sequential_tag>)
+    std::cout << " (Sequential)\n";
+  else
+    std::cout << " (Parallel)\n";
+
+  std::vector<K::Point_3> points;
+  std::vector< std::vector<std::size_t> > triangles;
+  if (!CGAL::IO::read_polygon_soup(fname, points, triangles))
+  {
+    std::cerr << "  Input mesh is not a valid file." << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+// Testing autorefine()
+  My_visitor visitor(triangles.size(), expected_nb_output);
+  PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::visitor(visitor).concurrency_tag(tag).apply_iterative_snap_rounding(true));
+  assert( nb_vertices_after_autorefine==points.size());
+  assert( expected_nb_output==triangles.size());
+  assert( !PMP::does_triangle_soup_self_intersect(points, triangles) );
+//  CGAL::IO::write_polygon_soup("/tmp/debug.off", points, triangles);
+}
+
+int main(int argc, const char** argv)
+{
+  // file expected_nb_of_vertices expected_nb_of_faces (after repair)
+  for (int i=0;i<(argc-1)/3; ++i)
+  {
+    test(argv[1+3*i], atoi(argv[1+3*i+1]), atoi(argv[1+3*i+2]), CGAL::Sequential_tag());
+#ifdef CGAL_LINKED_WITH_TBB
+    test(argv[1+3*i], atoi(argv[1+3*i+1]), atoi(argv[1+3*i+2]), CGAL::Parallel_tag());
+#endif
+  }
+}

Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding_full.cmd

@@ -0,0 +1,38 @@
+data-autoref/test_01.off 4 2
+data-autoref/test_02.off 10 17
+data-autoref/test_03.off 13 20
+data-autoref/test_04.off 12 20
+data-autoref/test_05.off 12 20
+data-autoref/test_06.off 89 248
+data-autoref/test_07.off 8 12
+data-autoref/test_08.off 57 148
+data-autoref/test_09.off 4 3
+data-autoref/test_10.off 10 13
+data-autoref/test_11.off 9 12
+data-autoref/test_12.off 9 6
+data-autoref/test_13.off 12 22
+data-autoref/test_14.off 12 22
+data-autoref/test_15.off 12 22
+data-autoref/test_16.off 4 2
+data-autoref/test_17.off 4 2
+data-autoref/triple_inter_exception/triple.off 15 18
+data-autoref/triple_inter_exception/cubes_cpln_1.off 66 206
+data-autoref/triple_inter_exception/cubes_cpln_2.off 54 158
+data-autoref/triple_inter_exception/cubes_cpln_3.off 61 188
+data-autoref/open_01.off 1313 2622
+data-autoref/open_02.off 562 1056
+data-autoref/cpln_01.off 30 78
+data-autoref/cpln_02.off 24 56
+data-autoref/cpln_03.off 27 68
+data-autoref/four_cubes.off 106 352
+data-autoref/spiral.off 19 31
+data-snap/intersection_when_rounded_1.off 9 6
+data-snap/intersection_when_rounded_2.off 10 7
+data-snap/intersection_when_rounded_3.off 16 14
+data-snap/collapse_1.off 6 4
+data-snap/collapse_2.off 13 16
+data-snap/collapse_3.off 11 10
+data-snap/no_collapse_1.off 8 4
+data-snap/no_collapse_2.off 12 4
+data-snap/cubes.off 5042 26036
+data-snap/coplanar_cubes.off 1514 5544

STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h

@@ -177,7 +177,8 @@ CGAL_add_named_parameter(weight_limiting_t, weight_limiting, weight_limiting)
 CGAL_add_named_parameter(progressive_t, progressive, progressive)
 CGAL_add_named_parameter(tiling_t, tiling, tiling)
 CGAL_add_named_parameter(dimension_t, dimension, dimension)
-
+CGAL_add_named_parameter(apply_iterative_snap_rounding_t, apply_iterative_snap_rounding, apply_iterative_snap_rounding)
+CGAL_add_named_parameter(snap_grid_size_t, snap_grid_size, snap_grid_size)
 
 // List of named parameters that we use in the package 'Surface Mesh Simplification'
 CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost)