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)
@@ -569,23 +570,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;
@@ -594,37 +608,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;
 
@@ -1097,8 +1097,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;
 
@@ -1112,7 +1116,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())
@@ -1143,9 +1149,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

@@ -1308,10 +1308,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)
+              && nb_incident_subdomains(vi, c3t3) > 1)
             return false;
         }
       }