Tetrahedral remeshing sequential refactoring

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/Elementary_operation.h

@@ -0,0 +1,98 @@
+// Copyright (c) 2025 GeometryFactory (France) and Telecom Paris (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)     : Iasonas Manolas, Jane Tournois
+
+#ifndef CGAL_TETRAHEDRAL_REMESHING_ELEMENTARY_OPERATIONS_H
+#define CGAL_TETRAHEDRAL_REMESHING_ELEMENTARY_OPERATIONS_H
+
+#include <CGAL/license/Tetrahedral_remeshing.h>
+
+#include <string>
+
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+#include <CGAL/Real_timer.h>
+#include <cstddef>
+#include <iostream>
+#endif
+
+namespace CGAL {
+namespace Tetrahedral_remeshing {
+namespace internal {
+
+template <typename C3t3_, typename ElementType, typename ElementRange>
+class Elementary_operation
+{
+public:
+  using C3t3 = C3t3_;
+  using Triangulation = typename C3t3::Triangulation;
+  using Element_type = ElementType;
+  using Element_range = ElementRange;
+
+  Elementary_operation() = default;
+  virtual ~Elementary_operation() = default;
+
+  virtual Element_range get_elements(const C3t3& c3t3) const = 0;
+  virtual bool execute_operation(const Element_type& e, C3t3& c3t3) = 0;
+  virtual std::string operation_name() const = 0;
+};
+
+template <typename Operation>
+class Elementary_operation_execution_sequential
+{
+public:
+  using C3t3 = typename Operation::C3t3;
+  using Element_range = typename Operation::Element_range;
+
+  bool execute(Operation& op, C3t3& c3t3) const
+  {
+    Element_range candidates = op.get_elements(c3t3);
+    if (candidates.empty())
+      return false;
+
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+    std::size_t nb_done = 0;
+    const std::size_t nb_candidates = candidates.size();
+    CGAL::Real_timer timer;
+    timer.start();
+#endif
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS
+    std::size_t nb_processed = 0;
+#endif
+    for (const auto& element : candidates)
+    {
+      if (op.execute_operation(element, c3t3))
+      {
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+        ++nb_done;
+#endif
+      }
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS
+      std::cout << "\r" << op.operation_name() << "... ("
+                << ++nb_processed << "/" << nb_candidates << ")";
+      std::cout.flush();
+#endif
+    }
+
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+    timer.stop();
+    std::cout << op.operation_name() << ": " << nb_done << "/"
+              << nb_candidates << " done ("
+              << timer.time() << " sec)." << std::endl;
+#endif
+    return true;
+  }
+};
+
+} // namespace internal
+} // namespace Tetrahedral_remeshing
+} // namespace CGAL
+
+#endif // CGAL_TETRAHEDRAL_REMESHING_ELEMENTARY_OPERATIONS_H

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h

@@ -21,6 +21,7 @@
 #include <boost/container/small_vector.hpp>
 #include <boost/container/flat_set.hpp>
 #include <boost/functional/hash.hpp>
+#include <boost/unordered_set.hpp>
 
 #include <vector>
 #include <algorithm>
@@ -31,6 +32,7 @@
 
 #include <CGAL/SMDS_3/tet_soup_to_c3t3.h>
 #include <CGAL/utility.h>
+#include <CGAL/Tetrahedral_remeshing/internal/Elementary_operation.h>
 #include <CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h>
 
 #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
