Tetrahedral_remeshing - make it determistic (#8945)

## Summary of Changes

Tetrahedral_remeshing is not deterministic (though outputs are very
close in 2 consecutive runs).
Making it deterministic...

## Release Management

* Affected package(s): Tetrahedral_remeshing
* License and copyright ownership: unchanged

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h

@@ -91,6 +91,13 @@ private:
   std::vector<bool> m_free_vertices{};
   bool m_flip_smooth_steps{false};
 
+  struct Move
+  {
+    Vector_3 move;
+    int neighbors;
+    FT mass;
+  };
+
 public:
   Tetrahedral_remeshing_smoother(const SizingFunction& sizing,
                                  const CellSelector& cell_selector,
@@ -627,8 +634,8 @@ private:
                            const C3t3& c3t3,
                            const bool boundary_edge = false) const
   {
-    const auto mwi = midpoint_with_info(e, boundary_edge, c3t3);
-    const FT s = sizing_at_midpoint(e, mwi.dim, mwi.index, m_sizing, c3t3, m_cell_selector);
+    const auto [pt, dim, index] = midpoint_with_info(e, boundary_edge, c3t3);
+    const FT s = sizing_at_midpoint(e, pt, dim, index, m_sizing, c3t3, m_cell_selector);
     const FT density = 1. / s; //density = 1 / size^(dimension)
                  //edge dimension is 1, so density = 1 / size
                  //to have mass = length * density with no dimension
@@ -655,9 +662,8 @@ private:
     auto& tr = c3t3.triangulation();
 
     const std::size_t nbv = tr.number_of_vertices();
-    std::vector<Vector_3> moves(nbv, CGAL::NULL_VECTOR);
-    std::vector<int> neighbors(nbv, 0);
-    std::vector<FT> masses(nbv, 0.);
+    const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/};
+    std::vector<Move> moves(nbv, default_move);
 
     //collect neighbors
     for (const Edge& e : c3t3.edges_in_complex())
@@ -683,33 +689,35 @@ private:
 
       if (vh0_moving)
       {
-        moves[i0] += density * Vector_3(p0, p1);
-        neighbors[i0]++;
-        masses[i0] += density;
+        moves[i0].move += density * Vector_3(p0, p1);
+        moves[i0].mass += density;
+        ++moves[i0].neighbors;
       }
       if (vh1_moving)
       {
-        moves[i1] += density * Vector_3(p1, p0);
-        neighbors[i1]++;
-        masses[i1] += density;
+        moves[i1].move += density * Vector_3(p1, p0);
+        moves[i1].mass += density;
+        ++moves[i1].neighbors;
       }
     }
 
