From 02730495418741d588199883d94aa50564de4b8c Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 12 Sep 2025 11:31:45 +0200 Subject: [PATCH 01/21] use if-constexpr instead of C++ preprocessor --- ...ing_constrained_Delaunay_triangulation_3.h | 49 ++++++------------- 1 file changed, 14 insertions(+), 35 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 60029fab9ce..7bf38f39750 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -71,6 +71,7 @@ #include #include #include +#include #include #include #include @@ -2105,18 +2106,15 @@ private: for(std::size_t i = 0; i < intersecting_edges.size(); ++i) { const auto intersecting_edge = intersecting_edges[i]; const auto [v_above, v_below] = tr().vertices(intersecting_edge); -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_regions()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { std::cerr << cdt_3_format("restore_subface_region face index: {}, region #{}, intersecting edge #{}: ({} {})\n", face_index, region_index, i, IO::oformat(v_above, with_point_and_info), IO::oformat(v_below, with_point_and_info)); dump_region(face_index, region_index, cdt_2); } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_regions()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { const auto p_above = this->point(v_above); const auto p_below = this->point(v_below); const auto edge_segment = typename Geom_traits::Segment_3{p_above, p_below}; @@ -2180,25 +2178,20 @@ private: } } } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT auto test_edge = [&](Cell_handle cell, Vertex_handle v0, int index_v0, Vertex_handle v1, int index_v1, [[maybe_unused]] int expected) { auto value_returned = [this](bool b) { CGAL_USE(this); -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_regions()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { std::cerr << cdt_3_format(" return {}\n", b); } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT return b; }; -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_regions()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { std::cerr << cdt_3_format(" test_edge {} {} ", IO::oformat(v0, with_point_and_info), IO::oformat(v1, with_point_and_info)); } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT auto [cached_value_it, not_visited] = new_edge(v0, v1, false); if(!not_visited) return value_returned(cached_value_it->second); int v0v1_intersects_region = (v0->ccdt_3_data().is_marked(Vertex_marker::REGION_INSIDE) || @@ -2593,14 +2586,12 @@ private: } return vertices; }); -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_regions()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { std::cerr << "region_border_vertices.size() = " << region_border_vertices.size() << "\n"; for(auto v : region_border_vertices) { std::cerr << cdt_3_format(" {}\n", IO::oformat(v, with_point)); } } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT for(auto v: region_border_vertices) { v->ccdt_3_data().set_mark(Vertex_marker::REGION_BORDER); } @@ -2669,7 +2660,6 @@ private: } } -#if CGAL_CAN_USE_CXX20_FORMAT if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { std::cerr << cdt_3_format ("NOTE: diagonal: {:.6} {:.6} {} in tr\n", @@ -2710,7 +2700,6 @@ private: } } } -#endif // CGAL_CAN_USE_CXX20_FORMAT } return false; }; @@ -2788,8 +2777,7 @@ private: return true; }; -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_regions()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { std::cerr << cdt_3_format("Cavity has {} cells and {} edges, " "{} vertices in upper cavity and {} in lower, " "{} facets in upper cavity and {} in lower\n", @@ -2814,7 +2802,6 @@ private: // dump_facets_of_cavity(face_index, region_index, "upper", original_facets_of_upper_cavity); } } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT auto register_internal_constrained_facet = [this](Facet f) { this->register_facet_to_be_constrained(f); }; if(this->debug_copy_triangulation_into_hole()) { @@ -3003,8 +2990,7 @@ private: const auto upper_inner_map = tr().create_triangulation_inner_map( upper_cavity_triangulation, map_upper_cavity_vertices_to_ambient_vertices, false); -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_copy_triangulation_into_hole()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_copy_triangulation_into_hole()) { std::cerr << "upper_inner_map:\n"; for(auto [vt, _] : upper_inner_map) { std::cerr << cdt_3_format(" {:.6}, {:.6}, {:.6})\n", @@ -3013,7 +2999,6 @@ private: IO::oformat(vt[2], with_point)); } } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT if(this->debug_copy_triangulation_into_hole()) { std::cerr << "# glu the lower triangulation of the cavity\n"; } @@ -3042,8 +3027,8 @@ private: { const auto lower_inner_map = tr().create_triangulation_inner_map( lower_cavity_triangulation, map_lower_cavity_vertices_to_ambient_vertices, false); -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_copy_triangulation_into_hole()) { +#if CGAL_CAN_USE_CXX20_FORMAT + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_copy_triangulation_into_hole()) { std::cerr << "outer_map:\n"; for(auto [vt, _] : outer_map) { std::cerr << cdt_3_format(" {:.6}, {:.6}, {:.6})\n", @@ -3056,7 +3041,7 @@ private: write_facets(out, *this, std::ranges::views::values(outer_map)); out.close(); } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT +#endif // CGAL_CAN_USE_CXX20_FORMAT this->copy_triangulation_into_hole(map_lower_cavity_vertices_to_ambient_vertices, std::move(outer_map), lower_inner_map, this->new_cells_output_iterator()); } @@ -3217,14 +3202,12 @@ private: tr().is_infinite(v) ? cavity_triangulation.infinite_vertex() : cavity_triangulation.insert(this->point(v)); map_ambient_vertices_to_cavity_vertices[v] = cavity_v; map_cavity_vertices_to_ambient_vertices[cavity_v] = v; -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_regions()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { std::cerr << cdt_3_format("inserted {}cavity vertex {:.6} -> {:.6}\n", extra, IO::oformat(cavity_v, with_point_and_info), IO::oformat(v, with_point_and_info)); } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT return cavity_v; }; @@ -3464,25 +3447,21 @@ private: bool restore_face(CDT_3_signed_index face_index) { CDT_2& non_const_cdt_2 = face_cdt_2[face_index]; const CDT_2& cdt_2 = non_const_cdt_2; -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_copy_triangulation_into_hole()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_copy_triangulation_into_hole()) { std::cerr << cdt_3_format("restore_face({}): CDT_2 has {} vertices\n", face_index, cdt_2.number_of_vertices()); } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT for(const auto& edge : cdt_2.finite_edges()) { const auto fh = edge.first; const auto i = edge.second; const auto va_3d = fh->vertex(cdt_2.cw(i))->info().vertex_handle_3d; const auto vb_3d = fh->vertex(cdt_2.ccw(i))->info().vertex_handle_3d; const bool is_3d = this->is_edge(va_3d, vb_3d); -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - if(this->debug_copy_triangulation_into_hole()) { + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_copy_triangulation_into_hole()) { std::cerr << cdt_3_format("Edge is 3D: {:6} ({} , {})\n", is_3d, IO::oformat(va_3d, with_point_and_info), IO::oformat(vb_3d, with_point_and_info)); } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT CGAL_assertion(is_3d || !cdt_2.is_constrained(edge)); fh->info().is_edge_also_in_3d_triangulation[unsigned(i)] = is_3d; const auto reverse_edge = cdt_2.mirror_edge(edge); From 0b2ebbc23e7c510cd6adc2ca848f8ea5c3b523d2 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 18 Sep 2025 16:57:54 +0200 Subject: [PATCH 02/21] extract debug output functions from construct_cavities --- .../Conforming_Delaunay_triangulation_3.h | 4 + ...ing_constrained_Delaunay_triangulation_3.h | 423 ++++++++++-------- 2 files changed, 241 insertions(+), 186 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h index a8b4b4edea4..b98eccbb60e 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h @@ -369,6 +369,10 @@ public: return debug_flags[static_cast(Debug_flags::use_older_cavity_algorithm)]; } + bool use_newer_cavity_algorithm() const { + return !debug_flags[static_cast(Debug_flags::use_older_cavity_algorithm)]; + } + void use_older_cavity_algorithm(bool b) { debug_flags.set(static_cast(Debug_flags::use_older_cavity_algorithm), b); } diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 7bf38f39750..8e40019e09c 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2015,90 +2015,9 @@ private: return visited_edges.emplace(CGAL::make_sorted_pair(v0, v1), does_intersect); }; -#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT - using Mesh = Surface_mesh; - using Face_index = typename Mesh::Face_index; - using EK = CGAL::Exact_predicates_exact_constructions_kernel; - const auto to_exact = CGAL::Cartesian_converter(); - const auto from_exact = CGAL::Cartesian_converter(); - if(this->debug_regions()) { - - Mesh tets_intersect_region_mesh; - auto [color_vpmap, _] = tets_intersect_region_mesh.template add_property_map("f:patch_id"); - - for(auto ch : tr().finite_cell_handles()) { - auto tetrahedron = typename Geom_traits::Tetrahedron_3{tr().point(ch->vertex(0)), tr().point(ch->vertex(1)), - tr().point(ch->vertex(2)), tr().point(ch->vertex(3))}; - if(!std::any_of(fh_region.begin(), fh_region.end(), [&](auto fh) { - const auto v0 = fh->vertex(0)->info().vertex_handle_3d; - const auto v1 = fh->vertex(1)->info().vertex_handle_3d; - const auto v2 = fh->vertex(2)->info().vertex_handle_3d; - const auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; - return does_tetrahedron_intersect_triangle_interior(tetrahedron, triangle, tr().geom_traits()); - })) - { - continue; - } - bool intersects = false; - for(int i = 0; i < 4; ++i) { - for(int j = i + 1; j < 4; ++j) { - int intersects_region = does_edge_intersect_region(ch, i, j, cdt_2, fh_region); - if(intersects_region != 0) { - intersects = true; - } - } - } - if(!intersects) { - std::cerr << "ERROR: tetrahedron #" << ch->time_stamp() << " has no edge intersecting the region\n"; - } - std::ofstream dump_tetrahedron( - cdt_3_format("dump_intersecting_{}_{}_tetrahedron_{}.off", face_index, region_index, ch->time_stamp())); - dump_tetrahedron.precision(17); - Mesh mesh; - CGAL::make_tetrahedron(tr().point(ch->vertex(0)), tr().point(ch->vertex(1)), tr().point(ch->vertex(2)), - tr().point(ch->vertex(3)), mesh); - dump_tetrahedron << mesh; - dump_tetrahedron.close(); - - auto exact_tetrahedron = to_exact(tetrahedron); - for(auto fh : fh_region) { - auto v0 = fh->vertex(0)->info().vertex_handle_3d; - auto v1 = fh->vertex(1)->info().vertex_handle_3d; - auto v2 = fh->vertex(2)->info().vertex_handle_3d; - auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; - - auto exact_triangle = to_exact(triangle); - auto tetrahedron_triangle_intersection_opt = CGAL::intersection(exact_tetrahedron, exact_triangle); - if(!tetrahedron_triangle_intersection_opt) { - continue; - } - if(const auto* tri = std::get_if(&tetrahedron_triangle_intersection_opt.value())) { - exact(*tri); - auto v0 = tets_intersect_region_mesh.add_vertex(from_exact((*tri)[0])); - auto v1 = tets_intersect_region_mesh.add_vertex(from_exact((*tri)[1])); - auto v2 = tets_intersect_region_mesh.add_vertex(from_exact((*tri)[2])); - std::array arr{v0, v1, v2}; - auto f = CGAL::Euler::add_face(arr, tets_intersect_region_mesh); - put(color_vpmap, f, static_cast(ch->time_stamp())); - } - if(const auto* vec = std::get_if>(&tetrahedron_triangle_intersection_opt.value())) - { - std::vector vec_of_indices; - for(const auto& p : *vec) { - exact(p); - vec_of_indices.push_back(tets_intersect_region_mesh.add_vertex(from_exact(p))); - } - CGAL::Euler::add_face(vec_of_indices, tets_intersect_region_mesh); - } - } - } - std::ofstream tets_intersect_region_out( - cdt_3_format("dump_tets_intersect_region_{}_{}.ply", face_index, region_index)); - tets_intersect_region_out.precision(17); - CGAL::IO::write_PLY(tets_intersect_region_out, tets_intersect_region_mesh); - tets_intersect_region_out.close(); + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { + debug_dump_tetrahedra_intersect_region(face_index, region_index, cdt_2, fh_region); } -#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT intersecting_edges.push_back(first_intersecting_edge); const auto [v0, v1] = tr().vertices(first_intersecting_edge); @@ -2107,76 +2026,7 @@ private: const auto intersecting_edge = intersecting_edges[i]; const auto [v_above, v_below] = tr().vertices(intersecting_edge); if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { - std::cerr << cdt_3_format("restore_subface_region face index: {}, region #{}, intersecting edge #{}: ({} {})\n", - face_index, region_index, i, - IO::oformat(v_above, with_point_and_info), - IO::oformat(v_below, with_point_and_info)); - dump_region(face_index, region_index, cdt_2); - } - - if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { - const auto p_above = this->point(v_above); - const auto p_below = this->point(v_below); - const auto edge_segment = typename Geom_traits::Segment_3{p_above, p_below}; - const auto exact_edge_segment = to_exact(edge_segment); - - std::ofstream intersect_out("dump_edge_region_intersection.xyz"); - intersect_out.precision(17); - for(auto fh: fh_region) { - auto v0 = fh->vertex(0)->info().vertex_handle_3d; - auto v1 = fh->vertex(1)->info().vertex_handle_3d; - auto v2 = fh->vertex(2)->info().vertex_handle_3d; - auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; - auto exact_triangle = to_exact(triangle); - if(auto edge_intersection_opt = CGAL::intersection(exact_edge_segment, exact_triangle)) { - const auto& edge_intersection = *edge_intersection_opt; - if(const auto* p = std::get_if(&edge_intersection)) { - exact(*p); - intersect_out << *p << '\n'; - } - } - } - intersect_out.close(); - - auto cells_around_intersecting_edge = Container_from_circulator{this->incident_cells(intersecting_edge)}; - for(const auto& cell: cells_around_intersecting_edge) { - CGAL_assertion(!cell.has_vertex(tr().infinite_vertex())); - auto tetrahedron = - typename Geom_traits::Tetrahedron_3{tr().point(cell.vertex(0)), tr().point(cell.vertex(1)), - tr().point(cell.vertex(2)), tr().point(cell.vertex(3))}; - for(auto fh: fh_region) { - auto v0 = fh->vertex(0)->info().vertex_handle_3d; - auto v1 = fh->vertex(1)->info().vertex_handle_3d; - auto v2 = fh->vertex(2)->info().vertex_handle_3d; - auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; - auto exact_triangle = to_exact(triangle); - } - - std::cerr << cdt_3_format("Test tetrahedron (#{}):\n {}\n {}\n {}\n {}\n", - cell.time_stamp(), - IO::oformat(cell.vertex(0), with_point_and_info), - IO::oformat(cell.vertex(1), with_point_and_info), - IO::oformat(cell.vertex(2), with_point_and_info), - IO::oformat(cell.vertex(3), with_point_and_info)); - if(!std::any_of(fh_region.begin(), fh_region.end(), [&](const auto fh) { - auto v0 = fh->vertex(0)->info().vertex_handle_3d; - auto v1 = fh->vertex(1)->info().vertex_handle_3d; - auto v2 = fh->vertex(2)->info().vertex_handle_3d; - auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; - bool b = does_tetrahedron_intersect_triangle_interior(tetrahedron, triangle, tr().geom_traits()); - if(b) { - std::cerr << " intersects the region\n"; - } - return b; - })) - { - std::cerr << cdt_3_format( - "ERROR: The following tetrahedron (#{}) does not intersect the region:\n {}\n {}\n {}\n {}\n", - cell.time_stamp(), - IO::oformat(cell.vertex(0), with_point_and_info), IO::oformat(cell.vertex(1), with_point_and_info), - IO::oformat(cell.vertex(2), with_point_and_info), IO::oformat(cell.vertex(3), with_point_and_info)); - } - } + debug_dump_edge_region_intersection(face_index, region_index, fh_region, i, v_above, v_below, intersecting_edge); } auto test_edge = [&](Cell_handle cell, Vertex_handle v0, int index_v0, Vertex_handle v1, int index_v1, @@ -2265,7 +2115,7 @@ private: __FILE__, __LINE__, face_index, region_index}; } } while(++facet_circ != facet_circ_end); - if(!this->use_older_cavity_algorithm() && i + 1 == intersecting_edges.size()) { + if(this->use_newer_cavity_algorithm() && i + 1 == intersecting_edges.size()) { for(auto ch: intersecting_cells) { if(this->debug_regions()) { std::cerr << "tetrahedron #" << ch->time_stamp() << " intersects the region\n"; @@ -2339,7 +2189,7 @@ private: std::set facets_of_border; Union_find vertices_of_cavity_union_find; - if(!this->use_older_cavity_algorithm()) { + if(this->use_newer_cavity_algorithm()) { for(auto c: intersecting_cells) { for(int i = 0; i < 4; ++i) { auto n = c->neighbor(i); @@ -2474,11 +2324,7 @@ private: } } // 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); + debug_dump_cavity_outputs(face_index, region_index, intersecting_edges, facets_of_border, facets_of_upper_cavity, facets_of_lower_cavity); } if(this->debug_regions()) { @@ -2488,18 +2334,10 @@ private: IO::oformat(v2, with_point_and_info)); } } - if(!this->use_older_cavity_algorithm()) { + if(this->use_newer_cavity_algorithm()) { for(auto facet: facets_of_border) { if(this->debug_regions()) { - std::cerr << " facet: "; - const auto facet_vertices = tr().vertices(facet); - for(auto v: facet_vertices) { - std::cerr << IO::oformat(v, with_point_and_info) << " "; - } - // This assertion is wrong, because there might be only one half-cavity and not a full cavity. - // CGAL_assertion(!std::all_of(facet_vertices.begin(), facet_vertices.end(), - // [](auto v) { return v->is_marked(Vertex_marker::REGION_BORDER); })); - std::cerr << "\n"; + debug_output_facet_vertices({facet}); } for(auto v: tr().vertices(facet)) { if(v->ccdt_3_data().is_marked(Vertex_marker::CAVITY_ABOVE)) { @@ -2519,19 +2357,6 @@ private: v->ccdt_3_data().clear_mark(Vertex_marker::CAVITY_BELOW); } } - if(this->debug_regions()) { - std::stringstream ss_filename; - ss_filename << "dump_facets_of_upper_cavity_region_" << face_index << "_" << region_index << "_border.off"; - 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; } @@ -2591,6 +2416,11 @@ private: for(auto v : region_border_vertices) { std::cerr << cdt_3_format(" {}\n", IO::oformat(v, with_point)); } + std::ofstream out(cdt_3_format("dump_border_edges_region_{}_{}.polylines.txt", face_index, region_index)); + out.precision(17); + for(auto edge : border_edges) { + write_segment(out, edge); + } } for(auto v: region_border_vertices) { v->ccdt_3_data().set_mark(Vertex_marker::REGION_BORDER); @@ -3027,8 +2857,8 @@ private: { const auto lower_inner_map = tr().create_triangulation_inner_map( lower_cavity_triangulation, map_lower_cavity_vertices_to_ambient_vertices, false); -#if CGAL_CAN_USE_CXX20_FORMAT - if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_copy_triangulation_into_hole()) { +#if CGAL_CDT_3_CAN_USE_CXX20_FORMAT + if(this->debug_copy_triangulation_into_hole()) { std::cerr << "outer_map:\n"; for(auto [vt, _] : outer_map) { std::cerr << cdt_3_format(" {:.6}, {:.6}, {:.6})\n", @@ -3041,7 +2871,7 @@ private: write_facets(out, *this, std::ranges::views::values(outer_map)); out.close(); } -#endif // CGAL_CAN_USE_CXX20_FORMAT +#endif // CGAL_CDT_3_CAN_USE_CXX20_FORMAT this->copy_triangulation_into_hole(map_lower_cavity_vertices_to_ambient_vertices, std::move(outer_map), lower_inner_map, this->new_cells_output_iterator()); } @@ -3773,6 +3603,227 @@ public: write_facets(out, tr, tr.finite_facets()); } + void debug_dump_cavity_outputs(CDT_3_signed_index face_index, + int region_index, + const std::vector& intersecting_edges, + const std::set& facets_of_border, + const std::vector& facets_of_upper_cavity, + const std::vector& facets_of_lower_cavity) + { + // Dump facets of cavity border + 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); + out.close(); + + // Dump intersecting edges information + for(auto edge : intersecting_edges) { + auto [v1, v2] = tr().vertices(edge); + std::cerr << cdt_3_format(" edge: {} {}\n", IO::oformat(v1, with_point_and_info), + IO::oformat(v2, with_point_and_info)); + } + + // Dump upper and lower cavity facets + ss_filename.str(""); + ss_filename << "dump_facets_of_upper_cavity_region_" << face_index << "_" << region_index << "_border.off"; + out.open(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(); + } + + void debug_output_facet_vertices(const std::set& facets_of_border) + { + for(auto facet: facets_of_border) { + std::cerr << " facet: "; + const auto facet_vertices = tr().vertices(facet); + for(auto v: facet_vertices) { + std::cerr << IO::oformat(v, with_point_and_info) << " "; + } + // This assertion is wrong, because there might be only one half-cavity and not a full cavity. + // CGAL_assertion(!std::all_of(facet_vertices.begin(), facet_vertices.end(), + // [](auto v) { return v->is_marked(Vertex_marker::REGION_BORDER); })); + std::cerr << "\n"; + } + } + + template + void debug_dump_edge_region_intersection(CDT_3_signed_index face_index, + int region_index, + const Fh_region& fh_region, + std::size_t edge_index, + Vertex_handle v_above, + Vertex_handle v_below, + Edge intersecting_edge) + { + using EK = CGAL::Exact_predicates_exact_constructions_kernel; + const auto to_exact = CGAL::Cartesian_converter(); + + std::cerr << cdt_3_format("restore_subface_region face index: {}, region #{}, intersecting edge #{}: ({} {})\n", + face_index, region_index, edge_index, + IO::oformat(v_above, with_point_and_info), + IO::oformat(v_below, with_point_and_info)); + dump_region(face_index, region_index, this->face_cdt_2[face_index]); + + const auto p_above = this->point(v_above); + const auto p_below = this->point(v_below); + const auto edge_segment = typename Geom_traits::Segment_3{p_above, p_below}; + const auto exact_edge_segment = to_exact(edge_segment); + + std::ofstream intersect_out("dump_edge_region_intersection.xyz"); + intersect_out.precision(17); + for(auto fh: fh_region) { + auto v0 = fh->vertex(0)->info().vertex_handle_3d; + auto v1 = fh->vertex(1)->info().vertex_handle_3d; + auto v2 = fh->vertex(2)->info().vertex_handle_3d; + auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; + auto exact_triangle = to_exact(triangle); + if(auto edge_intersection_opt = CGAL::intersection(exact_edge_segment, exact_triangle)) { + const auto& edge_intersection = *edge_intersection_opt; + if(const auto* p = std::get_if(&edge_intersection)) { + exact(*p); + intersect_out << *p << '\n'; + } + } + } + intersect_out.close(); + + auto cells_around_intersecting_edge = Container_from_circulator{this->incident_cells(intersecting_edge)}; + for(const auto& cell: cells_around_intersecting_edge) { + CGAL_assertion(!cell.has_vertex(tr().infinite_vertex())); + auto tetrahedron = + typename Geom_traits::Tetrahedron_3{tr().point(cell.vertex(0)), tr().point(cell.vertex(1)), + tr().point(cell.vertex(2)), tr().point(cell.vertex(3))}; + for(auto fh: fh_region) { + auto v0 = fh->vertex(0)->info().vertex_handle_3d; + auto v1 = fh->vertex(1)->info().vertex_handle_3d; + auto v2 = fh->vertex(2)->info().vertex_handle_3d; + auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; + auto exact_triangle = to_exact(triangle); + } + + std::cerr << cdt_3_format("Test tetrahedron (#{}):\n {}\n {}\n {}\n {}\n", + cell.time_stamp(), + IO::oformat(cell.vertex(0), with_point_and_info), + IO::oformat(cell.vertex(1), with_point_and_info), + IO::oformat(cell.vertex(2), with_point_and_info), + IO::oformat(cell.vertex(3), with_point_and_info)); + if(!std::any_of(fh_region.begin(), fh_region.end(), [&](const auto fh) { + auto v0 = fh->vertex(0)->info().vertex_handle_3d; + auto v1 = fh->vertex(1)->info().vertex_handle_3d; + auto v2 = fh->vertex(2)->info().vertex_handle_3d; + auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; + bool b = does_tetrahedron_intersect_triangle_interior(tetrahedron, triangle, tr().geom_traits()); + if(b) { + std::cerr << " intersects the region\n"; + } + return b; + })) + { + std::cerr << cdt_3_format( + "ERROR: The following tetrahedron (#{}) does not intersect the region:\n {}\n {}\n {}\n {}\n", + cell.time_stamp(), + IO::oformat(cell.vertex(0), with_point_and_info), IO::oformat(cell.vertex(1), with_point_and_info), + IO::oformat(cell.vertex(2), with_point_and_info), IO::oformat(cell.vertex(3), with_point_and_info)); + } + } + } + + template + void debug_dump_tetrahedra_intersect_region(CDT_3_signed_index face_index, + int region_index, + const CDT_2& cdt_2, + const Fh_region& fh_region) + { + using Mesh = Surface_mesh; + using Face_index = typename Mesh::Face_index; + using EK = CGAL::Exact_predicates_exact_constructions_kernel; + const auto to_exact = CGAL::Cartesian_converter(); + const auto from_exact = CGAL::Cartesian_converter(); + + Mesh tets_intersect_region_mesh; + auto [color_vpmap, _] = tets_intersect_region_mesh.template add_property_map("f:patch_id"); + + for(auto ch : tr().finite_cell_handles()) { + auto tetrahedron = typename Geom_traits::Tetrahedron_3{tr().point(ch->vertex(0)), tr().point(ch->vertex(1)), + tr().point(ch->vertex(2)), tr().point(ch->vertex(3))}; + if(!std::any_of(fh_region.begin(), fh_region.end(), [&](auto fh) { + const auto v0 = fh->vertex(0)->info().vertex_handle_3d; + const auto v1 = fh->vertex(1)->info().vertex_handle_3d; + const auto v2 = fh->vertex(2)->info().vertex_handle_3d; + const auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; + return does_tetrahedron_intersect_triangle_interior(tetrahedron, triangle, tr().geom_traits()); + })) + { + continue; + } + bool intersects = false; + for(int i = 0; i < 4; ++i) { + for(int j = i + 1; j < 4; ++j) { + int intersects_region = does_edge_intersect_region(ch, i, j, cdt_2, fh_region); + if(intersects_region != 0) { + intersects = true; + } + } + } + if(!intersects) { + std::cerr << "ERROR: tetrahedron #" << ch->time_stamp() << " has no edge intersecting the region\n"; + } + std::ofstream dump_tetrahedron( + cdt_3_format("dump_intersecting_{}_{}_tetrahedron_{}.off", face_index, region_index, ch->time_stamp())); + dump_tetrahedron.precision(17); + Mesh mesh; + CGAL::make_tetrahedron(tr().point(ch->vertex(0)), tr().point(ch->vertex(1)), tr().point(ch->vertex(2)), + tr().point(ch->vertex(3)), mesh); + dump_tetrahedron << mesh; + dump_tetrahedron.close(); + + auto exact_tetrahedron = to_exact(tetrahedron); + for(auto fh : fh_region) { + auto v0 = fh->vertex(0)->info().vertex_handle_3d; + auto v1 = fh->vertex(1)->info().vertex_handle_3d; + auto v2 = fh->vertex(2)->info().vertex_handle_3d; + auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; + + auto exact_triangle = to_exact(triangle); + auto tetrahedron_triangle_intersection_opt = CGAL::intersection(exact_tetrahedron, exact_triangle); + if(!tetrahedron_triangle_intersection_opt) { + continue; + } + if(const auto* tri = std::get_if(&tetrahedron_triangle_intersection_opt.value())) { + exact(*tri); + auto v0 = tets_intersect_region_mesh.add_vertex(from_exact((*tri)[0])); + auto v1 = tets_intersect_region_mesh.add_vertex(from_exact((*tri)[1])); + auto v2 = tets_intersect_region_mesh.add_vertex(from_exact((*tri)[2])); + std::array arr{v0, v1, v2}; + auto f = CGAL::Euler::add_face(arr, tets_intersect_region_mesh); + put(color_vpmap, f, static_cast(ch->time_stamp())); + } + if(const auto* vec = std::get_if>(&tetrahedron_triangle_intersection_opt.value())) + { + std::vector vec_of_indices; + for(const auto& p : *vec) { + exact(p); + vec_of_indices.push_back(tets_intersect_region_mesh.add_vertex(from_exact(p))); + } + CGAL::Euler::add_face(vec_of_indices, tets_intersect_region_mesh); + } + } + } + std::ofstream tets_intersect_region_out( + cdt_3_format("dump_tets_intersect_region_{}_{}.ply", face_index, region_index)); + tets_intersect_region_out.precision(17); + CGAL::IO::write_PLY(tets_intersect_region_out, tets_intersect_region_mesh); + tets_intersect_region_out.close(); + } + void dump_3d_triangulation(CDT_3_signed_index face_index, int region_index, std::string type, From 2a8a32d6ad485772640de1532c71c451dbb3e3dd Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 19 Sep 2025 11:53:02 +0200 Subject: [PATCH 03/21] fix Prevent_deref `std::iterator_traits` requires that the iterator type defines all five nested types. I have also modified the implementation of ```c++ make_prevent_deref_range(Range) ``` --- STL_Extension/include/CGAL/iterator.h | 14 ++++++++------ Union_find/include/CGAL/Union_find.h | 3 +++ 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/STL_Extension/include/CGAL/iterator.h b/STL_Extension/include/CGAL/iterator.h index b93190ded73..8f8a2228b07 100644 --- a/STL_Extension/include/CGAL/iterator.h +++ b/STL_Extension/include/CGAL/iterator.h @@ -69,12 +69,6 @@ private: } }; -template -Iterator_range > make_prevent_deref_range(const Iterator_range& range) -{ - return Iterator_range >(make_prevent_deref(range.first), make_prevent_deref(range.second)); -} - template Prevent_deref make_prevent_deref(const I& i) { @@ -87,6 +81,14 @@ Iterator_range > make_prevent_deref_range(const I& begin, const return Iterator_range >(make_prevent_deref(begin), make_prevent_deref(end)); } +template +auto make_prevent_deref_range(const R& range) +{ + using std::begin; + using std::end; + return make_range(make_prevent_deref(begin(range)), make_prevent_deref(end(range))); +} + namespace cpp98 { template #include #include +#include namespace CGAL { @@ -39,6 +40,7 @@ public: typedef R_ reference; typedef P_ pointer; typedef std::forward_iterator_tag iterator_category; + typedef std::ptrdiff_t difference_type; UF_forward_iterator() : m_p(0) {} UF_forward_iterator(PTR_ p) : m_p(p) {} @@ -104,6 +106,7 @@ public: typedef T value_type; typedef T& reference; typedef const T& const_reference; + typedef std::forward_iterator_tag iterator_category; typedef internal::UF_forward_iterator< pointer, T, T&, T*> iterator; typedef iterator handle; From c5ad5bb1135bf16670faf494df47913ed271c9d5 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 19 Sep 2025 12:11:28 +0200 Subject: [PATCH 04/21] more concise code --- ...ing_constrained_Delaunay_triangulation_3.h | 90 +++++++++---------- 1 file changed, 42 insertions(+), 48 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 8e40019e09c..c036dd75bc9 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2030,7 +2030,7 @@ private: } auto test_edge = [&](Cell_handle cell, Vertex_handle v0, int index_v0, Vertex_handle v1, int index_v1, - [[maybe_unused]] int expected) { + int expected) { auto value_returned = [this](bool b) { CGAL_USE(this); if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { @@ -2130,22 +2130,13 @@ private: if(tr().is_infinite(n_ch)) continue; if(new_cell(n_ch)) { - auto tetrahedron = - typename Geom_traits::Tetrahedron_3{tr().point(n_ch->vertex(0)), tr().point(n_ch->vertex(1)), - tr().point(n_ch->vertex(2)), tr().point(n_ch->vertex(3))}; - auto tet_bbox = tetrahedron.bbox(); + const auto tetrahedron = tr().tetrahedron(n_ch); + const auto tet_bbox = tetrahedron.bbox(); if(std::any_of(fh_region.begin(), fh_region.end(), [&](auto fh) { - const auto v0 = fh->vertex(0)->info().vertex_handle_3d; - const auto v1 = fh->vertex(1)->info().vertex_handle_3d; - const auto v2 = fh->vertex(2)->info().vertex_handle_3d; - const auto triangle = - typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; + const auto triangle = cdt_2.triangle(fh); const auto tri_bbox = triangle.bbox(); - if(CGAL::do_overlap(tet_bbox, tri_bbox)) { - return does_tetrahedron_intersect_triangle_interior(tetrahedron, triangle, tr().geom_traits()); - } else { - return false; - } + return CGAL::do_overlap(tet_bbox, tri_bbox) && + does_tetrahedron_intersect_triangle_interior(tetrahedron, triangle, tr().geom_traits()); })) { intersecting_cells.insert(n_ch); @@ -2166,25 +2157,7 @@ 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); - } + process_older_cavity_algorithm(intersecting_edges, intersecting_cells, facets_of_upper_cavity, facets_of_lower_cavity); } // older algorithm std::set facets_of_border; @@ -2322,19 +2295,6 @@ private: } }); } - } // new algorithm - if(this->debug_regions()) { - debug_dump_cavity_outputs(face_index, region_index, intersecting_edges, facets_of_border, facets_of_upper_cavity, facets_of_lower_cavity); - } - - if(this->debug_regions()) { - for(auto edge : intersecting_edges) { - auto [v1, v2] = tr().vertices(edge); - std::cerr << cdt_3_format(" edge: {} {}\n", IO::oformat(v1, with_point_and_info), - IO::oformat(v2, with_point_and_info)); - } - } - if(this->use_newer_cavity_algorithm()) { for(auto facet: facets_of_border) { if(this->debug_regions()) { debug_output_facet_vertices({facet}); @@ -2356,8 +2316,16 @@ private: for(auto v: vertices_of_lower_cavity) { v->ccdt_3_data().clear_mark(Vertex_marker::CAVITY_BELOW); } - } + } // new algorithm + if(this->debug_regions()) { + debug_dump_cavity_outputs(face_index, region_index, intersecting_edges, facets_of_border, facets_of_upper_cavity, facets_of_lower_cavity); + for(auto edge : intersecting_edges) { + auto [v1, v2] = tr().vertices(edge); + std::cerr << cdt_3_format(" edge: {} {}\n", IO::oformat(v1, with_point_and_info), + IO::oformat(v2, with_point_and_info)); + } + } return outputs; } @@ -3603,6 +3571,32 @@ public: write_facets(out, tr, tr.finite_facets()); } + void process_older_cavity_algorithm(const std::vector& intersecting_edges, + const std::set& intersecting_cells, + std::vector& facets_of_upper_cavity, + std::vector& facets_of_lower_cavity) + { + 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); + } + } + void debug_dump_cavity_outputs(CDT_3_signed_index face_index, int region_index, const std::vector& intersecting_edges, From 92fab371296f0efda5c95961a87e9b1cf9788d50 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 19 Sep 2025 16:20:21 +0200 Subject: [PATCH 05/21] add member fcts to handle vertex marks --- ...ing_constrained_Delaunay_triangulation_3.h | 101 ++++++++++-------- 1 file changed, 59 insertions(+), 42 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index c036dd75bc9..84d662f8e16 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -1571,6 +1571,36 @@ public: private: + void set_mark(Vertex_handle v, Vertex_marker m) { + v->ccdt_3_data().set_mark(m); + } + + template + void set_marks(Range_of_vertices&& vertices, Vertex_marker m) { + for(auto v : vertices) { + set_mark(v, m); + } + } + + void clear_mark(Vertex_handle v, Vertex_marker m) { + v->ccdt_3_data().clear_mark(m); + } + + template + void clear_marks(Range_of_vertices&& vertices, Vertex_marker m) { + for(auto v : vertices) { + clear_mark(v, m); + } + } + + bool is_marked(Vertex_handle v, Vertex_marker m) const { + return v->ccdt_3_data().is_marked(m); + } + + bool is_marked(Vertex_handle v) const { + return v->ccdt_3_data().is_marked(); + } + class Non_planar_plc_facet_exception : public std::exception { const char* what() const throw() @@ -1938,8 +1968,8 @@ private: const auto index_vd = this->next_around_edge(index_vb, index_va); //write_segment(dump_edges_around, cell_circ->vertex(index_vc), cell_circ->vertex(index_vd)); - if(cell_circ->vertex(index_vc)->ccdt_3_data().is_marked(Vertex_marker::REGION_BORDER)) continue; - if(cell_circ->vertex(index_vd)->ccdt_3_data().is_marked(Vertex_marker::REGION_BORDER)) continue; + if(is_marked(cell_circ->vertex(index_vc), Vertex_marker::REGION_BORDER)) continue; + if(is_marked(cell_circ->vertex(index_vd), Vertex_marker::REGION_BORDER)) continue; int cd_intersects_region = does_edge_intersect_region(cell_circ, index_vc, index_vd, cdt_2, fh_region); if(cd_intersects_region == 1) { return Search_first_intersection_result_type{ Edge{cell_circ, index_vc, index_vd}, border_edge }; @@ -2044,8 +2074,8 @@ private: } auto [cached_value_it, not_visited] = new_edge(v0, v1, false); if(!not_visited) return value_returned(cached_value_it->second); - int v0v1_intersects_region = (v0->ccdt_3_data().is_marked(Vertex_marker::REGION_INSIDE) || - v1->ccdt_3_data().is_marked(Vertex_marker::REGION_INSIDE)) + int v0v1_intersects_region = (is_marked(v0, Vertex_marker::REGION_INSIDE) || + is_marked(v1, Vertex_marker::REGION_INSIDE)) ? expected : does_edge_intersect_region(cell, index_v0, index_v1, cdt_2, fh_region); if(v0v1_intersects_region != 0) { @@ -2174,8 +2204,8 @@ private: Unique_hash_map::handle> vertices_of_cavity_handles; for(auto c: intersecting_cells) { for(auto v : tr().vertices(c)) { - if(!v->ccdt_3_data().is_marked()) { - v->ccdt_3_data().set_mark(Vertex_marker::CAVITY); + if(!is_marked(v)) { + set_mark(v, Vertex_marker::CAVITY); vertices_of_cavity_handles[v] = vertices_of_cavity_union_find.make_set(v); } } @@ -2185,7 +2215,7 @@ private: for(int i = 0; i < 3; ++i) { auto v1 = vertices[i]; auto v2 = vertices[(i + 1) % 3]; - if(v1->ccdt_3_data().is_marked(Vertex_marker::CAVITY) && v2->ccdt_3_data().is_marked(Vertex_marker::CAVITY)) { + if(is_marked(v1, Vertex_marker::CAVITY) && is_marked(v2, Vertex_marker::CAVITY)) { vertices_of_cavity_union_find.unify_sets(vertices_of_cavity_handles[v1], vertices_of_cavity_handles[v2]); } @@ -2198,8 +2228,7 @@ private: for(int j = i + 1; j < 4; ++j) { const auto v1 = c->vertex(i); const auto v2 = c->vertex(j); - if(v1->ccdt_3_data().is_marked(Vertex_marker::CAVITY) && - v2->ccdt_3_data().is_marked(Vertex_marker::CAVITY) && + if(is_marked(v1, Vertex_marker::CAVITY) && is_marked(v2, Vertex_marker::CAVITY) && non_intersecting_edges_set.count(make_sorted_pair(v1, v2)) > 0) { vertices_of_cavity_union_find.unify_sets(vertices_of_cavity_handles[v1], @@ -2252,13 +2281,13 @@ private: handle != end; ++handle) { auto v = *handle; - v->ccdt_3_data().clear_mark(Vertex_marker::CAVITY); + clear_mark(v, Vertex_marker::CAVITY); if(vertices_of_cavity_union_find.same_set(handle, vertex_above_handle)) { vertices_of_upper_cavity.push_back(v); - v->ccdt_3_data().set_mark(Vertex_marker::CAVITY_ABOVE); + set_mark(v, Vertex_marker::CAVITY_ABOVE); } else if(vertices_of_cavity_union_find.same_set(handle, vertex_below_handle)) { vertices_of_lower_cavity.push_back(v); - v->ccdt_3_data().set_mark(Vertex_marker::CAVITY_BELOW); + set_mark(v, Vertex_marker::CAVITY_BELOW); } else { CGAL_error(); } @@ -2266,7 +2295,7 @@ private: while(std::any_of(intersecting_cells.begin(), intersecting_cells.end(), [&](Cell_handle c) { const auto vs = tr().vertices(c); return std::any_of(vs.begin(), vs.end(), [&](auto v) { - if(!v->ccdt_3_data().is_marked()) { + if(!is_marked(v)) { std::cerr << "INFO: Vertex " << IO::oformat(v, with_point_and_info) << " is not marked\n"; return true; } @@ -2279,16 +2308,16 @@ private: for(int j = i + 1; j < 4; ++j) { auto v1 = c->vertex(i); auto v2 = c->vertex(j); - if(v1->ccdt_3_data().is_marked() != v2->ccdt_3_data().is_marked()) { - if(v2->ccdt_3_data().is_marked()) { + if(is_marked(v1) != is_marked(v2)) { + if(is_marked(v2)) { std::swap(v1, v2); } // here v1 is marked and v2 is not - if(v1->ccdt_3_data().is_marked(Vertex_marker::CAVITY_ABOVE)) { + if(is_marked(v1, Vertex_marker::CAVITY_ABOVE)) { vertices_of_upper_cavity_.push_back(v2); - v2->ccdt_3_data().set_mark(Vertex_marker::CAVITY_ABOVE); - } else if(v1->ccdt_3_data().is_marked(Vertex_marker::CAVITY_BELOW)) { + set_mark(v2, Vertex_marker::CAVITY_ABOVE); + } else if(is_marked(v1, Vertex_marker::CAVITY_BELOW)) { vertices_of_lower_cavity_.push_back(v2); - v2->ccdt_3_data().set_mark(Vertex_marker::CAVITY_BELOW); + set_mark(v2, Vertex_marker::CAVITY_BELOW); } } } @@ -2300,22 +2329,18 @@ private: debug_output_facet_vertices({facet}); } for(auto v: tr().vertices(facet)) { - if(v->ccdt_3_data().is_marked(Vertex_marker::CAVITY_ABOVE)) { + if(is_marked(v, Vertex_marker::CAVITY_ABOVE)) { facets_of_upper_cavity.push_back(facet); break; } - if(v->ccdt_3_data().is_marked(Vertex_marker::CAVITY_BELOW)) { + if(is_marked(v, Vertex_marker::CAVITY_BELOW)) { facets_of_lower_cavity.push_back(facet); break; } } } - for(auto v: vertices_of_upper_cavity) { - v->ccdt_3_data().clear_mark(Vertex_marker::CAVITY_ABOVE); - } - for(auto v: vertices_of_lower_cavity) { - v->ccdt_3_data().clear_mark(Vertex_marker::CAVITY_BELOW); - } + clear_marks(vertices_of_upper_cavity, Vertex_marker::CAVITY_ABOVE); + clear_marks(vertices_of_lower_cavity, Vertex_marker::CAVITY_BELOW); } // new algorithm if(this->debug_regions()) { @@ -2390,13 +2415,9 @@ private: write_segment(out, edge); } } - for(auto v: region_border_vertices) { - v->ccdt_3_data().set_mark(Vertex_marker::REGION_BORDER); - } + set_marks(region_border_vertices, Vertex_marker::REGION_BORDER); const auto found_edge_opt = search_first_intersection(face_index, cdt_2, fh_region, border_edges); - for(auto v: region_border_vertices) { - v->ccdt_3_data().clear_mark(Vertex_marker::REGION_BORDER); - } + clear_marks(region_border_vertices, Vertex_marker::REGION_BORDER); [[maybe_unused]] auto try_flip_region_size_4 = [&] { if(region_border_vertices.size() == 4) { @@ -2534,20 +2555,16 @@ private: } CGAL_assertion(found_edge_opt != std::nullopt); - for(auto v : region_border_vertices) { - v->ccdt_3_data().set_mark(Vertex_marker::REGION_BORDER); - } + set_marks(region_border_vertices, Vertex_marker::REGION_BORDER); for(auto v : region_vertices) { - if(v->ccdt_3_data().is_marked(Vertex_marker::REGION_BORDER)) + if(is_marked(v, Vertex_marker::REGION_BORDER)) continue; - v->ccdt_3_data().set_mark(Vertex_marker::REGION_INSIDE); + set_mark(v, Vertex_marker::REGION_INSIDE); } Scope_exit guard{[&] { - for(auto v : region_vertices) { - v->ccdt_3_data().clear_mark(Vertex_marker::REGION_BORDER); - v->ccdt_3_data().clear_mark(Vertex_marker::REGION_INSIDE); - } + clear_marks(region_vertices, Vertex_marker::REGION_BORDER); + clear_marks(region_vertices, Vertex_marker::REGION_INSIDE); }}; const auto [first_intersecting_edge, _] = *found_edge_opt; From 489a9675dd11218518c44113912c49b3aae4b7de Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 19 Sep 2025 16:39:32 +0200 Subject: [PATCH 06/21] fix iwyu warnings --- TDS_3/include/CGAL/Triangulation_data_structure_3.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/TDS_3/include/CGAL/Triangulation_data_structure_3.h b/TDS_3/include/CGAL/Triangulation_data_structure_3.h index 3c2d36a1b59..7f845200500 100644 --- a/TDS_3/include/CGAL/Triangulation_data_structure_3.h +++ b/TDS_3/include/CGAL/Triangulation_data_structure_3.h @@ -24,8 +24,6 @@ #include #include -#include -#include #include #include #include From 9d08d4087d2894c044394cb17a93aa12cf79437c Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 19 Sep 2025 16:40:12 +0200 Subject: [PATCH 07/21] missing is_marked --- .../CGAL/Conforming_constrained_Delaunay_triangulation_3.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 84d662f8e16..df3b0c64b0b 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2257,7 +2257,7 @@ private: if(facets_of_border.count(Facet{circ, face_index}) > 0) { const auto other_vertex_index = 6 - index_va - index_vb - face_index; const auto other_vertex = circ->vertex(other_vertex_index); - if(other_vertex->ccdt_3_data().is_marked(Vertex_marker::CAVITY)) { + if(is_marked(other_vertex, Vertex_marker::CAVITY)) { vertex_above = circ->vertex(other_vertex_index); break; } From 7267b3a347e138e6ea055346369a0134a70c01b1 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 19 Sep 2025 16:40:58 +0200 Subject: [PATCH 08/21] refactor the definition of vertex_below_handle --- ...ming_constrained_Delaunay_triangulation_3.h | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index df3b0c64b0b..1c4c26be852 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2268,15 +2268,17 @@ private: CGAL_assume(vertex_above != Vertex_handle{}); const auto vertex_above_handle = vertices_of_cavity_handles[vertex_above]; - auto it = vertices_of_cavity_union_find.begin(); - while(it != vertices_of_cavity_union_find.end() && - vertices_of_cavity_union_find.same_set(it, vertex_above_handle)) - { - ++it; - } - CGAL_assertion((it == vertices_of_cavity_union_find.end()) == + + const auto vertex_below_handle = std::invoke([&] { + auto [b, e] = make_prevent_deref_range(vertices_of_cavity_union_find); + return *std::find_if_not( + b, e, [&](auto handle) { return vertices_of_cavity_union_find.same_set(handle, vertex_above_handle); }); + }); + CGAL_assertion(vertex_below_handle == vertices_of_cavity_union_find.end() || + !vertices_of_cavity_union_find.same_set(vertex_below_handle, vertex_above_handle)); + CGAL_assertion((vertex_below_handle == vertices_of_cavity_union_find.end()) == (vertices_of_cavity_union_find.number_of_sets() == 1)); - const auto vertex_below_handle = it; + for(auto handle = vertices_of_cavity_union_find.begin(), end = vertices_of_cavity_union_find.end(); handle != end; ++handle) { From 9726fb59bcc5d89819be66529fc6aeda77cb9ab2 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 23 Sep 2025 09:41:54 +0200 Subject: [PATCH 09/21] const-construct vertex_above --- ...ing_constrained_Delaunay_triangulation_3.h | 49 ++++++++++--------- 1 file changed, 25 insertions(+), 24 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 1c4c26be852..84da0438018 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2240,31 +2240,32 @@ private: break; } } - Vertex_handle vertex_above{}; - Edges_container all_border_edges{border_edges.begin(), border_edges.end()}; - std::for_each(border_edges.begin(), border_edges.end(), [&](auto edge) { - all_border_edges.emplace_back(edge.first, edge.third, edge.second); - }); - for(const auto& border_edge: all_border_edges) { - const auto [border_edge_va, border_edge_vb] = tr().vertices(border_edge); - auto circ = tr().incident_cells(border_edge); - CGAL_assertion(circ != nullptr); - const auto end = circ; - do { - const auto index_va = circ->index(border_edge_va); - const auto index_vb = circ->index(border_edge_vb); - const auto face_index = tr().next_around_edge(index_va, index_vb); - if(facets_of_border.count(Facet{circ, face_index}) > 0) { - const auto other_vertex_index = 6 - index_va - index_vb - face_index; - const auto other_vertex = circ->vertex(other_vertex_index); - if(is_marked(other_vertex, Vertex_marker::CAVITY)) { - vertex_above = circ->vertex(other_vertex_index); - break; + + const Vertex_handle vertex_above = std::invoke([&] { + Edges_container all_border_edges{border_edges.begin(), border_edges.end()}; + std::for_each(border_edges.begin(), border_edges.end(), [&](auto edge) { + all_border_edges.emplace_back(edge.first, edge.third, edge.second); + }); + for(const auto& border_edge: all_border_edges) { + const auto [border_edge_va, border_edge_vb] = tr().vertices(border_edge); + auto circ = tr().incident_cells(border_edge); + CGAL_assertion(circ != nullptr); + const auto end = circ; + do { + const auto index_va = circ->index(border_edge_va); + const auto index_vb = circ->index(border_edge_vb); + const auto face_index = tr().next_around_edge(index_va, index_vb); + if(facets_of_border.count(Facet{circ, face_index}) > 0) { + const auto other_vertex_index = 6 - index_va - index_vb - face_index; + const auto other_vertex = circ->vertex(other_vertex_index); + if(is_marked(other_vertex, Vertex_marker::CAVITY)) { + return circ->vertex(other_vertex_index); + } } - } - } while(++circ != end); - if(vertex_above != Vertex_handle{}) break; - } + } while(++circ != end); + } + return Vertex_handle{}; + }); CGAL_assume(vertex_above != Vertex_handle{}); const auto vertex_above_handle = vertices_of_cavity_handles[vertex_above]; From 9a241b2a9d1581a13f7f3c4f40327f1ebb785e03 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 23 Sep 2025 15:54:04 +0200 Subject: [PATCH 10/21] const-construct facets_of_border --- ...ing_constrained_Delaunay_triangulation_3.h | 22 ++++++++++++------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 84da0438018..9e8fb5f0a05 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2190,18 +2190,24 @@ private: process_older_cavity_algorithm(intersecting_edges, intersecting_cells, facets_of_upper_cavity, facets_of_lower_cavity); } // older algorithm - std::set facets_of_border; - Union_find vertices_of_cavity_union_find; - if(this->use_newer_cavity_algorithm()) { - 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(n, n->index(c)); + // those facets are viewed from the outside of the cavity + const std::set facets_of_border = std::invoke([&] { + std::set facets_of_border; + if(this->use_newer_cavity_algorithm()) { + 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(n, n->index(c)); + } } } } + return facets_of_border; + }); + if(this->use_newer_cavity_algorithm()) { Unique_hash_map::handle> vertices_of_cavity_handles; + Union_find vertices_of_cavity_union_find; for(auto c: intersecting_cells) { for(auto v : tr().vertices(c)) { if(!is_marked(v)) { From 538c2721b341c6a0bcdf2696ed8558728b692756 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 23 Sep 2025 19:06:53 +0200 Subject: [PATCH 11/21] add border_facet_above --- .../Conforming_constrained_Delaunay_triangulation_3.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 9e8fb5f0a05..228e8024d7a 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2247,7 +2247,7 @@ private: } } - const Vertex_handle vertex_above = std::invoke([&] { + const auto [vertex_above, border_facet_above] = std::invoke([&] { Edges_container all_border_edges{border_edges.begin(), border_edges.end()}; std::for_each(border_edges.begin(), border_edges.end(), [&](auto edge) { all_border_edges.emplace_back(edge.first, edge.third, edge.second); @@ -2261,16 +2261,17 @@ private: const auto index_va = circ->index(border_edge_va); const auto index_vb = circ->index(border_edge_vb); const auto face_index = tr().next_around_edge(index_va, index_vb); - if(facets_of_border.count(Facet{circ, face_index}) > 0) { + Facet facet{circ, face_index}; + if(facets_of_border.count(facet) > 0) { const auto other_vertex_index = 6 - index_va - index_vb - face_index; const auto other_vertex = circ->vertex(other_vertex_index); if(is_marked(other_vertex, Vertex_marker::CAVITY)) { - return circ->vertex(other_vertex_index); + return std::make_pair(circ->vertex(other_vertex_index), facet); } } } while(++circ != end); } - return Vertex_handle{}; + return std::make_pair(Vertex_handle{}, Facet{}); }); CGAL_assume(vertex_above != Vertex_handle{}); From 5e80ca60bb7af7e6fd6c91c6d1840c7d75a7ccfc Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 23 Sep 2025 21:37:55 +0200 Subject: [PATCH 12/21] fix when CGAL_CDT_3_CAN_USE_CXX20_FORMAT is false --- .../include/CGAL/Constrained_triangulation_3/internal/config.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Constrained_triangulation_3/include/CGAL/Constrained_triangulation_3/internal/config.h b/Constrained_triangulation_3/include/CGAL/Constrained_triangulation_3/internal/config.h index 99c39f2555b..47f53d42ac5 100644 --- a/Constrained_triangulation_3/include/CGAL/Constrained_triangulation_3/internal/config.h +++ b/Constrained_triangulation_3/include/CGAL/Constrained_triangulation_3/internal/config.h @@ -48,7 +48,7 @@ decltype(auto) cdt_3_format(std::string_view fmt, const Args&... args) { template constexpr decltype(auto) cdt_3_format(Args&&...) { - return ""; + return std::string{}; } constexpr bool cdt_3_can_use_cxx20_format() { From f6ebe208e9352843662f96336bed1e84b2bdac11 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 24 Sep 2025 16:59:32 +0200 Subject: [PATCH 13/21] extract a member function `detect_edges_and_cells_intersecting_region` --- ...ing_constrained_Delaunay_triangulation_3.h | 335 +++++++++--------- 1 file changed, 177 insertions(+), 158 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 228e8024d7a..7e900f2c954 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2020,174 +2020,21 @@ private: facets_of_upper_cavity, facets_of_lower_cavity] = outputs; // to avoid "warning: captured structured bindings are a C++20 extension [-Wc++20-extensions]"" - auto& intersecting_edges_ = intersecting_edges; auto& vertices_of_upper_cavity_ = vertices_of_upper_cavity; auto& vertices_of_lower_cavity_ = vertices_of_lower_cavity; std::set> non_intersecting_edges_set; - // marker for already visited elements - std::set visited_vertices; - std::map, bool> visited_edges; - std::set visited_cells; - - auto make__new_element_functor = [](auto& visited_set) { - return [&visited_set](auto... e) { - const auto [_, not_already_visited] = visited_set.emplace(e...); - return not_already_visited; - }; - }; - - 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) { - CGAL_assertion(v0 != Vertex_handle{}); - return visited_edges.emplace(CGAL::make_sorted_pair(v0, v1), does_intersect); - }; - if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { debug_dump_tetrahedra_intersect_region(face_index, region_index, cdt_2, fh_region); } - intersecting_edges.push_back(first_intersecting_edge); - const auto [v0, v1] = tr().vertices(first_intersecting_edge); - (void)new_edge(v0, v1, true); - for(std::size_t i = 0; i < intersecting_edges.size(); ++i) { - const auto intersecting_edge = intersecting_edges[i]; - const auto [v_above, v_below] = tr().vertices(intersecting_edge); - if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { - debug_dump_edge_region_intersection(face_index, region_index, fh_region, i, v_above, v_below, intersecting_edge); - } - - auto test_edge = [&](Cell_handle cell, Vertex_handle v0, int index_v0, Vertex_handle v1, int index_v1, - int expected) { - auto value_returned = [this](bool b) { - CGAL_USE(this); - if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { - std::cerr << cdt_3_format(" return {}\n", b); - } - return b; - }; - if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { - std::cerr << cdt_3_format(" test_edge {} {} ", IO::oformat(v0, with_point_and_info), - IO::oformat(v1, with_point_and_info)); - } - auto [cached_value_it, not_visited] = new_edge(v0, v1, false); - if(!not_visited) return value_returned(cached_value_it->second); - int v0v1_intersects_region = (is_marked(v0, Vertex_marker::REGION_INSIDE) || - is_marked(v1, Vertex_marker::REGION_INSIDE)) - ? expected - : does_edge_intersect_region(cell, index_v0, index_v1, cdt_2, fh_region); - if(v0v1_intersects_region != 0) { - if(this->use_older_cavity_algorithm()) { - if(v0v1_intersects_region != expected) { - 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); - } - intersecting_edges_.emplace_back(cell, index_v0, index_v1); - cached_value_it->second = true; - return value_returned(true); - } else { - non_intersecting_edges_set.insert(make_sorted_pair(v0, v1)); - cached_value_it->second = false; - return value_returned(false); - } - }; - - 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] - CGAL_assertion(false == this->is_infinite(*facet_circ)); - const auto cell = facet_circ->first; - const auto facet_index = facet_circ->second; - CGAL_assertion_msg(!cell->ccdt_3_data().is_facet_constrained(facet_index), - std::invoke([&]() { - this->dump_triangulation_to_off(); - return std::string("intersecting polygons!"); - }).c_str()); - if(new_cell(cell)) { - intersecting_cells.insert(cell); - } - const auto index_v_above = cell->index(v_above); - const auto index_v_below = cell->index(v_below); - const auto index_vc = 6 - index_v_above - index_v_below - facet_index; - const auto vc = cell->vertex(index_vc); - 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) && - this->use_older_cavity_algorithm()) - { - dump_triangulation(); - dump_region(face_index, region_index, cdt_2); - { - std::ofstream out(std::string("dump_two_edges_") + std::to_string(face_index) + ".polylines.txt"); - out.precision(17); - write_segment(out, Edge{cell, index_v_above, index_vc}); - write_segment(out, Edge{cell, index_v_below, index_vc}); - } - throw PLC_error{"PLC error: !test_edge(v_above..) && !test_edge(v_below..)" , - __FILE__, __LINE__, face_index, region_index}; - } - } while(++facet_circ != facet_circ_end); - if(this->use_newer_cavity_algorithm() && i + 1 == intersecting_edges.size()) { - for(auto ch: intersecting_cells) { - if(this->debug_regions()) { - std::cerr << "tetrahedron #" << ch->time_stamp() << " intersects the region\n"; - } - for(int i = 0; i < 4; ++i) { - for(int j = i + 1; j < 4; ++j) { - test_edge(ch, ch->vertex(i), i, ch->vertex(j), j, 1); - } - } - for(int i = 0; i < 4; ++i) { - auto n_ch = ch->neighbor(i); - if(tr().is_infinite(n_ch)) - continue; - if(new_cell(n_ch)) { - const auto tetrahedron = tr().tetrahedron(n_ch); - const auto tet_bbox = tetrahedron.bbox(); - if(std::any_of(fh_region.begin(), fh_region.end(), [&](auto fh) { - const auto triangle = cdt_2.triangle(fh); - const auto tri_bbox = triangle.bbox(); - return CGAL::do_overlap(tet_bbox, tri_bbox) && - does_tetrahedron_intersect_triangle_interior(tetrahedron, triangle, tr().geom_traits()); - })) - { - intersecting_cells.insert(n_ch); - if(this->debug_regions()) { - std::cerr << "tetrahedron #" << n_ch->time_stamp() << " intersects the region\n"; - } - } else if(this->debug_regions()) { - std::cerr << "NO, tetrahedron #" << n_ch->time_stamp() << " does not intersect the region\n"; - } - for(int i = 0; i < 4; ++i) { - for(int j = i + 1; j < 4; ++j) { - test_edge(n_ch, n_ch->vertex(i), i, n_ch->vertex(j), j, 1); - } - } - } - } - } - } // last intersecting edge, and new algorithm - } // end loop on intersecting_edges + detect_edges_and_cells_intersecting_region(face_index, region_index, cdt_2, fh_region, region_border_vertices, + first_intersecting_edge, intersecting_edges, intersecting_cells, + non_intersecting_edges_set); if(this->use_older_cavity_algorithm()) { - process_older_cavity_algorithm(intersecting_edges, intersecting_cells, facets_of_upper_cavity, facets_of_lower_cavity); + process_older_cavity_algorithm(intersecting_edges, intersecting_cells, vertices_of_upper_cavity, + vertices_of_lower_cavity, facets_of_upper_cavity, facets_of_lower_cavity); } // older algorithm // those facets are viewed from the outside of the cavity @@ -3600,12 +3447,23 @@ public: void process_older_cavity_algorithm(const std::vector& intersecting_edges, const std::set& intersecting_cells, + std::vector& vertices_of_upper_cavity, + std::vector& vertices_of_lower_cavity, std::vector& facets_of_upper_cavity, std::vector& facets_of_lower_cavity) { + auto new_vertex = make__new_element_functor(); + for(auto intersecting_edge: intersecting_edges) { const auto [v_above, v_below] = tr().vertices(intersecting_edge); + 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 cell_circ = this->incident_cells(intersecting_edge), end = cell_circ; CGAL_assume(cell_circ != nullptr); do { @@ -3624,6 +3482,167 @@ public: } } + template + static auto make__new_element_functor() { + return [visited_set = std::set()](auto... e) mutable { + const auto [_, not_already_visited] = visited_set.emplace(e...); + return not_already_visited; + }; + }; + + template + void detect_edges_and_cells_intersecting_region( + CDT_3_signed_index face_index, + int region_index, + const CDT_2& cdt_2, + const Fh_region& fh_region, + const Vertices_container& region_border_vertices, + Edge first_intersecting_edge, + std::vector& intersecting_edges, + std::set& intersecting_cells, + std::set>& non_intersecting_edges_set) + { + // Create visitor functors + auto new_cell = make__new_element_functor(); + std::map, bool> visited_edges; + auto new_edge = [&](Vertex_handle v0, Vertex_handle v1, bool does_intersect) { + CGAL_assertion(v0 != Vertex_handle{}); + return visited_edges.emplace(CGAL::make_sorted_pair(v0, v1), does_intersect); + }; + + // Alias for C++20 extension warning workaround + auto& intersecting_edges_ = intersecting_edges; + + intersecting_edges.push_back(first_intersecting_edge); + const auto [v0, v1] = tr().vertices(first_intersecting_edge); + (void)new_edge(v0, v1, true); + + for(std::size_t i = 0; i < intersecting_edges.size(); ++i) { + const auto intersecting_edge = intersecting_edges[i]; + const auto [v_above, v_below] = tr().vertices(intersecting_edge); + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { + debug_dump_edge_region_intersection(face_index, region_index, fh_region, i, v_above, v_below, intersecting_edge); + } + + auto test_edge = [&](Cell_handle cell, Vertex_handle v0, int index_v0, Vertex_handle v1, int index_v1, + int expected) { + auto value_returned = [this](bool b) { + CGAL_USE(this); + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { + std::cerr << cdt_3_format(" return {}\n", b); + } + return b; + }; + if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { + std::cerr << cdt_3_format(" test_edge {} {} ", IO::oformat(v0, with_point_and_info), + IO::oformat(v1, with_point_and_info)); + } + auto [cached_value_it, not_visited] = new_edge(v0, v1, false); + if(!not_visited) return value_returned(cached_value_it->second); + int v0v1_intersects_region = (is_marked(v0, Vertex_marker::REGION_INSIDE) || + is_marked(v1, Vertex_marker::REGION_INSIDE)) + ? expected + : does_edge_intersect_region(cell, index_v0, index_v1, cdt_2, fh_region); + if(v0v1_intersects_region != 0) { + if(this->use_older_cavity_algorithm()) { + if(v0v1_intersects_region != expected) { + 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); + } + intersecting_edges_.emplace_back(cell, index_v0, index_v1); + cached_value_it->second = true; + return value_returned(true); + } else { + non_intersecting_edges_set.insert(make_sorted_pair(v0, v1)); + cached_value_it->second = false; + return value_returned(false); + } + }; + + auto facet_circ = this->incident_facets(intersecting_edge); + const auto facet_circ_end = facet_circ; + do { // loop facets around [v_above, v_below] + CGAL_assertion(false == this->is_infinite(*facet_circ)); + const auto cell = facet_circ->first; + const auto facet_index = facet_circ->second; + CGAL_assertion_msg(!cell->ccdt_3_data().is_facet_constrained(facet_index), + std::invoke([&]() { + this->dump_triangulation_to_off(); + return std::string("intersecting polygons!"); + }).c_str()); + if(new_cell(cell)) { + intersecting_cells.insert(cell); + } + const auto index_v_above = cell->index(v_above); + const auto index_v_below = cell->index(v_below); + const auto index_vc = 6 - index_v_above - index_v_below - facet_index; + const auto vc = cell->vertex(index_vc); + 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) && + this->use_older_cavity_algorithm()) + { + if(this->debug_regions()) { + dump_triangulation(); + dump_region(face_index, region_index, cdt_2); + std::ofstream out(std::string("dump_two_edges_") + std::to_string(face_index) + ".polylines.txt"); + out.precision(17); + write_segment(out, Edge{cell, index_v_above, index_vc}); + write_segment(out, Edge{cell, index_v_below, index_vc}); + } + throw PLC_error{"PLC error: !test_edge(v_above..) && !test_edge(v_below..)" , + __FILE__, __LINE__, face_index, region_index}; + } + } while(++facet_circ != facet_circ_end); + if(this->use_newer_cavity_algorithm() && i + 1 == intersecting_edges.size()) { + for(auto ch: intersecting_cells) { + if(this->debug_regions()) { + std::cerr << "tetrahedron #" << ch->time_stamp() << " intersects the region\n"; + } + for(int i = 0; i < 4; ++i) { + for(int j = i + 1; j < 4; ++j) { + test_edge(ch, ch->vertex(i), i, ch->vertex(j), j, 1); + } + } + for(int i = 0; i < 4; ++i) { + auto n_ch = ch->neighbor(i); + if(tr().is_infinite(n_ch)) + continue; + if(new_cell(n_ch)) { + const auto tetrahedron = tr().tetrahedron(n_ch); + const auto tet_bbox = tetrahedron.bbox(); + if(std::any_of(fh_region.begin(), fh_region.end(), [&](auto fh) { + const auto triangle = cdt_2.triangle(fh); + const auto tri_bbox = triangle.bbox(); + return CGAL::do_overlap(tet_bbox, tri_bbox) && + does_tetrahedron_intersect_triangle_interior(tetrahedron, triangle, tr().geom_traits()); + })) + { + intersecting_cells.insert(n_ch); + if(this->debug_regions()) { + std::cerr << "tetrahedron #" << n_ch->time_stamp() << " intersects the region\n"; + } + } else if(this->debug_regions()) { + std::cerr << "NO, tetrahedron #" << n_ch->time_stamp() << " does not intersect the region\n"; + } + for(int i = 0; i < 4; ++i) { + for(int j = i + 1; j < 4; ++j) { + test_edge(n_ch, n_ch->vertex(i), i, n_ch->vertex(j), j, 1); + } + } + } + } + } + } // last intersecting edge, and new algorithm + } // end loop on intersecting_edges + } + void debug_dump_cavity_outputs(CDT_3_signed_index face_index, int region_index, const std::vector& intersecting_edges, From 1a2c0613646e88f09c2f45c08a17769d296ea71f Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 24 Sep 2025 17:10:26 +0200 Subject: [PATCH 14/21] add another way to reduce to two sets with union-find --- ...ing_constrained_Delaunay_triangulation_3.h | 81 +++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 7e900f2c954..add29daf565 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -34,6 +34,8 @@ #include #include #include +#include +#include #include #include @@ -74,6 +76,7 @@ #include #include #include +#include #include #if CGAL_CXX20 && __has_include() # include @@ -2122,6 +2125,84 @@ private: }); CGAL_assume(vertex_above != Vertex_handle{}); + if(vertices_of_cavity_union_find.number_of_sets() > 2) { + const auto border_edges_set = std::invoke([&] { + using Ordered_pair = std::pair; + using Hash = boost::hash; + CGAL::unordered_flat_set border_edges_set; + for(auto edge : border_edges) { + auto [va, vb] = tr().vertices(edge); + border_edges_set.insert(make_sorted_pair(va, vb)); + } + return border_edges_set; + }); + CGAL::unordered_flat_set> remaining_facets_of_border(facets_of_border.begin(), + facets_of_border.end()); + + std::stack stack; + std::optional::handle> reference_handle_of_the_connected_component; + stack.push(border_facet_above); + remaining_facets_of_border.erase(border_facet_above); + while(!stack.empty()) { + const auto facet = stack.top(); + stack.pop(); + const auto [cell, facet_index] = facet; // border facet seen from the outside of the cavity + CGAL_assertion(intersecting_cells.count(cell) == 0); //REMOVE + CGAL_assertion(intersecting_cells.count(cell->neighbor(facet_index)) > 0); //REMOVE + const auto vertices = tr().vertices(facet); + for(auto v : vertices) { + if(is_marked(v, Vertex_marker::CAVITY)) { + if(!reference_handle_of_the_connected_component) { + reference_handle_of_the_connected_component = vertices_of_cavity_handles[v]; + } else { + vertices_of_cavity_union_find.unify_sets(*reference_handle_of_the_connected_component, + vertices_of_cavity_handles[v]); + } + } + } + for(int i = 0; i < 3; ++i) { + const auto va = vertices[i]; + const auto vb = vertices[tr().ccw(i)]; + + if((is_marked(va, Vertex_marker::CAVITY) || is_marked(vb, Vertex_marker::CAVITY)) && + border_edges_set.count(make_sorted_pair(va, vb)) == 0) + { + // loop around the edge [va, vb] to get another facet on the border of the cavity + auto previous_cell = cell; + auto other_cell = cell->neighbor(facet_index); + do { + CGAL_assertion(intersecting_cells.count(other_cell) >= 0); // REMOVE + auto index_va = other_cell->index(va); + auto index_vb = other_cell->index(vb); + auto other_facet_index = tr().next_around_edge(index_vb, index_va); + previous_cell = other_cell; + other_cell = previous_cell->neighbor(other_facet_index); + + } while(intersecting_cells.count(other_cell) > 0); + const Facet neighbor_facet{other_cell, other_cell->index(previous_cell)}; + CGAL_assertion(facets_of_border.count(neighbor_facet) > 0); + if(remaining_facets_of_border.erase(neighbor_facet) > 0) { + stack.push(neighbor_facet); + } + } + } + + // if the stack is empty but there are still facets to process, we start again to recover + // another connected component of the cavity border + if(stack.empty() && !remaining_facets_of_border.empty()) { + stack.push(*remaining_facets_of_border.begin()); + remaining_facets_of_border.erase(remaining_facets_of_border.begin()); + reference_handle_of_the_connected_component.reset(); + } + } + } + CGAL_assertion_msg(vertices_of_cavity_union_find.number_of_sets() <= 2, + std::invoke([&] { + std::stringstream ss; + ss << "Error: cavity has " << vertices_of_cavity_union_find.number_of_sets() + << " sub-cavities (should be <=2)\n"; + return ss.str(); + }).c_str()); const auto vertex_above_handle = vertices_of_cavity_handles[vertex_above]; const auto vertex_below_handle = std::invoke([&] { From 44e3994c5180204d63f3323139114e9912acaccb Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 24 Sep 2025 18:25:59 +0200 Subject: [PATCH 15/21] careful when the union-find has only one set --- ...rming_constrained_Delaunay_triangulation_3.h | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index add29daf565..a0813ae75e2 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2056,8 +2056,10 @@ private: return facets_of_border; }); if(this->use_newer_cavity_algorithm()) { - Unique_hash_map::handle> vertices_of_cavity_handles; + Union_find vertices_of_cavity_union_find; + using Union_find_handle = typename Union_find::handle; + Unique_hash_map vertices_of_cavity_handles; for(auto c: intersecting_cells) { for(auto v : tr().vertices(c)) { if(!is_marked(v)) { @@ -2140,7 +2142,7 @@ private: facets_of_border.end()); std::stack stack; - std::optional::handle> reference_handle_of_the_connected_component; + std::optional reference_handle_of_the_connected_component; stack.push(border_facet_above); remaining_facets_of_border.erase(border_facet_above); while(!stack.empty()) { @@ -2207,8 +2209,13 @@ private: const auto vertex_below_handle = std::invoke([&] { auto [b, e] = make_prevent_deref_range(vertices_of_cavity_union_find); - return *std::find_if_not( + auto it = std::find_if_not( b, e, [&](auto handle) { return vertices_of_cavity_union_find.same_set(handle, vertex_above_handle); }); + if(it != e) { + return *it; + } else { + return std::as_const(vertices_of_cavity_union_find).end(); + } }); CGAL_assertion(vertex_below_handle == vertices_of_cavity_union_find.end() || !vertices_of_cavity_union_find.same_set(vertex_below_handle, vertex_above_handle)); @@ -2223,7 +2230,9 @@ private: if(vertices_of_cavity_union_find.same_set(handle, vertex_above_handle)) { vertices_of_upper_cavity.push_back(v); set_mark(v, Vertex_marker::CAVITY_ABOVE); - } else if(vertices_of_cavity_union_find.same_set(handle, vertex_below_handle)) { + } else if(vertex_below_handle != vertices_of_cavity_union_find.end() && + vertices_of_cavity_union_find.same_set(handle, vertex_below_handle)) + { vertices_of_lower_cavity.push_back(v); set_mark(v, Vertex_marker::CAVITY_BELOW); } else { From 3f088282c5437a1847ec39fb0b93b4d8f6babf35 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 24 Sep 2025 17:42:08 +0200 Subject: [PATCH 16/21] add comments add comments --- ...ing_constrained_Delaunay_triangulation_3.h | 20 +++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index a0813ae75e2..729fd9d4fc1 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2057,6 +2057,7 @@ private: }); if(this->use_newer_cavity_algorithm()) { + // create a union-find of the vertices of the cavity (but those on the region border) Union_find vertices_of_cavity_union_find; using Union_find_handle = typename Union_find::handle; Unique_hash_map vertices_of_cavity_handles; @@ -2068,6 +2069,8 @@ private: } } } + + // use edges of the border facets to unify sets for(auto facet: facets_of_border) { auto vertices = tr().vertices(facet); for(int i = 0; i < 3; ++i) { @@ -2079,6 +2082,8 @@ private: } } } + + // use non-intersecting edges to unify sets, until we have at most 2 sets if(vertices_of_cavity_union_find.number_of_sets() > 2) { for(auto c : intersecting_cells) { @@ -2099,6 +2104,7 @@ private: } } + // find a vertex and a border-facet above the region const auto [vertex_above, border_facet_above] = std::invoke([&] { Edges_container all_border_edges{border_edges.begin(), border_edges.end()}; std::for_each(border_edges.begin(), border_edges.end(), [&](auto edge) { @@ -2127,6 +2133,8 @@ private: }); CGAL_assume(vertex_above != Vertex_handle{}); + // if there are still more than 2 sets, we need to propagate the information + // using a DFS on the border facets if(vertices_of_cavity_union_find.number_of_sets() > 2) { const auto border_edges_set = std::invoke([&] { using Ordered_pair = std::pair; @@ -2198,6 +2206,7 @@ private: } } } + CGAL_assertion_msg(vertices_of_cavity_union_find.number_of_sets() <= 2, std::invoke([&] { std::stringstream ss; @@ -2207,6 +2216,7 @@ private: }).c_str()); const auto vertex_above_handle = vertices_of_cavity_handles[vertex_above]; + // find a vertex below the region (if any) const auto vertex_below_handle = std::invoke([&] { auto [b, e] = make_prevent_deref_range(vertices_of_cavity_union_find); auto it = std::find_if_not( @@ -2222,8 +2232,8 @@ private: CGAL_assertion((vertex_below_handle == vertices_of_cavity_union_find.end()) == (vertices_of_cavity_union_find.number_of_sets() == 1)); - for(auto handle = vertices_of_cavity_union_find.begin(), end = vertices_of_cavity_union_find.end(); - handle != end; ++handle) + // use the union-find sets to mark vertices as above or below the region + for(auto handle : make_prevent_deref_range(vertices_of_cavity_union_find)) { auto v = *handle; clear_mark(v, Vertex_marker::CAVITY); @@ -2239,6 +2249,10 @@ private: CGAL_error(); } } + + // if any vertex is still unmarked, it means that the union-find did not + // connect all the vertices of the cavity. Then propagate the information + // using the intersecting cells. while(std::any_of(intersecting_cells.begin(), intersecting_cells.end(), [&](Cell_handle c) { const auto vs = tr().vertices(c); return std::any_of(vs.begin(), vs.end(), [&](auto v) { @@ -2271,6 +2285,8 @@ private: } }); } + + // classify the facets of the border of the cavity for(auto facet: facets_of_border) { if(this->debug_regions()) { debug_output_facet_vertices({facet}); From 88f9f009226546f61618fae9bc4ab2f97ec1d217 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 25 Sep 2025 16:25:37 +0200 Subject: [PATCH 17/21] fix warning -Wstringop-overflow MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fix that warning, due to a copy of `tuple` in compare operators for `std::sort` and `std::unique`, in `CGAL::Polygon_mesh_processing::autorefine_impl::collect_intersections`. ``` In member function ‘std::__atomic_base<_IntTp>::__int_type std::__atomic_base<_IntTp>::fetch_add(__int_type, std::memory_order) [with _ITp = int]’, inlined from ‘void CGAL::Handle::incref() const’ at /mnt/testsuite/include/CGAL/Handle.h:87:29, inlined from ‘CGAL::Handle::Handle(const CGAL::Handle&)’ at /mnt/testsuite/include/CGAL/Handle.h:55:13, inlined from ‘CGAL::Lazy > >, CGAL::Point_3 > >, CGAL::Cartesian_converter >, CGAL::Simple_cartesian >, CGAL::NT_converter<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, CGAL::Interval_nt > > >::Lazy(const CGAL::Lazy > >, CGAL::Point_3 > >, CGAL::Cartesian_converter >, CGAL::Simple_cartesian >, CGAL::NT_converter<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, CGAL::Interval_nt > > >&)’ at /mnt/testsuite/include/CGAL/Lazy.h:877:7, inlined from ‘CGAL::Point_3::Point_3(const CGAL::Point_3&)’ at /mnt/testsuite/include/CGAL/Point_3.h:30:7, inlined from ‘std::_Head_base<_Idx, _Head, false>::_Head_base(const std::_Head_base<_Idx, _Head, false>&) [with long unsigned int _Idx = 0; _Head = CGAL::Point_3]’ at /usr/include/c++/15/tuple:208:17, inlined from ‘std::_Tuple_impl<_Idx, _Head, _Tail ...>::_Tuple_impl(const std::_Tuple_impl<_Idx, _Head, _Tail ...>&) [with long unsigned int _Idx = 0; _Head = CGAL::Point_3; _Tail = {int, int}]’ at /usr/include/c++/15/tuple:318:17, inlined from ‘std::tuple< >::tuple(const std::tuple< >&) [with _Elements = {CGAL::Point_3, int, int}]’ at /usr/include/c++/15/tuple:1502:17, inlined from ‘bool __gnu_cxx::__ops::_Val_comp_iter<_Compare>::operator()(_Value&, _Iterator) [with _Value = std::tuple, int, int>; _Iterator = __gnu_cxx::__normal_iterator, int, int>*, std::vector, int, int> > >; _Compare = CGAL::Polygon_mesh_processing::autorefine_impl::collect_intersections(const std::array, 3>&, const std::array, 3>&, std::vector, int, int> >&)::]’ at /usr/include/c++/15/bits/predefined_ops.h:240:23, inlined from ‘void std::__unguarded_linear_insert(_RandomAccessIterator, _Compare) [with _RandomAccessIterator = __gnu_cxx::__normal_iterator, int, int>*, vector, int, int> > >; _Compare = __gnu_cxx::__ops::_Val_comp_iter(const std::array, 3>&, const std::array, 3>&, std::vector, int, int> >&):: >]’ at /usr/include/c++/15/bits/stl_algo.h:1758:20: /usr/include/c++/15/bits/atomic_base.h:631:34: warning: ‘unsigned int __atomic_fetch_add_4(volatile void*, unsigned int, int)’ writing 4 bytes into a region of size 0 overflows the destination [-Wstringop-overflow=] 631 | { return __atomic_fetch_add(&_M_i, __i, int(__m)); } | ~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~ ``` See for example https://cgal.geometryfactory.com/CGAL/testsuite/CGAL-6.2-Ic-6/Constrained_triangulation_3_Examples/TestReport_cgaltest_Fedora-rawhide-Release.gz --- .../include/CGAL/Polygon_mesh_processing/autorefinement.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index 54fce583918..9432e29a421 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -475,8 +475,8 @@ bool collect_intersections(const std::array& t1, // #warning TODO get rid of sort and unique calls // because we don't handle intersection type and can have edge-edge edge-vertex duplicates - std::sort(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return std::get<0>(p)(q);}); - auto last = std::unique(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return std::get<0>(p)==std::get<0>(q);}); + std::sort(inter_pts.begin(), inter_pts.end(), [](const auto& p, const auto& q){return std::get<0>(p)(q);}); + auto last = std::unique(inter_pts.begin(), inter_pts.end(), [](const auto& p, const auto& q){return std::get<0>(p)==std::get<0>(q);}); inter_pts.erase(last, inter_pts.end()); #ifdef CGAL_AUTOREF_DEBUG_DEPTH From 7df5c3a7bb3285ec597391798fcc4580580797b1 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 25 Sep 2025 16:26:13 +0200 Subject: [PATCH 18/21] fix warning: captured structured bindings are a C++20 extension [-Wc++20-extensions] --- ...ing_constrained_Delaunay_triangulation_3.h | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 729fd9d4fc1..61c398d6448 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2025,6 +2025,7 @@ private: // to avoid "warning: captured structured bindings are a C++20 extension [-Wc++20-extensions]"" auto& vertices_of_upper_cavity_ = vertices_of_upper_cavity; auto& vertices_of_lower_cavity_ = vertices_of_lower_cavity; + const auto& cr_intersecting_cells = intersecting_cells; std::set> non_intersecting_edges_set; @@ -2036,7 +2037,7 @@ private: first_intersecting_edge, intersecting_edges, intersecting_cells, non_intersecting_edges_set); if(this->use_older_cavity_algorithm()) { - process_older_cavity_algorithm(intersecting_edges, intersecting_cells, vertices_of_upper_cavity, + process_older_cavity_algorithm(intersecting_edges, cr_intersecting_cells, vertices_of_upper_cavity, vertices_of_lower_cavity, facets_of_upper_cavity, facets_of_lower_cavity); } // older algorithm @@ -2044,10 +2045,10 @@ private: const std::set facets_of_border = std::invoke([&] { std::set facets_of_border; if(this->use_newer_cavity_algorithm()) { - for(auto c : intersecting_cells) { + for(auto c : cr_intersecting_cells) { for(int i = 0; i < 4; ++i) { auto n = c->neighbor(i); - if(intersecting_cells.count(n) == 0) { + if(cr_intersecting_cells.count(n) == 0) { facets_of_border.emplace(n, n->index(c)); } } @@ -2061,7 +2062,7 @@ private: Union_find vertices_of_cavity_union_find; using Union_find_handle = typename Union_find::handle; Unique_hash_map vertices_of_cavity_handles; - for(auto c: intersecting_cells) { + for(auto c: cr_intersecting_cells) { for(auto v : tr().vertices(c)) { if(!is_marked(v)) { set_mark(v, Vertex_marker::CAVITY); @@ -2085,7 +2086,7 @@ private: // use non-intersecting edges to unify sets, until we have at most 2 sets if(vertices_of_cavity_union_find.number_of_sets() > 2) { - for(auto c : intersecting_cells) { + for(auto c : cr_intersecting_cells) { for(int i = 0; i < 4; ++i) { for(int j = i + 1; j < 4; ++j) { @@ -2157,8 +2158,8 @@ private: const auto facet = stack.top(); stack.pop(); const auto [cell, facet_index] = facet; // border facet seen from the outside of the cavity - CGAL_assertion(intersecting_cells.count(cell) == 0); //REMOVE - CGAL_assertion(intersecting_cells.count(cell->neighbor(facet_index)) > 0); //REMOVE + CGAL_assertion(cr_intersecting_cells.count(cell) == 0); //REMOVE + CGAL_assertion(cr_intersecting_cells.count(cell->neighbor(facet_index)) > 0); //REMOVE const auto vertices = tr().vertices(facet); for(auto v : vertices) { if(is_marked(v, Vertex_marker::CAVITY)) { @@ -2181,14 +2182,14 @@ private: auto previous_cell = cell; auto other_cell = cell->neighbor(facet_index); do { - CGAL_assertion(intersecting_cells.count(other_cell) >= 0); // REMOVE + CGAL_assertion(cr_intersecting_cells.count(other_cell) >= 0); // REMOVE auto index_va = other_cell->index(va); auto index_vb = other_cell->index(vb); auto other_facet_index = tr().next_around_edge(index_vb, index_va); previous_cell = other_cell; other_cell = previous_cell->neighbor(other_facet_index); - } while(intersecting_cells.count(other_cell) > 0); + } while(cr_intersecting_cells.count(other_cell) > 0); const Facet neighbor_facet{other_cell, other_cell->index(previous_cell)}; CGAL_assertion(facets_of_border.count(neighbor_facet) > 0); if(remaining_facets_of_border.erase(neighbor_facet) > 0) { @@ -2253,7 +2254,7 @@ private: // if any vertex is still unmarked, it means that the union-find did not // connect all the vertices of the cavity. Then propagate the information // using the intersecting cells. - while(std::any_of(intersecting_cells.begin(), intersecting_cells.end(), [&](Cell_handle c) { + while(std::any_of(cr_intersecting_cells.begin(), cr_intersecting_cells.end(), [&](Cell_handle c) { const auto vs = tr().vertices(c); return std::any_of(vs.begin(), vs.end(), [&](auto v) { if(!is_marked(v)) { @@ -2264,7 +2265,7 @@ private: }); })) { - std::for_each(intersecting_cells.begin(), intersecting_cells.end(), [&](Cell_handle c) { + std::for_each(cr_intersecting_cells.begin(), cr_intersecting_cells.end(), [&](Cell_handle c) { for(int i = 0; i < 4; ++i) { for(int j = i + 1; j < 4; ++j) { auto v1 = c->vertex(i); From 79cb90d5e705b52a2c179a3abf9ebc1b332b7f3a Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 25 Sep 2025 16:28:13 +0200 Subject: [PATCH 19/21] fix warning MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ``` In constructor ‘constexpr std::pair<_T1, _T2>::pair(_U1&&, _U2&&) [with _U1 = CGAL::Triangle_3; _U2 = bool; _T1 = CGAL::Triangle_3; _T2 = bool]’, inlined from ‘constexpr std::pair::type>::__type, typename std::__strip_reference_wrapper::type>::__type> std::make_pair(_T1&&, _T2&&) [with _T1 = CGAL::Triangle_3; _T2 = bool]’ at /usr/include/c++/15/bits/stl_pair.h:1169:72, inlined from ‘std::pair, bool> CGAL::Epic_converter::operator()(const typename IK::Triangle_3&) const [with IK = CGAL::Simple_cartesian >]’ at /mnt/testsuite/include/CGAL/Epic_converter.h:224:28: /usr/include/c++/15/bits/stl_pair.h:464:11: warning: ‘’ may be used uninitialized [-Wmaybe-uninitialized] 464 | : first(std::forward<_U1>(__x)), second(std::forward<_U2>(__y)) | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ``` -> Use the default constructor of the pair, instead of `make_pair`. --- Filtered_kernel/include/CGAL/Epic_converter.h | 76 +++++++++---------- 1 file changed, 38 insertions(+), 38 deletions(-) diff --git a/Filtered_kernel/include/CGAL/Epic_converter.h b/Filtered_kernel/include/CGAL/Epic_converter.h index b8dd072176d..b7d88e77564 100644 --- a/Filtered_kernel/include/CGAL/Epic_converter.h +++ b/Filtered_kernel/include/CGAL/Epic_converter.h @@ -66,7 +66,7 @@ public: if(fit_in_double(n,d)){ return std::make_pair(d,true); } - return std::make_pair(0,false); + return {}; } std::pair operator()(const Bbox_2 b) const @@ -86,7 +86,7 @@ public: if(fit_in_double(p.x(),x) && fit_in_double(p.y(),y)){ return std::make_pair(Point_2(x,y),true); } - return std::make_pair(ORIGIN,false); + return {}; } std::pair operator()(const typename IK::Vector_2& v) const @@ -96,7 +96,7 @@ public: if(fit_in_double(v.x(),x) && fit_in_double(v.y(),y)){ return std::make_pair(Vector_2(x,y),true); } - return std::make_pair(Vector_2(),false); + return {}; } std::pair operator()(const typename IK::Direction_2& d) const @@ -106,7 +106,7 @@ public: if(fit_in_double(d.dx(),x) && fit_in_double(d.dy(),y)){ return std::make_pair(Direction_2(x,y),true); } - return std::make_pair(Direction_2(),false); + return {}; } std::pair operator()(const typename IK::Weighted_point_2& wp) const @@ -116,18 +116,18 @@ public: if(sp.second && w.second){ return std::make_pair(Weighted_point_2(sp.first,w.first),true); } - return std::make_pair(Weighted_point_2(),false); + return {}; } std::pair operator()(const typename IK::Segment_2& s) const { std::pair sp = operator()(s.source()); if(! sp.second){ - return std::make_pair(Segment_2(),false); + return {}; } std::pair tp = operator()(s.target()); if(! tp.second){ - return std::make_pair(Segment_2(),false); + return {}; } return std::make_pair(Segment_2(sp.first,tp.first), true); } @@ -138,18 +138,18 @@ public: if(a.second && b.second && c.second){ return std::make_pair(Line_2(a.first, b.first, c.first),true); } - return std::make_pair(Line_2(), false); + return {}; } std::pair operator()(const typename IK::Ray_2& r) const { std::pair sp = operator()(r.source()); if(! sp.second){ - return std::make_pair(Ray_2(),false); + return {}; } std::pair tp = operator()(r.second_point()); if(! tp.second){ - return std::make_pair(Ray_2(),false); + return {}; } return std::make_pair(Ray_2(sp.first,tp.first), true); } @@ -158,15 +158,15 @@ public: { std::pair v0 = operator()(t.vertex(0)); if(! v0.second){ - return std::make_pair(Triangle_2(),false); + return {}; } std::pair v1 = operator()(t.vertex(1)); if(! v1.second){ - return std::make_pair(Triangle_2(),false); + return {}; } std::pair v2 = operator()(t.vertex(2)); if(! v2.second){ - return std::make_pair(Triangle_2(),false); + return {}; } return std::make_pair(Triangle_2(v0.first,v1.first, v2.first), true); } @@ -178,18 +178,18 @@ public: if(c.second && sr.second){ return std::make_pair(Circle_2(c.first, sr.first, ci.orientation()),true); } - return std::make_pair(Circle_2(), false); + return {}; } std::pair operator()(const typename IK::Iso_rectangle_2& ir) const { std::pair sp = operator()((ir.min)()); if(! sp.second){ - return std::make_pair(Iso_rectangle_2(),false); + return {}; } std::pair tp = operator()((ir.max)()); if(! tp.second){ - return std::make_pair(Iso_rectangle_2(),false); + return {}; } return std::make_pair(Iso_rectangle_2(sp.first,tp.first), true); } @@ -199,11 +199,11 @@ public: { std::pair sp = operator()(li.point()); if(! sp.second){ - return std::make_pair(Line_3(),false); + return {}; } std::pair tp = operator()(li.to_vector()); if(! tp.second){ - return std::make_pair(Line_3(),false); + return {}; } return std::make_pair(Line_3(sp.first,tp.first), true); } @@ -214,22 +214,22 @@ public: if(a.second && b.second && c.second && d.second){ return std::make_pair(Plane_3(a.first, b.first, c.first, d.first),true); } - return std::make_pair(Plane_3(), false); + return {}; } std::pair operator()(const typename IK::Triangle_3& t) const { std::pair v0 = operator()(t.vertex(0)); if(! v0.second){ - return std::make_pair(Triangle_3(),false); + return {}; } std::pair v1 = operator()(t.vertex(1)); if(! v1.second){ - return std::make_pair(Triangle_3(),false); + return {}; } std::pair v2 = operator()(t.vertex(2)); if(! v2.second){ - return std::make_pair(Triangle_3(),false); + return {}; } return std::make_pair(Triangle_3(v0.first,v1.first, v2.first), true); } @@ -238,19 +238,19 @@ public: { std::pair v0 = operator()(t.vertex(0)); if(! v0.second){ - return std::make_pair(Tetrahedron_3(),false); + return {}; } std::pair v1 = operator()(t.vertex(1)); if(! v1.second){ - return std::make_pair(Tetrahedron_3(),false); + return {}; } std::pair v2 = operator()(t.vertex(2)); if(! v2.second){ - return std::make_pair(Tetrahedron_3(),false); + return {}; } std::pair v3 = operator()(t.vertex(3)); if(! v3.second){ - return std::make_pair(Tetrahedron_3(),false); + return {}; } return std::make_pair(Tetrahedron_3(v0.first,v1.first, v2.first, v3.first), true); } @@ -259,11 +259,11 @@ public: { std::pair sp = operator()(r.source()); if(! sp.second){ - return std::make_pair(Ray_3(),false); + return {}; } std::pair tp = operator()(r.second_point()); if(! tp.second){ - return std::make_pair(Ray_3(),false); + return {}; } return std::make_pair(Ray_3(sp.first,tp.first), true); } @@ -275,7 +275,7 @@ public: if(fit_in_double(p.x(),x) && fit_in_double(p.y(),y) && fit_in_double(p.z(),z)){ return std::make_pair(Point_3(x,y,z),true); } - return std::make_pair(ORIGIN,false); + return {}; } std::pair operator()(const typename IK::Vector_3& v) const @@ -285,7 +285,7 @@ public: if(fit_in_double(v.x(),x) && fit_in_double(v.y(),y) && fit_in_double(v.z(),z)){ return std::make_pair(Vector_3(x,y,z),true); } - return std::make_pair(Vector_3(),false); + return {}; } std::pair operator()(const typename IK::Direction_3& d) const @@ -295,18 +295,18 @@ public: if(fit_in_double(d.dx(),x) && fit_in_double(d.dy(),y) && fit_in_double(d.dz(),z)){ return std::make_pair(Direction_3(x,y,z),true); } - return std::make_pair(Direction_3(),false); + return {}; } std::pair operator()(const typename IK::Segment_3& s) const { std::pair sp = operator()(s.source()); if(! sp.second){ - return std::make_pair(Segment_3(),false); + return {}; } std::pair tp = operator()(s.target()); if(! tp.second){ - return std::make_pair(Segment_3(),false); + return {}; } return std::make_pair(Segment_3(sp.first,tp.first), true); } @@ -318,7 +318,7 @@ public: if(sp.second && w.second){ return std::make_pair(Weighted_point_3(sp.first,w.first),true); } - return std::make_pair(Weighted_point_3(),false); + return {}; } std::pair operator()(const typename IK::Sphere_3& s) const @@ -328,7 +328,7 @@ public: if(c.second && sr.second){ return std::make_pair(Sphere_3(c.first, sr.first, s.orientation()),true); } - return std::make_pair(Sphere_3(), false); + return {}; } std::pair operator()(const typename IK::Circle_3& ci) const @@ -338,18 +338,18 @@ public: if(c.second && sr.second){ return std::make_pair(Circle_3(sr.first, c.first),true); } - return std::make_pair(Circle_3(), false); + return {}; } std::pair operator()(const typename IK::Iso_cuboid_3& ic) const { std::pair sp = operator()((ic.min)()); if(! sp.second){ - return std::make_pair(Iso_cuboid_3(),false); + return {}; } std::pair tp = operator()((ic.max)()); if(! tp.second){ - return std::make_pair(Iso_cuboid_3(),false); + return {}; } return std::make_pair(Iso_cuboid_3(sp.first,tp.first), true); } From 757cc8525d34b1fc5ea5f6c18fb651a3d89e9f0f Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 25 Sep 2025 16:36:37 +0200 Subject: [PATCH 20/21] fix a warning MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ``` warning: array subscript 5 is outside array bounds of ‘std::array*, 3> [1]’ [-Warray-bounds=] ``` --- .../CGAL/Conforming_constrained_Delaunay_triangulation_3.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 61c398d6448..241f0e0eb30 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -190,8 +190,10 @@ does_first_triangle_intersect_second_triangle_interior(const typename K::Triangl auto comp = k.compare_xyz_3_object(); auto sort_ptrs = [&comp](const Point_3* p1, const Point_3* p2) { return comp(*p1, *p2) == SMALLER; }; auto intersection_is_a_vertex_or_a_common_edge = [&]() { + CGAL_assume(nb_of_t1_vertices_in_the_line >= 0 && nb_of_t1_vertices_in_the_line <= 3); std::sort(t1_vertices_in_the_line.data(), t1_vertices_in_the_line.data() + nb_of_t1_vertices_in_the_line, sort_ptrs); + CGAL_assume(nb_of_t2_vertices_in_the_line >= 0 && nb_of_t2_vertices_in_the_line <= 3); std::sort(t2_vertices_in_the_line.data(), t2_vertices_in_the_line.data() + nb_of_t2_vertices_in_the_line, sort_ptrs); std::size_t nb_of_common_vertices = 0; From b5a180d9df940ea12e9ac9174ebf568f00216fa6 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 26 Sep 2025 10:22:15 +0200 Subject: [PATCH 21/21] fix the confusion between `handle` and `const_handle` --- .../Conforming_constrained_Delaunay_triangulation_3.h | 2 +- STL_Extension/include/CGAL/iterator.h | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index 241f0e0eb30..21b63c5eb3a 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -2227,7 +2227,7 @@ private: if(it != e) { return *it; } else { - return std::as_const(vertices_of_cavity_union_find).end(); + return vertices_of_cavity_union_find.end(); } }); CGAL_assertion(vertex_below_handle == vertices_of_cavity_union_find.end() || diff --git a/STL_Extension/include/CGAL/iterator.h b/STL_Extension/include/CGAL/iterator.h index 8f8a2228b07..9deb7a73ea0 100644 --- a/STL_Extension/include/CGAL/iterator.h +++ b/STL_Extension/include/CGAL/iterator.h @@ -82,8 +82,14 @@ Iterator_range > make_prevent_deref_range(const I& begin, const } template -auto make_prevent_deref_range(const R& range) +auto make_prevent_deref_range(R&& range) { + static_assert( !std::is_rvalue_reference_v, + "make_prevent_deref_range cannot be used with" + " rvalue references to avoid dangling references"); + // Note: If CGAL were allowed to use C++20, we could use `std::ranges::begin/end`. + // That would allow this to work with rvalue ranges when `std::borrowed_range` is `true`. + // See https://en.cppreference.com/w/cpp/ranges/begin.html#Notes using std::begin; using std::end; return make_range(make_prevent_deref(begin(range)), make_prevent_deref(end(range)));