add option --use-new-cavity-algorithm

That option allows to use, or not, the previous implementation of the cavity construction, to compare.

Triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h

@@ -382,6 +382,14 @@ public:
     debug_flags.set(static_cast<int>(Debug_flags::copy_triangulation_into_hole), b);
   }
 
+  bool use_older_cavity_algorithm() const {
+    return debug_flags[static_cast<int>(Debug_flags::use_older_cavity_algorithm)];
+  }
+
+  void use_older_cavity_algorithm(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::use_older_cavity_algorithm), b);
+  }
+
   Vertex_handle insert(const Point &p, Locate_type lt, Cell_handle c,
                        int li, int lj)
   {
@@ -896,6 +904,7 @@ protected:
     missing_region,
     regions,
     copy_triangulation_into_hole,
+    use_older_cavity_algorithm,
     nb_of_flags
   };
   std::bitset<static_cast<int>(Debug_flags::nb_of_flags)> debug_flags{};

Triangulation_3/include/CGAL/Constrained_Delaunay_triangulation_3.h

@@ -1238,6 +1238,7 @@ private:
       };
     };
 
+    auto new_vertex = make__new_element_functor(visited_vertices);
     auto new_cell = make__new_element_functor(visited_cells);
     auto new_edge = [&](Vertex_handle v0, Vertex_handle v1, bool does_intersect) {
       return visited_edges.emplace(CGAL::make_sorted_pair(v0, v1), does_intersect);
@@ -1436,17 +1437,18 @@ private:
         }
 #endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT
         auto [cached_value_it, not_visited] = new_edge(v0, v1, false);
-        if(!not_visited)
+        if(!not_visited) return value_returned(cached_value_it->second);
-          return value_returned(cached_value_it->second);
         int v0v1_intersects_region =
             (v0->is_marked(Vertex_marker::REGION_INSIDE) || v1->is_marked(Vertex_marker::REGION_INSIDE))
                 ? expected
                 : does_edge_intersect_region(cell, index_v0, index_v1, cdt_2, fh_region);
         if(v0v1_intersects_region != 0) {
-          // if(v0v1_intersects_region != expected) {
+          if(this->use_older_cavity_algorithm()) {
-          //   throw PLC_error{"PLC error: v0v1_intersects_region != expected" ,
+            if(v0v1_intersects_region != expected) {
-          //         __FILE__, __LINE__, face_index, region_index};
+              throw PLC_error{"PLC error: v0v1_intersects_region != expected" ,
-          // }
+                    __FILE__, __LINE__, face_index, region_index};
+            }
+          }
           // report the edge with first vertex above the region
           if(v0v1_intersects_region < 0) {
             std::swap(index_v0, index_v1);
@@ -1460,6 +1462,16 @@ private:
         }
       };
 
+      if(this->use_older_cavity_algorithm()) {
+        CGAL_assertion(0 == region_border_vertices.count(v_above));
+        CGAL_assertion(0 == region_border_vertices.count(v_below));
+        if(new_vertex(v_above)) {
+          vertices_of_upper_cavity.push_back(v_above);
+        }
+        if(new_vertex(v_below)) {
+          vertices_of_lower_cavity.push_back(v_below);
+        }
+      }
       auto facet_circ = this->incident_facets(intersecting_edge);
       const auto facet_circ_end = facet_circ;
       do { // loop facets around [v_above, v_below]
@@ -1477,7 +1489,8 @@ private:
         if(region_border_vertices.count(vc) > 0) continue; // intersecting edges cannot touch the border
 
         if(!test_edge(cell, v_above, index_v_above, vc, index_vc, 1) &&
-           !test_edge(cell, v_below, index_v_below, vc, index_vc, -1) && false)
+           !test_edge(cell, v_below, index_v_below, vc, index_vc, -1) &&
+           !this->use_older_cavity_algorithm())
         {
           dump_triangulation();
           dump_region(face_index, region_index, cdt_2);
@@ -1491,7 +1504,7 @@ private:
                 __FILE__, __LINE__, face_index, region_index};
         }
       } while(++facet_circ != facet_circ_end);
