Flip chord topology reorder

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h

@@ -25,6 +25,7 @@
 
 #include <unordered_map>
 #include <unordered_set>
+#include <algorithm>
 #include <limits>
 #include <queue>
 
@@ -88,7 +89,7 @@ void update_c3t3_facets(C3t3& c3t3,
 template<typename C3t3, typename IncCellsVectorMap, typename CellSelector>
 Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge,
                                   C3t3& c3t3,
-                                  const std::vector<typename C3t3::Vertex_handle>& vertices_around_edge,
+                                  const boost::container::small_vector<typename C3t3::Vertex_handle, 16>& vertices_around_edge,
                                   const Flip_Criterion& criterion,
                                   IncCellsVectorMap& inc_cells,
                                   CellSelector& cell_selector)
@@ -578,23 +579,36 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
   Facet_circulator curr_fcirc = tr.incident_facets(edge);
   Facet_circulator curr_fdone = curr_fcirc;
 
-  //Only keep the possible flips
-  std::vector<Vertex_handle> opposite_vertices;
-  int nb_cells_around_edge = 0;
+  //Collect the vertex opposite to the edge in each facet around it, in
+  //circulation order. The chord test below used to re-circulate the facets
+  //around the edge once per such vertex to enumerate the other apices; those
+  //are the same vertices collected here, so the tests can be indexed on the
+  //ring instead. is_edge_uv is read-only, so stopping at the first chord
+  //found yields the same verdict as running the ring to its end.
+  boost::container::small_vector<Vertex_handle, 32> ring_apices;
   do
   {
-    Vertex_handle vh;
-    //Get the ids of the opposite vertices
+    //Get the id of the opposite vertex
     for (int i = 0; i < 3; ++i)
     {
       Vertex_handle curr_vertex = curr_fcirc->first->vertex(
                                     indices(curr_fcirc->second, i));
       if (curr_vertex != vh0 && curr_vertex != vh1)
       {
-        vh = curr_vertex;
+        ring_apices.push_back(curr_vertex);
         break;
       }
     }
+  }
+  while (++curr_fcirc != curr_fdone);
+
+  //Only keep the possible flips
+  std::vector<Vertex_handle> opposite_vertices;
+  int nb_cells_around_edge = 0;
+  const int n_apices = static_cast<int>(ring_apices.size());
+  for (int p = 0; p < n_apices; ++p)
+  {
+    const Vertex_handle vh = ring_apices[p];
 
     if(tr.is_infinite(vh))
       continue;
@@ -603,37 +617,23 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
     if (o_inc_vh.empty())
       tr.incident_cells(vh, std::back_inserter(o_inc_vh));
 
-    Facet_circulator facet_circulator = curr_fcirc;
-    Facet_circulator facet_done = curr_fcirc;
-
-    facet_done--;
-    facet_circulator++;
-    facet_circulator++;
+    //a chord is an edge joining vh to an apex that is not one of its two
+    //neighbors on the ring (positions p-1 and p+1)
     bool is_edge = false;
-    do
+    for (int j = p + 2; j <= p + n_apices - 2; ++j)
     {
-      //Get the ids of the opposite vertices
-      for (int i = 0; i < 3; ++i)
+      if (is_edge_uv(vh, ring_apices[j % n_apices], o_inc_vh))
       {
-        Vertex_handle curr_vertex = facet_circulator->first->vertex(
-                                      indices(facet_circulator->second, i));
-        if (curr_vertex != vh0  && curr_vertex != vh1)
-        {
-          if (is_edge_uv(vh, curr_vertex, o_inc_vh))
-          {
-            is_edge = true;
-            break;
-          }
-        }
+        is_edge = true;
+        break;
       }
-    } while (++facet_circulator != facet_done);
+    }
 
     if (!is_edge)
       opposite_vertices.push_back(vh);
 
     nb_cells_around_edge++;
   }
-  while (++curr_fcirc != curr_fdone);
   if (nb_cells_around_edge < 4)
     return;
 