-    // iterate over map of <vertex, id>
-    for (auto [v, vid] : m_vertex_id)
+    // iterate over vertices and move
+    for(Vertex_handle v : tr.finite_vertex_handles())
     {
+      const std::size_t vid = vertex_id(v);
+
       if (!is_free(vid) || !is_on_feature(v))
         continue;
 
       const Point_3 current_pos = point(v->point());
 
-      const std::size_t nb_neighbors = neighbors[vid];
+      const std::size_t nb_neighbors = moves[vid].neighbors;
       if(nb_neighbors == 0)
         continue;
 
-      CGAL_assertion(masses[vid] > 0);
+      CGAL_assertion(moves[vid].mass > 0);
       const Vector_3 move = (nb_neighbors > 0)
-                          ? moves[vid] / masses[vid]
+                          ? moves[vid].move / moves[vid].mass
                           : CGAL::NULL_VECTOR;
 
       const Point_3 smoothed_position = current_pos + move;
@@ -771,9 +779,8 @@ std::size_t smooth_vertices_on_surfaces(C3t3& c3t3,
   auto& tr = c3t3.triangulation();
 
   const std::size_t nbv = tr.number_of_vertices();
-  std::vector<Vector_3> moves(nbv, CGAL::NULL_VECTOR);
-  std::vector<int> neighbors(nbv, 0);
-  std::vector<FT> masses(nbv, 0.);
+  const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/};
+  std::vector<Move> moves(nbv, default_move);
 
   for (const Edge& e : tr.finite_edges())
   {
@@ -797,26 +804,28 @@ std::size_t smooth_vertices_on_surfaces(C3t3& c3t3,
 
       if (vh0_moving)
       {
-        moves[i0] += density * Vector_3(p0, p1);
-        neighbors[i0]++;
-        masses[i0] += density;
+        moves[i0].move += density * Vector_3(p0, p1);
+        moves[i0].mass += density;
+        ++moves[i0].neighbors;
       }
       if (vh1_moving)
       {
-        moves[i1] += density * Vector_3(p1, p0);
-        neighbors[i1]++;
-        masses[i1] += density;
+        moves[i1].move += density * Vector_3(p1, p0);
+        moves[i1].mass += density;
+        ++moves[i1].neighbors;
       }
     }
   }
 
-  // iterate over map of <vertex, id>
-  for (auto [v, vid] : m_vertex_id)
+  // iterate over vertices and move
+  for(Vertex_handle v : tr.finite_vertex_handles())
   {
+    const std::size_t vid = vertex_id(v);
+
     if (!is_free(vid) || v->in_dimension() != 2)
       continue;
 
-    const std::size_t nb_neighbors = neighbors[vid];
+    const std::size_t nb_neighbors = moves[vid].neighbors;
     const Point_3 current_pos = point(v->point());
 
     const auto& incident_surface_patches = vertices_surface_indices.at(v);
@@ -830,7 +839,7 @@ std::size_t smooth_vertices_on_surfaces(C3t3& c3t3,
 
     if (nb_neighbors > 1)
     {
-      const Vector_3 move = moves[vid] / masses[vid];
+      const Vector_3 move = moves[vid].move / moves[vid].mass;
       const Point_3 smoothed_position = point(v->point()) + move;
 
 #ifdef CGAL_TET_REMESHING_SMOOTHING_WITH_MLS
@@ -969,9 +978,9 @@ std::size_t smooth_internal_vertices(C3t3& c3t3,
   auto& tr = c3t3.triangulation();
 
   const std::size_t nbv = tr.number_of_vertices();
-  std::vector<Vector_3> moves(nbv, CGAL::NULL_VECTOR);
-  std::vector<int> neighbors(nbv, 0);/*for dim 3 vertices, start counting directly from 0*/
-  std::vector<FT> masses(nbv, 0.);
+  const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/};
+  std::vector<Move> moves(nbv, default_move);
+  /*for dim 3 vertices, start counting neighbors directly from 0*/
 
   for (const Edge& e : tr.finite_edges())
   {
@@ -979,8 +988,7 @@ std::size_t smooth_internal_vertices(C3t3& c3t3,
       continue;
     else
     {
-      const Vertex_handle vh0 = e.first->vertex(e.second);
-      const Vertex_handle vh1 = e.first->vertex(e.third);
+      const auto [vh0, vh1] =  make_vertex_pair(e);
 
       const std::size_t& i0 = vertex_id(vh0);
       const std::size_t& i1 = vertex_id(vh1);
@@ -997,31 +1005,32 @@ std::size_t smooth_internal_vertices(C3t3& c3t3,
 
       if (vh0_moving)
       {
-        moves[i0] += density * Vector_3(p0, p1);
-        neighbors[i0]++;
-        masses[i0] += density;
+        moves[i0].move += density * Vector_3(p0, p1);
+        moves[i0].mass += density;
+        ++moves[i0].neighbors;
       }
       if (vh1_moving)
       {
-        moves[i1] += density * Vector_3(p1, p0);
-        neighbors[i1]++;
-        masses[i1] += density;
+        moves[i1].move += density * Vector_3(p1, p0);
+        moves[i1].mass += density;
+        ++moves[i1].neighbors;
       }
     }
   }
 
-  // iterate over map of <vertex, id>
-  for (auto [v, vid] : m_vertex_id)
+  // iterate over vertices and move
+  for(Vertex_handle v : tr.finite_vertex_handles())
   {
+    const std::size_t vid = vertex_id(v);
     if (!is_free(vid))
       continue;
 
-    if (c3t3.in_dimension(v) == 3 && neighbors[vid] > 1)
+    if (c3t3.in_dimension(v) == 3 && moves[vid].neighbors > 1)
     {
 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
       os_vol << "2 " << point(v->point());
 #endif
-      const Vector_3 move = moves[vid] / masses[vid];// static_cast<FT>(neighbors[vid]);
+      const Vector_3 move = moves[vid].move / moves[vid].mass;// static_cast<FT>(neighbors[vid]);
       Point_3 new_pos = point(v->point()) + move;
       if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
         nb_done_3d++;

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h

@@ -1439,6 +1439,7 @@ auto sizing_at_vertex(const Vertex_handle v,
 
 template<typename Sizing, typename C3t3, typename Cell_selector>
 auto sizing_at_midpoint(const typename C3t3::Edge& e,
+                        const typename C3t3::Triangulation::Geom_traits::Point_3& m, // edge midpoint
                         const int dim,
                         const typename C3t3::Index& index,
                         const Sizing& sizing,
@@ -1446,17 +1447,13 @@ auto sizing_at_midpoint(const typename C3t3::Edge& e,
                         const Cell_selector& cell_selector)
 {
   using FT = typename C3t3::Triangulation::Geom_traits::FT;
-  using Point_3 = typename C3t3::Triangulation::Geom_traits::Point_3;
-
   auto cp = c3t3.triangulation().geom_traits().construct_point_3_object();
-  const Point_3 m = CGAL::midpoint(cp(e.first->vertex(e.second)->point()),
-                                   cp(e.first->vertex(e.third)->point()));
+
   const FT size = sizing(m, dim, index);
 
   if (dim < 3 && size == 0)
   {
-    const auto u = e.first->vertex(e.second);
-    const auto v = e.first->vertex(e.third);
+    const auto [u, v] = make_vertex_pair(e);
 
     const FT size_at_u = sizing(cp(u->point()), u->in_dimension(), u->index());
     const FT size_at_v = sizing(cp(v->point()), v->in_dimension(), v->index());
@@ -1572,7 +1569,6 @@ auto midpoint_with_info(const typename C3t3::Edge& e,
                         const C3t3& c3t3)
 {
   using Tr = typename C3t3::Triangulation;
-  using Vertex_handle = typename Tr::Vertex_handle;
   using Gt = typename Tr::Geom_traits;
   using Point_3 = typename Gt::Point_3;
   using Index = typename C3t3::Index;
@@ -1584,9 +1580,7 @@ auto midpoint_with_info(const typename C3t3::Edge& e,
     Index index;
   };
 
-  const auto vs = c3t3.triangulation().vertices(e);
-  const Vertex_handle u = vs[0];
-  const Vertex_handle v = vs[1];
+  const auto [u, v] = make_vertex_pair(e);
 
   const auto& gt = c3t3.triangulation().geom_traits();
   auto cp = gt.construct_point_3_object();

Tetrahedral_remeshing/test/Tetrahedral_remeshing/CMakeLists.txt

@@ -16,6 +16,7 @@ create_single_source_cgal_program("test_tetrahedral_remeshing_with_features.cpp"
 create_single_source_cgal_program("test_tetrahedral_remeshing_of_one_subdomain.cpp")
 create_single_source_cgal_program("test_tetrahedral_remeshing_io.cpp")
 create_single_source_cgal_program("test_tetrahedral_remeshing_from_mesh_file.cpp")
+create_single_source_cgal_program("test_tetrahedral_remeshing_determinism.cpp")
 
 # Test MLS projection
 add_executable(test_tetrahedral_remeshing_mls

Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_determinism.cpp

@@ -0,0 +1,104 @@
+//#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+//#define CGAL_TETRAHEDRAL_REMESHING_DEBUG
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation_3.h>
+#include <CGAL/tetrahedral_remeshing.h>
+#include <CGAL/Tetrahedral_remeshing/tetrahedral_remeshing_io.h>
+#include <CGAL/Random.h>
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <cassert>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K> Remeshing_triangulation;
+
+template<typename T3>
+void generate_input_one_subdomain(const std::size_t nbv, T3& tr)
+{
+  CGAL::Random rng;
+
+  typedef typename T3::Point Point;
+  std::vector<Point> pts;
+  while (pts.size() < nbv)
+  {
+    const float x = rng.uniform_real(-10.f, 10.f);
+    const float y = rng.uniform_real(-10.f, 10.f);
+    const float z = rng.uniform_real(-10.f, 10.f);
+
+    pts.push_back(Point(x, y, z));
+  }
+  tr.insert(pts.begin(), pts.end());
+
+  for (typename T3::Cell_handle c : tr.finite_cell_handles())
+    c->set_subdomain_index(1);
+
+  assert(tr.is_valid(true));
+
+#ifdef CGAL_TETRAHEDRAL_REMESHING_GENERATE_INPUT_FILES
+  std::ofstream out("data/triangulation_one_subdomain.binary.cgal",
+                    std::ios_base::out | std::ios_base::binary);
+  CGAL::save_binary_triangulation(out, tr);
+  out.close();
+#endif
+}
+
+int main(int argc, char* argv[])
+{
+  std::cout << "CGAL Random seed = " << CGAL::get_default_random().get_seed() << std::endl;
+
+  // Collect options
+  std::size_t nb_runs = 5;
+
+  Remeshing_triangulation tr;
+  generate_input_one_subdomain(1000, tr);
+
+  const int target_edge_length = (argc > 1) ? atoi(argv[1]) : 1;
+
+  std::ofstream ofs0("in.mesh");
+  CGAL::IO::write_MEDIT(ofs0, tr);
+  ofs0.close();
+
+  const Remeshing_triangulation tr_ref = tr; // Keep a reference for comparison
+
+  std::vector<std::string> output_tr;
+  output_tr.reserve(nb_runs);
+
+  for(std::size_t i = 0; i < nb_runs; ++i)
+  {
+    std::cout << "Run " << i << " of " << nb_runs << std::endl;
+
+    tr = tr_ref; // Reset triangulation to reference
+    CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length);
+
+    std::ostringstream oss;
+    CGAL::IO::write_MEDIT(oss, tr);
+    output_tr.push_back(oss.str());
+    oss.clear();
+
+    if(i == 0)
+      continue; // skip first run, it is the reference
+    else
+    {
+      if(0 != output_tr[i-1].compare(output_tr[i]))
+      {
+        std::cerr << "******************************************" << std::endl;
+        std::cerr << "*** Run " << i << " differs from run " << i-1 << std::endl;
+        assert(false); // This should not happen
+      }
+      else
+      {
+        std::cout << "******************************************" << std::endl;
+        std::cout << "*** Run " << i << " is identical to run " << i - 1 << std::endl;
+        std::cout << "******************************************" << std::endl;
+      }
+    }
+  }
+
+  return EXIT_SUCCESS;
+}