-      if(i + 1 == intersecting_edges.size()) {
+      if(!this->use_older_cavity_algorithm() && i + 1 == intersecting_edges.size()) {
         for(auto ch: intersecting_cells) {
           std::cerr << "tetrahedron #" << ch->time_stamp() << " intersects the region\n";
           for(int i = 0; i < 4; ++i) {
@@ -1528,87 +1541,142 @@ private:
             }
           }
         }
+      } // last intersecting edge, and new algorithm
+    } // end loop on intersecting_edges
+    if(this->use_older_cavity_algorithm()) {
+      for(auto intersecting_edge: intersecting_edges) {
+        const auto [v_above, v_below] = tr.vertices(intersecting_edge);
+
+        auto cell_circ = this->incident_cells(intersecting_edge), end = cell_circ;
+        CGAL_assume(cell_circ != nullptr);
+        do {
+          const Cell_handle cell = cell_circ;
+          const auto index_v_above = cell->index(v_above);
+          const auto index_v_below = cell->index(v_below);
+          const auto cell_above = cell->neighbor(index_v_below);
+          const auto cell_below = cell->neighbor(index_v_above);
+          if(0 == intersecting_cells.count(cell_above)) {
+            facets_of_upper_cavity.emplace_back(cell_above, cell_above->index(cell));
+          }
+          if(0 == intersecting_cells.count(cell_below)) {
+            facets_of_lower_cavity.emplace_back(cell_below, cell_below->index(cell));
+          }
+        } while(++cell_circ != end);
       }
-    }
+    } // older algorithm
+
     std::set<Facet> facets_of_border;
-    for(auto c: intersecting_cells) {
-      for(int i = 0; i < 4; ++i) {
-        auto n = c->neighbor(i);
-        if(intersecting_cells.count(n) == 0) {
-          facets_of_border.emplace(c, i);
-        }
-      }
-    }
     Union_find<Vertex_handle> vertices_of_border_union_find;