@@ -1115,8 +1115,12 @@ Sliver_removal_result find_best_flip(typename C3t3::Edge& edge,
   Facet_circulator circ = tr.incident_facets(edge);
   Facet_circulator done = circ;
 
-  //Identify the vertices around this edge
-  std::unordered_set<Vertex_handle> vertices_around_edge;
+  //Identify the vertices around this edge. The ring of apices around an edge
+  //holds a handful of distinct vertices, so they are kept inline and scanned
+  //rather than hashed. flip_3_to_2 does not depend on their order : it picks
+  //vh2/vh3 by testing each vertex against ch0/ch1 individually, and is_facet
+  //is symmetric in its three vertices.
+  boost::container::small_vector<Vertex_handle, 16> vertices_around_edge;
   bool boundary_edge = false;
   bool hull_edge = false;
 
@@ -1130,7 +1134,9 @@ Sliver_removal_result find_best_flip(typename C3t3::Edge& edge,
       Vertex_handle vi = circ->first->vertex(indices(circ->second, i));
       if (vi != v0 && vi != v1)
       {
-        vertices_around_edge.insert(vi);
+        if (std::find(vertices_around_edge.begin(), vertices_around_edge.end(), vi)
+            == vertices_around_edge.end())
+          vertices_around_edge.push_back(vi);
 
         if ( circ->first->subdomain_index()
              != circ->first->neighbor(circ->second)->subdomain_index())
@@ -1161,9 +1167,7 @@ Sliver_removal_result find_best_flip(typename C3t3::Edge& edge,
   {
     if (!boundary_edge && !hull_edge)
     {
-      std::vector<Vertex_handle> vertices;
-      vertices.insert(vertices.end(), vertices_around_edge.begin(), vertices_around_edge.end());
-      res = flip_3_to_2(edge, c3t3, vertices, criterion, inc_cells, cell_selector);
+      res = flip_3_to_2(edge, c3t3, vertices_around_edge, criterion, inc_cells, cell_selector);
     }
   }
   else

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

@@ -976,6 +976,63 @@ std::size_t nb_incident_subdomains(const typename C3t3::Vertex_handle v,
   return indices.size();
 }
 
+// `nb_incident_subdomains(v, c3t3) > 1`, without counting the whole star. The
+// traversal is the one `TDS_3::incident_cells_3()` performs - from v's cell,
+// across the facets that contain v, marking cells as it goes - so it sees the
+// same cells in the same order, and stops as soon as a second index appears.
+// The marks are the TDS's own conflict flags and are cleared before
+// returning, so this must not be called from inside another marking
+// traversal; topology_test, which calls it, is not.
+template<typename C3t3>
+bool has_several_incident_subdomains(const typename C3t3::Vertex_handle v,
+                                     const C3t3& c3t3)
+{
+  typedef typename C3t3::Triangulation::Cell_handle Cell_handle;
+  CGAL_USE(c3t3);
+
+  const Cell_handle start = v->cell();
+  const auto si0 = start->subdomain_index();
+
+  boost::container::small_vector<Cell_handle, 128> marked;
+  boost::container::small_vector<Cell_handle, 128> to_visit;
+
+  start->tds_data().mark_in_conflict();
+  marked.push_back(start);
+  to_visit.push_back(start);
+
+  bool several = false;
+  while (!to_visit.empty())
+  {
+    const Cell_handle c = to_visit.back();
+    to_visit.pop_back();
+
+    if (c->subdomain_index() != si0)
+    {
+      several = true;
+      break;
+    }
+
+    for (int i = 0; i < 4; ++i)
+    {
+      if (c->vertex(i) == v)
+        continue;
+
+      const Cell_handle n = c->neighbor(i);
+      if (!n->tds_data().is_clear())
+        continue;
+
+      n->tds_data().mark_in_conflict();
+      marked.push_back(n);
+      to_visit.push_back(n);
+    }
+  }
+
+  for (const Cell_handle c : marked)
+    c->tds_data().clear();
+
+  return several;
+}
+
 template<typename C3t3>
 std::size_t nb_incident_subdomains(const typename C3t3::Edge& e,
                                    const C3t3& c3t3)
@@ -1312,10 +1369,15 @@ bool topology_test(const typename C3t3::Edge& edge,
       for (int i = 1; i < 4; i++)
       {
         Vertex_handle vi = f.first->vertex((f.second + i) % 4);
-        if (vi != v0 && vi != v1 && nb_incident_subdomains(vi, c3t3) > 1)
+        if (vi != v0 && vi != v1)
         {
-          if (is_edge_in_complex(v0, vi, c3t3)
-              && is_edge_in_complex(v1, vi, c3t3))
+          //(v0,vi) and (v1,vi) are edges of f, so testing them for the
+          //complex needs no is_edge() star walk, and feature edges are rare
+          //enough that the subdomain star walk is skipped almost always.
+          //The three tests are pure, so the conjunction is unchanged.
+          if (c3t3.is_in_complex(v0, vi)
+              && c3t3.is_in_complex(v1, vi)
+              && has_several_incident_subdomains(vi, c3t3))
             return false;
         }
       }