Add do_snap parameter to PMP::autorefine_triangle_soup

Documentation/doc/biblio/geom.bib

@@ -152085,3 +152085,16 @@ @article{ledoux2014triangulation
   year={2014},
   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},
+}

Installation/CHANGES.md

@@ -9,7 +9,7 @@
 ### [Polygon Mesh Processing](https://doc.cgal.org/6.1/Manual/packages.html#PkgPolygonMeshProcessing)
 - Added the function `CGAL::Polygon_mesh_processing::discrete_mean_curvature` and `CGAL::Polygon_mesh_processing::discrete_Guassian_curvature` to evaluate the discrete curvature at a vertex of a mesh.
 - Added the function `CGAL::Polygon_mesh_processing::angle_sum` to compute the sum of the angles around a vertex.
-
+- 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 remains free of self-intersections.
 
 ### [Algebraic Kernel](https://doc.cgal.org/6.1/Manual/packages.html#PkgAlgebraicKernelD)
 

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\"}"