-    Unique_hash_map<Vertex_handle, typename Union_find<Vertex_handle>::handle> vertices_of_border_handles;
+    if(!this->use_older_cavity_algorithm()) {
-    for(auto facet: facets_of_border) {
+      for(auto c: intersecting_cells) {
-      for(auto v : tr.vertices(facet)) {
+        for(int i = 0; i < 4; ++i) {
-        if(!v->is_marked()) {
+          auto n = c->neighbor(i);
-          v->set_mark(Vertex_marker::CAVITY);
+          if(intersecting_cells.count(n) == 0) {
-          vertices_of_border_handles[v] = vertices_of_border_union_find.make_set(v);
+            facets_of_border.emplace(n, n->index(c));
+          }
         }
       }
-    }
+      Unique_hash_map<Vertex_handle, typename Union_find<Vertex_handle>::handle> vertices_of_border_handles;
-    for(auto facet: facets_of_border) {
+      for(auto facet: facets_of_border) {
-      auto vertices = tr.vertices(facet);
+        for(auto v : tr.vertices(facet)) {
-      for(int i = 0; i < 3; ++i) {
+          if(!v->is_marked()) {
-        auto v1 = vertices[i];
+            v->set_mark(Vertex_marker::CAVITY);
-        auto v2 = vertices[(i + 1) % 3];
+            vertices_of_border_handles[v] = vertices_of_border_union_find.make_set(v);
-        if(v1->is_marked(Vertex_marker::CAVITY) && v2->is_marked(Vertex_marker::CAVITY)) {
+          }
-          vertices_of_border_union_find.unify_sets(vertices_of_border_handles[v1],
-                                                   vertices_of_border_handles[v2]);
         }
       }
+      for(auto facet: facets_of_border) {
+        auto vertices = tr.vertices(facet);
+        for(int i = 0; i < 3; ++i) {
+          auto v1 = vertices[i];
+          auto v2 = vertices[(i + 1) % 3];
+          if(v1->is_marked(Vertex_marker::CAVITY) && v2->is_marked(Vertex_marker::CAVITY)) {
+            vertices_of_border_union_find.unify_sets(vertices_of_border_handles[v1],
+                                                    vertices_of_border_handles[v2]);
+          }
+        }
+      }
+      CGAL_assertion(vertices_of_border_union_find.number_of_sets() <= 2);
+      auto it = vertices_of_border_union_find.begin();
+      const auto vertex_above_handle = it;
+      do {
+        ++it;
+      } while(it != vertices_of_border_union_find.end() &&
+              vertices_of_border_union_find.same_set(vertex_above_handle, it));
+      const auto vertex_below_handle = it;
+      CGAL_assertion((it == vertices_of_border_union_find.end()) ==
+                    (vertices_of_border_union_find.number_of_sets() == 1));
+      for(auto handle = vertices_of_border_union_find.begin(), end = vertices_of_border_union_find.end();
+          handle != end; ++handle)
+      {
+        auto v = *handle;
+        v->clear_mark(Vertex_marker::CAVITY);
+        if(vertices_of_border_union_find.same_set(handle, vertex_above_handle)) {
+          vertices_of_upper_cavity.push_back(v);
+          v->set_mark(Vertex_marker::CAVITY_ABOVE);
+        } else if(vertices_of_border_union_find.same_set(handle, vertex_below_handle)) {
+          vertices_of_lower_cavity.push_back(v);
+          v->set_mark(Vertex_marker::CAVITY_BELOW);
+        } else {
+          CGAL_error();
+        }
+      }
+    } // new algorithm
+    if(this->debug_regions()) {
+      std::stringstream ss_filename;
+      ss_filename << "dump_facets_of_cavity_region_" << face_index << "_" << region_index << "_border.off";
+      std::ofstream out(ss_filename.str());
+      out.precision(17);
+      write_facets(out, tr, facets_of_border);
     }
 
-    std::ofstream out("dump_facets_of_border.off");
+    if(this->debug_regions()) {
-    out.precision(17);
+      for(auto edge : intersecting_edges) {
-    write_facets(out, tr, facets_of_border);
+        auto [v1, v2] = tr.vertices(edge);
-
+        std::cerr << std::format("  edge: {}   {}\n", IO::oformat(v1, with_point_and_info),
-    CGAL_assertion(vertices_of_border_union_find.number_of_sets() <= 2);
+                                IO::oformat(v2, with_point_and_info));
-    auto it = vertices_of_border_union_find.begin();
-    const auto vertex_above_handle = it;
-    do {
-      ++it;
-    } while(it != vertices_of_border_union_find.end() &&
-            vertices_of_border_union_find.same_set(vertex_above_handle, it));
-    const auto vertex_below_handle = it;
-    CGAL_assertion((it == vertices_of_border_union_find.end()) ==
-                   (vertices_of_border_union_find.number_of_sets() == 1));
-    for(auto handle = vertices_of_border_union_find.begin(), end = vertices_of_border_union_find.end();
-        handle != end; ++handle)
-    {
-      auto v = *handle;
-      v->clear_mark(Vertex_marker::CAVITY);
-      if(vertices_of_border_union_find.same_set(handle, vertex_above_handle)) {
-        vertices_of_upper_cavity.push_back(v);
-        v->set_mark(Vertex_marker::CAVITY_ABOVE);
-      } else if(vertices_of_border_union_find.same_set(handle, vertex_below_handle)) {
-        vertices_of_lower_cavity.push_back(v);
-        v->set_mark(Vertex_marker::CAVITY_BELOW);
-      } else {
-        CGAL_error();
       }
     }
-    for(auto facet: facets_of_border) {
+    if(!this->use_older_cavity_algorithm()) {
-      for(auto v: tr.vertices(facet)) {
+      for(auto facet: facets_of_border) {
-        if(v->is_marked(Vertex_marker::CAVITY_ABOVE)) {
+        if(this->debug_regions()) {
-          facets_of_upper_cavity.push_back(facet);
+          std::cerr << "  facet:  ";
-          break;
+          for(auto v: tr.vertices(facet)) {
+            std::cerr << IO::oformat(v, with_point_and_info) << "  ";
+          }
+          std::cerr << "\n";
         }
-        if(v->is_marked(Vertex_marker::CAVITY_BELOW)) {
+        for(auto v: tr.vertices(facet)) {
-          facets_of_lower_cavity.push_back(facet);
+          if(v->is_marked(Vertex_marker::CAVITY_ABOVE)) {
-          break;
+            facets_of_upper_cavity.push_back(facet);
+            break;
+          }
+          if(v->is_marked(Vertex_marker::CAVITY_BELOW)) {
+            facets_of_lower_cavity.push_back(facet);
+            break;
+          }
         }
       }
+      for(auto v: vertices_of_border_union_find)
+      {
+        v->clear_mark(Vertex_marker::CAVITY_ABOVE);
+        v->clear_mark(Vertex_marker::CAVITY_BELOW);
+      }
     }
-    for(auto v: vertices_of_border_union_find)
+    if(this->debug_regions()) {
-    {
+      std::stringstream ss_filename;
-      v->clear_mark(Vertex_marker::CAVITY_ABOVE);
+      ss_filename << "dump_facets_of_upper_cavity_region_" << face_index << "_" << region_index << "_border.off";
-      v->clear_mark(Vertex_marker::CAVITY_BELOW);
+      std::ofstream out(ss_filename.str());
+      out.precision(17);
+      write_facets(out, tr, facets_of_upper_cavity);
+      out.close();
+      ss_filename.str("");
+      ss_filename << "dump_facets_of_lower_cavity_region_" << face_index << "_" << region_index << "_border.off";
+      out.open(ss_filename.str());
+      write_facets(out, tr, facets_of_lower_cavity);
+      out.close();
     }
 
-
     return outputs;
   }
 
@@ -1835,8 +1903,10 @@ private:
     });
 
     auto is_facet_of_polygon = [&](const auto& tr, Facet f) {
-      for(const auto vh: tr.vertices(f)) {
+      const auto [c, facet_index] = f;
-        if(!vh->is_marked()) {
+      for(int i = 0; i < 3; ++i) {
+        const auto vh = c->vertex(T_3::vertex_triple_index(facet_index, i));
+        if(0 == polygon_points.count(tr.point(vh))) {
           return false;
         }
       }

Triangulation_3/test/Triangulation_3/cdt_3_from_off.cpp

@@ -68,6 +68,7 @@ struct CDT_options
   bool debug_missing_regions = false;
   bool debug_regions = false;
   bool debug_copy_triangulation_into_hole = false;
+  bool use_new_cavity_algorithm = false;
   double ratio = 0.;
   double vertex_vertex_epsilon = 1e-6;
   double segment_vertex_epsilon = 1e-8;
@@ -119,6 +120,10 @@ int main(int argc, char* argv[])
       options.merge_facets = true;
     } else if(arg == "--no-merge-facets") {
       options.merge_facets = false;
+    } else if(arg == "--use-new-cavity-algorithm") {
+      options.use_new_cavity_algorithm = true;
+    } else if(arg == "--use-old-cavity-algorithm") {
+      options.use_new_cavity_algorithm = false;
     } else if(arg == "--no-repair") {
       options.repair_mesh = false;
     } else if(arg == "--merge-facets-old") {
@@ -316,18 +321,11 @@ auto segment_soup_to_polylines(Range_of_segments&& segment_soup) {
 
 int go(Mesh mesh, CDT_options options) {
   CDT cdt;
-  if(options.verbose > 0) {
+  cdt.debug_Steiner_points(options.verbose > 0);
-    cdt.debug_Steiner_points(true);
+  cdt.debug_missing_region(options.verbose > 1 || options.debug_missing_regions);
-  }
+  cdt.debug_regions(options.debug_regions);
-  if(options.verbose > 1 || options.debug_missing_regions) {
+  cdt.debug_copy_triangulation_into_hole(options.debug_copy_triangulation_into_hole);
-    cdt.debug_missing_region(true);
+  cdt.use_older_cavity_algorithm(!options.use_new_cavity_algorithm);
-  }
-  if(options.debug_regions) {
-    cdt.debug_regions(true);
-  }
-  if(options.debug_copy_triangulation_into_hole) {
-    cdt.debug_copy_triangulation_into_hole(true);
-  }
   cdt.set_segment_vertex_epsilon(options.segment_vertex_epsilon);
 
   const auto bbox = CGAL::Polygon_mesh_processing::bbox(mesh);