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},
+}

Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h

@@ -23,5 +23,8 @@ class PMPAutorefinementVisitor{
   /// 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 triangle delete because it was or became degenerated, `src_id` was its position before calling `autorefine_triangle_soup()`.
+  /// A triangle may degenerate inside `autorefine_triangle_soup()` if and only if `apply_iterative_snap_rounding` is set to `true`.
+  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, and Ilker %O. Yaz
+\authors David Coeurjolly, Jaques-Olivier Lachaud, Sylvain Lazard, Konstantinos Katrioplas, Sébastien Loriot, Ivan Pađen, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, Léo Valque and Ilker %O. Yaz
 
 \image html neptun_head.jpg
 \image latex neptun_head.jpg
@@ -889,8 +889,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.
-To guarantee that the triangle soup is free from self-intersections, a kernel with exact constructions must be used.
+might be introduced due to rounding issues of points coordinates. The `apply_iterative_snap_rounding` option can be used to resolve this issue.
+When set to `true`, it ensures 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
@@ -1450,5 +1452,7 @@ used as a reference during the project.
 
 The curvature-based sizing field version of isotropic remeshing was added by Ivan Pađen during GSoC 2023, under the supervision of Sébastien Loriot and Jane Tournois.
 
+The `apply_iterative_snap_rounding` option for autorefinement were implemented during 2025. This was implemented by Léo Valque and Sylvain Lazard. The implementation is based on \cgalCite{lazard:hal-04907149}.
+
 */
 } /* namespace CGAL */

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)
@@ -138,6 +139,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,125 @@
+#define PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE
+// #define CGAL_PMP_AUTOREFINE_USE_DEFAULT_VERBOSE
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/IO/polygon_soup_io.h>
+#include <CGAL/Real_timer.h>
+
+#include <boost/container/small_vector.hpp>
+
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK;
+typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK;
+typedef CGAL::Simple_cartesian<double> Cartesian;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+struct Track_visitor
+{
+  // inline void number_of_output_triangles(std::size_t nbt) { std::cout << "total number of triangles: " << nbt << std::endl;}
+  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) { std::cout << src_id << " goes to " << tgt_id << std::endl;}
+  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) { std::cout << src_id << " subdivides to " << tgt_id << std::endl;}
+  inline void new_subtriangle(std::size_t /*tgt_id*/, std::size_t /*src_id*/) {}
+
+  // inline void delete_triangle(std::size_t src_id) { std::cout << src_id << " deleted" << std::endl;}
+  inline void delete_triangle(std::size_t /*src_id*/) {}
+};
+
+int main(int argc, char** argv)
+{
+  const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
+                                         : std::string(argv[1]);
+
+  // const std::string out_file = argc <= 2 ? "rounded_soup.off"
+  //                                        : std::string(argv[2]);
+  const std::string out_file = "rounded_soup.off";
+
+  const int grid_size = argc <= 2 ? 23
+                                  : std::stoi(std::string(argv[2]));
+
+  const bool epeck = argc <= 3 ? false
+                               : std::string(argv[3])=="EPECK";
+
+  if(epeck)
+  {
+    std::vector<EPECK::Point_3> input_points;
+    std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
+
+    std::cout << "Snap rounding apply on " << filename << " with EPECK\n";
+    if (!CGAL::IO::read_polygon_soup(filename, input_points, input_triangles))
+    {
+      std::cerr << "Cannot read " << filename << "\n";
+      return 1;
+    }
+    std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << std::endl;
+    std::cout << "Is 2-manifold: " << PMP::is_polygon_soup_a_polygon_mesh(input_triangles) << std::endl;
+
+    PMP::repair_polygon_soup(input_points, input_triangles);
+    PMP::triangulate_polygons(input_points, input_triangles);
+
+    CGAL::Real_timer t;
+    t.start();
+    PMP::autorefine_triangle_soup(input_points, input_triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size));
+    t.stop();
+    std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << " in " << t.time() << " sec." << std::endl;
+    std::cout << "Does self-intersect: " << PMP::does_triangle_soup_self_intersect(input_points, input_triangles) << std::endl;
+    std::cout << "Is 2-manifold: " << PMP::orient_polygon_soup(input_points, input_triangles) << "\n\n"  << std::endl;
+
+    std::vector<Cartesian::Point_3> output_points;
+    for(auto &p: input_points)
+      output_points.emplace_back(CGAL::to_double(p.x()),CGAL::to_double(p.y()),CGAL::to_double(p.z()));
+
+    CGAL::IO::write_polygon_soup(out_file, output_points, input_triangles, CGAL::parameters::stream_precision(17));
+  }
+  else
+  {
+    std::vector<EPICK::Point_3> input_points;
+    std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
+    std::cout << "Snap rounding on " << filename << " with EPICK\n";
+    if (!CGAL::IO::read_polygon_soup(filename, input_points, input_triangles))
+    {
+      std::cerr << "Cannot read " << filename << "\n";
+      return 1;
+    }
+    std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << std::endl;
+    std::cout << "Is 2-manifold: " << PMP::is_polygon_soup_a_polygon_mesh(input_triangles) << std::endl;
+
+    std::vector<std::pair<std::size_t, std::size_t>> pairs_of_intersecting_triangles;
+    PMP::triangle_soup_self_intersections(input_points, input_triangles, std::back_inserter(pairs_of_intersecting_triangles));
+    std::cout << "Nb of pairs of intersecting triangles: " << pairs_of_intersecting_triangles.size() << std::endl;
+
+    PMP::repair_polygon_soup(input_points, input_triangles);
+    PMP::triangulate_polygons(input_points, input_triangles);
+
+    CGAL::Real_timer t;
+    t.start();
+#if 1
+    Track_visitor visitor;
+    bool success=PMP::autorefine_triangle_soup(input_points, input_triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).visitor(visitor));
+#else
+    bool success=PMP::autorefine_triangle_soup(input_points, input_triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).number_of_iterations(15));
+#endif
+    t.stop();
+    std::cout << "\nOutput:" << std::endl;
+    std::cout << "#points = " << input_points.size() << " and #triangles = " << input_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>(input_points, input_triangles) << std::endl;
+    else
+      std::cout << "ROUNDING FAILED" << std::endl;
+    CGAL::IO::write_polygon_soup(out_file, input_points, input_triangles, CGAL::parameters::stream_precision(17));
+    std::cout << "Is 2-manifold: " << PMP::orient_polygon_soup(input_points, input_triangles) << "\n\n"  << std::endl;
+  }
+  return 0;
+}