@@ -977,7 +979,7 @@ collapse(const typename C3t3::Cell_handle ch,
 
 
 template<typename C3t3, typename CellSelector, typename ShortEdgesBimap>
-typename C3t3::Vertex_handle collapse(typename C3t3::Edge& edge,
+typename C3t3::Vertex_handle collapse(const typename C3t3::Edge& edge,
                                       const Collapse_type& collapse_type,
                                       CellSelector& cell_selector,
                                       C3t3& c3t3,
@@ -1084,7 +1086,7 @@ template<typename C3t3,
          typename CellSelector,
          typename ShortEdgesBimap,
          typename Visitor>
-typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge,
+typename C3t3::Vertex_handle collapse_edge(const typename C3t3::Edge& edge,
     C3t3& c3t3,
     const Sizing& sizing,
     const bool /* protect_boundaries */,
@@ -1237,125 +1239,163 @@ auto can_be_collapsed(const typename C3T3::Edge& e,
   return Collapsible {true, boundary};
 }
 
-template<typename C3T3,
-         typename Sizing,
+// The short edges left to collapse, shortest first. Edges are compared by
+// their vertex pair, but stored with their orientation : `collapse_edge()`
+// reads it to decide which extremity survives, so an edge already in the map
+// keeps the orientation it entered with, and only its length is updated.
+template<typename C3t3>
+using Short_edges_bimap = boost::bimap<
+    boost::bimaps::set_of<typename C3t3::Triangulation::Edge,
+                          Compare_edges<typename C3t3::Triangulation::Edge> >,
+    boost::bimaps::multiset_of<typename C3t3::Triangulation::Geom_traits::FT,
+                               std::less<typename C3t3::Triangulation::Geom_traits::FT> > >;
+
+template<typename C3t3,
+         typename SizingFunction,
          typename CellSelector,
          typename Visitor>
-void collapse_short_edges(C3T3& c3t3,
-                          const Sizing& sizing,
-                          const bool protect_boundaries,
-                          CellSelector cell_selector,
-                          Visitor& visitor)
+class Edge_collapse_operation
+    : public Elementary_operation<C3t3,
+                                 typename C3t3::Triangulation::Edge,
+                                 Short_edges_bimap<C3t3> >
 {
-  typedef typename C3T3::Triangulation       T3;
-  typedef typename T3::Edge                  Edge;
-  typedef typename T3::Vertex_handle         Vertex_handle;
-
-  typedef typename T3::Geom_traits::FT FT;
-  typedef boost::bimap<
-        boost::bimaps::set_of<Edge, Compare_edges<Edge> >,
-        boost::bimaps::multiset_of<FT, std::less<FT> > >  Boost_bimap;
-  typedef typename Boost_bimap::value_type            short_edge;
-
-  T3& tr = c3t3.triangulation();
-
-#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
-  std::cout << "Collapse short edges...";
-  std::cout.flush();
-  std::size_t nb_collapses = 0;
-  CGAL::Real_timer timer;
-  timer.start();
-#endif
+public:
+  using Tr = typename C3t3::Triangulation;
+  using Vertex_handle = typename Tr::Vertex_handle;
+  using Edge = typename Tr::Edge;
+  using FT = typename Tr::Geom_traits::FT;
+
+  using Short_edges = Short_edges_bimap<C3t3>;
+  using Base_operation = Elementary_operation<C3t3, Edge, Short_edges>;
+  using Element_type = typename Base_operation::Element_type;
+  using Element_range = typename Base_operation::Element_range;
+
+private:
+  const SizingFunction& m_sizing;
+  const CellSelector& m_cell_selector;
+  bool m_protect_boundaries;
+  Visitor& m_visitor;
 
-  //collect long edges
-  Boost_bimap short_edges;
-  for (const Edge& e : tr.finite_edges())
+public:
+  Edge_collapse_operation(const SizingFunction& sizing,
+                        const CellSelector& cell_selector,
+                        const bool protect_boundaries,
+                        Visitor& visitor)
+      : m_sizing(sizing)
+      , m_cell_selector(cell_selector)
+      , m_protect_boundaries(protect_boundaries)
+      , m_visitor(visitor) {}
+
+  Element_range get_elements(const C3t3& c3t3) const override
   {
-    auto [collapsible, boundary] = can_be_collapsed(e, c3t3, protect_boundaries, cell_selector);
-    if (!collapsible)
-      continue;
-
-    const auto sqlen = is_too_short(e, boundary, sizing, c3t3, cell_selector);
-    if(sqlen != std::nullopt)
-      short_edges.insert(short_edge(e, sqlen.value()));
+    Short_edges short_edges;
+    for (const Edge& e : c3t3.triangulation().finite_edges())
+    {
+      auto [collapsible, boundary]
+        = can_be_collapsed(e, c3t3, m_protect_boundaries, m_cell_selector);
+      if (!collapsible)
+        continue;
+
+      const auto sqlen = is_too_short(e, boundary, m_sizing, c3t3, m_cell_selector);
+      if (sqlen != std::nullopt)
+        short_edges.insert(typename Short_edges::value_type(e, sqlen.value()));
+    }
+    return short_edges;
   }
 
-#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
-  debug::dump_edges(short_edges, "short_edges.polylines.txt");
-
-  std::ofstream short_success("short_collapse_success.polylines.txt");
-  std::ofstream short_fail("short_collapse_fail.polylines.txt");
-  std::ofstream short_cancel("short_collapse_canceled.polylines.txt");
-#endif
+  bool execute_operation(const Element_type& edge, C3t3& c3t3) override
+  {
+    Short_edges no_short_edges; // no work list to keep up to date
+    return execute_operation(edge, c3t3, no_short_edges);
+  }
 
-  while(!short_edges.empty())
+  /**
+  * Collapses `edge`, and keeps `short_edges` up to date : `collapse_edge()`
+  * removes from it the edges it destroys, and the edges incident to the
+  * vertex it keeps are re-evaluated here, since their length has changed.
+  */
+  bool execute_operation(const Element_type& edge, C3t3& c3t3,
+                         Short_edges& short_edges)
   {
-    //the edge with shortest length
-    typename Boost_bimap::right_map::iterator eit = short_edges.right.begin();
-    Edge e = eit->second;
-
-#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS
-    FT sqlen = eit->first;
-    std::cout << "\rCollapse... (" << short_edges.left.size() << " short edges, ";
-    std::cout << std::sqrt(sqlen) << ", ";
-    std::cout << nb_collapses << " collapses)";
-    std::cout.flush();
-#endif
+    const Vertex_handle vh = collapse_edge(edge, c3t3, m_sizing, m_protect_boundaries,
+                                           m_cell_selector, short_edges, m_visitor);
+    if (vh == Vertex_handle())
+      return false;
 
-    short_edges.right.erase(eit);
+    std::vector<Edge> incident_short;
+    c3t3.triangulation().finite_incident_edges(vh, std::back_inserter(incident_short));
+    for (const Edge& eshort : incident_short)
+    {
+      const auto [collapsible, boundary]
+        = can_be_collapsed(eshort, c3t3, m_protect_boundaries, m_cell_selector);
 
-    CGAL_expensive_assertion_code(const bool bd = is_boundary_edge(e));
-    CGAL_expensive_assertion(!!is_too_short(e, bd, sizing, c3t3, cell_selector));
-    CGAL_expensive_assertion(can_be_collapsed(e, c3t3, protect_boundaries, cell_selector));
+      // an edge that can no longer be collapsed leaves the work list, rather
+      // than being taken out of it later and refused by collapse_edge()
+      std::optional<FT> sqlen;
+      if (collapsible)
+        sqlen = is_too_short(eshort, boundary, m_sizing, c3t3, m_cell_selector);
 
-#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
-    const auto p1 = e.first->vertex(e.second)->point();
-    const auto p2 = e.first->vertex(e.third)->point();
-#endif
+      update_bimap(eshort, short_edges, sqlen);
+    }
+    return true;
+  }
 
-    Vertex_handle vh = collapse_edge(e, c3t3, sizing,
-                                     protect_boundaries, cell_selector,
-                                     short_edges,
-                                     visitor);
-    if (vh != Vertex_handle())
-    {
-      std::vector<Edge> incident_short;
-      c3t3.triangulation().finite_incident_edges(vh,
-          std::back_inserter(incident_short));
-      for (const Edge& eshort : incident_short)
-      {
-        const auto [collapsible, boundary]
-          = can_be_collapsed(eshort, c3t3, protect_boundaries, cell_selector);
-        if (!collapsible)
-          continue;
+  std::string operation_name() const override { return "Collapse short edges"; }
+};
 
-        const auto sqlen = is_too_short(eshort, boundary, sizing, c3t3, cell_selector);
-        update_bimap(eshort, short_edges, sqlen);
-      }
+/**
+* Collapse is the only operation whose elements change as it runs : collapsing
+* an edge shortens the ones around the vertex it keeps, and destroys others.
+* Its elements are therefore taken from a work list that `execute_operation()`
+* keeps up to date, shortest first, rather than from a list collected once.
+*/
+template<typename C3t3,
+         typename SizingFunction,
+         typename CellSelector,
+         typename Visitor>
+class Elementary_operation_execution_sequential<
+        Edge_collapse_operation<C3t3, SizingFunction, CellSelector, Visitor> >
+{
+  using Operation = Edge_collapse_operation<C3t3, SizingFunction, CellSelector, Visitor>;
+  using Short_edges = typename Operation::Short_edges;
+  using Edge = typename Operation::Edge;
+
+public:
+  bool execute(Operation& op, C3t3& c3t3) const
+  {
+    Short_edges short_edges = op.get_elements(c3t3);
+    if (short_edges.empty())
+      return false;
 
 #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
-      ++nb_collapses;
+    std::size_t nb_done = 0;
+    CGAL::Real_timer timer;
+    timer.start();
 #endif
+    while (!short_edges.empty())
+    {
+      // the edge with shortest length
+      typename Short_edges::right_map::iterator eit = short_edges.right.begin();
+      const Edge e = eit->second;
+      short_edges.right.erase(eit);
 
-#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
-      if (vh != Vertex_handle())
-        short_success << "2 " << point(p1) << " " << point(p2) << std::endl;
-      else
-        short_fail << "2 " << point(p1) << " " << point(p2) << std::endl;
+      if (op.execute_operation(e, c3t3, short_edges))
+      {
+#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
+        ++nb_done;
 #endif
+      }
     }
-  }//end loop on short_edges
-#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
-  short_success.close();
-  short_fail.close();
-#endif
 
 #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
-  timer.stop();
-  std::cout << " done (" << nb_collapses << " collapses, in "
-            << timer.time() << " seconds)." << std::endl;
+    timer.stop();
+    std::cout << op.operation_name() << ": " << nb_done << " done ("
+              << timer.time() << " sec)." << std::endl;
 #endif
-}
+    return true;
+  }
+};
+
 }
 }
 }