From fccbeaed5984fb6b3d1eacb65384976a77e7b134 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 6 Mar 2024 16:26:45 +0100 Subject: [PATCH] change the distance check between constraint segment and vertex Now it is done up-front, and not during the execution of the process. --- .../Conforming_Delaunay_triangulation_3.h | 83 ++++++++++--------- .../test/Triangulation_3/cdt_3_from_off.cpp | 70 ++++++++++++++-- 2 files changed, 108 insertions(+), 45 deletions(-) diff --git a/Triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h b/Triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h index 2add2c9db2a..53fd471d450 100644 --- a/Triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h @@ -248,20 +248,9 @@ protected: Visitor&) { if(va != vb) { if(segment_vertex_epsilon != 0.) { - if(!max_bbox_edge_length) { - update_max_bbox_edge_length(); - } auto [min_dist, min_vertex] = min_distance_and_vertex_between_constraint_and_encroaching_vertex(va, vb); - if(min_dist < segment_vertex_epsilon * *max_bbox_edge_length) { - std::stringstream ss; - ss.precision(std::cerr.precision()); - ss << "A constrained segment is too close to a vertex.\n"; - ss << " -> vertex " << display_vert(min_vertex) << '\n'; - ss << " -> constrained segment " << display_vert(va) << " - " << display_vert(vb) << '\n'; - ss << " -> distance = " << min_dist << '\n'; - ss << " -> max_bbox_edge_length = " << *max_bbox_edge_length << '\n'; - throw std::runtime_error(ss.str()); - } + check_segment_vertex_distance_or_throw(va, vb, min_vertex, CGAL::to_double(min_dist), + Check_distance::NON_SQUARED_DISTANCE); } const Constraint_id c_id = constraint_hierarchy.insert_constraint(va, vb); pair_of_vertices_to_cid.emplace(make_sorted_pair(va, vb), c_id); @@ -439,6 +428,32 @@ public: }); } + enum class Check_distance { SQUARED_DISTANCE, NON_SQUARED_DISTANCE }; + + void check_segment_vertex_distance_or_throw(Vertex_handle va, + Vertex_handle vb, + Vertex_handle min_vertex, + double min_dist, + Check_distance option) + { + if(!max_bbox_edge_length) { + update_max_bbox_edge_length(); + } + if((option == Check_distance::NON_SQUARED_DISTANCE && min_dist < segment_vertex_epsilon * *max_bbox_edge_length) || + (option == Check_distance::SQUARED_DISTANCE && + min_dist < CGAL::square(segment_vertex_epsilon * *max_bbox_edge_length))) + { + std::stringstream ss; + ss.precision(std::cerr.precision()); + ss << "A constrained segment is too close to a vertex.\n"; + ss << " -> vertex " << display_vert(min_vertex) << '\n'; + ss << " -> constrained segment " << display_vert(va) << " - " << display_vert(vb) << '\n'; + ss << " -> distance = " << min_dist << '\n'; + ss << " -> max_bbox_edge_length = " << *max_bbox_edge_length << '\n'; + throw std::runtime_error(ss.str()); + } + } + auto ancestors_of_Steiner_vertex_on_edge(Vertex_handle v) const { std::pair result; CGAL_precondition(v->is_Steiner_vertex_on_edge()); @@ -471,6 +486,21 @@ public: return result; } + auto min_distance_and_vertex_between_constraint_and_encroaching_vertex(Vertex_handle va, Vertex_handle vb) const { + struct Result { + typename T_3::Geom_traits::FT min_dist = std::numeric_limits::max(); + Vertex_handle v = {}; + } result; + const auto vector_of_encroaching_vertices = encroaching_vertices(va, vb); + for(auto v: vector_of_encroaching_vertices) { + const auto dist = CGAL::approximate_sqrt(squared_distance(tr.point(v), Line{tr.point(va), tr.point(vb)})); + if(dist < result.min_dist) { + result = Result{dist, v}; + } + } + return result; + } + bool write_missing_segments_file(std::ostream &out) { bool any_missing_segment = false; std::for_each( @@ -753,21 +783,6 @@ protected: return vector_of_encroaching_vertices; } - auto min_distance_and_vertex_between_constraint_and_encroaching_vertex(Vertex_handle va, Vertex_handle vb) const { - struct Result { - typename T_3::Geom_traits::FT min_dist = std::numeric_limits::max(); - Vertex_handle v = {}; - } result; - const auto vector_of_encroaching_vertices = encroaching_vertices(va, vb); - for(auto v: vector_of_encroaching_vertices) { - const auto dist = CGAL::approximate_sqrt(squared_distance(tr.point(v), Line{tr.point(va), tr.point(vb)})); - if(dist < result.min_dist) { - result = Result{dist, v}; - } - } - return result; - } - Construct_Steiner_point_return_type construct_Steiner_point(Constraint_id constraint_id, Subconstraint subconstraint) { @@ -867,16 +882,8 @@ protected: update_max_bbox_edge_length(); } auto sq_dist = squared_distance(reference_point, Line{orig_pa, orig_pb}); - if(sq_dist < CGAL::square(segment_vertex_epsilon * *max_bbox_edge_length)) { - std::stringstream ss; - ss.precision(std::cerr.precision()); - ss << "A constrained segment is too close to a vertex.\n"; - ss << " -> vertex " << display_vert(reference_vertex) << '\n'; - ss << " -> constrained segment " << display_vert(orig_va) << " - " << display_vert(orig_vb) << '\n'; - ss << " -> distance = " << CGAL::approximate_sqrt(sq_dist) << '\n'; - ss << " -> max_bbox_edge_length = " << *max_bbox_edge_length << '\n'; - throw std::runtime_error(ss.str()); - } + check_segment_vertex_distance_or_throw(orig_va, orig_vb, reference_vertex, CGAL::to_double(sq_dist), + Check_distance::SQUARED_DISTANCE); } } // compute the projection of the reference point diff --git a/Triangulation_3/test/Triangulation_3/cdt_3_from_off.cpp b/Triangulation_3/test/Triangulation_3/cdt_3_from_off.cpp index d160c36d5e2..db73586c264 100644 --- a/Triangulation_3/test/Triangulation_3/cdt_3_from_off.cpp +++ b/Triangulation_3/test/Triangulation_3/cdt_3_from_off.cpp @@ -71,8 +71,8 @@ struct CDT_options bool use_new_cavity_algorithm = true; bool debug_validity = false; double ratio = 0.1; - double vertex_vertex_epsilon = 1e-6; - double segment_vertex_epsilon = 1e-8; + double vertex_vertex_epsilon = 1e-14; + double segment_vertex_epsilon = 1e-14; std::string failure_assertion_expression{}; std::string input_filename = CGAL::data_file_path("meshes/mpi.off"); std::string output_filename{"dump.off"}; @@ -357,10 +357,10 @@ int go(Mesh mesh, CDT_options options) { auto pmap = get(CGAL::vertex_point, mesh); - auto [patch_id_map, ok] = mesh.add_property_map("f:patch_id", -2); - assert(ok); CGAL_USE(ok); - auto [v_selected_map, ok2] = mesh.add_property_map("v:selected", false); - assert(ok2); CGAL_USE(ok2); + auto [patch_id_map, patch_id_map_ok] = mesh.add_property_map("f:patch_id", -2); + assert(patch_id_map_ok); CGAL_USE(patch_id_map_ok); + auto [v_selected_map, v_selected_map_ok] = mesh.add_property_map("v:selected", false); + assert(v_selected_map_ok); CGAL_USE(v_selected_map_ok); int nb_patches = 0; std::vector>> patch_edges; if(options.merge_facets) { @@ -490,10 +490,14 @@ int go(Mesh mesh, CDT_options options) { } }; + auto [tr_vertex_pmap, tr_vertex_pmap_ok] = mesh.add_property_map("tr_vertex"); + assert(tr_vertex_pmap_ok); CGAL_USE(tr_vertex_pmap_ok); + auto start_time = std::chrono::high_resolution_clock::now(); for(auto v: vertices(mesh)) { if(options.merge_facets && false == get(v_selected_map, v)) continue; - cdt.insert(get(pmap, v)); + auto vh = cdt.insert(get(pmap, v)); + put(tr_vertex_pmap, v, vh); } if(!options.quiet) { std::cout << "[timings] inserted vertices in " << std::chrono::duration_cast( @@ -517,6 +521,7 @@ int go(Mesh mesh, CDT_options options) { cdt.insert(Point(bbox.xmax() + d_x, bbox.ymin() - d_y, bbox.zmax() + d_z)); cdt.insert(Point(bbox.xmax() + d_x, bbox.ymax() + d_y, bbox.zmax() + d_z)); } + start_time = std::chrono::high_resolution_clock::now(); { auto [min_sq_distance, min_edge] = std::ranges::min( cdt.finite_edges() | std::views::transform([&](auto edge) { return std::make_pair(cdt.segment(edge).squared_length(), edge); })); @@ -533,6 +538,57 @@ int go(Mesh mesh, CDT_options options) { return exit_code; } } + { + double min_distance = std::numeric_limits::max(); + CDT::Vertex_handle min_va, min_vb, min_vertex; + if(options.merge_facets) { + for(int i = 0; i < nb_patches; ++i) { + const auto& edges = patch_edges[i]; + for(auto [vda, vdb]: edges) { + auto va = get(tr_vertex_pmap, vda); + auto vb = get(tr_vertex_pmap, vdb); + auto [min_dist, min_v] = cdt.min_distance_and_vertex_between_constraint_and_encroaching_vertex(va, vb); + if(min_dist < min_distance) { + min_distance = CGAL::to_double(min_dist); + min_va = va; + min_vb = vb; + min_vertex = min_v; + } + } + } + } else { + for(auto face_descriptor : faces(mesh)) { + auto he = halfedge(face_descriptor, mesh); + const auto end = he; + do { + auto va = get(tr_vertex_pmap, source(he, mesh)); + auto vb = get(tr_vertex_pmap, target(he, mesh)); + auto [min_dist, min_v] = cdt.min_distance_and_vertex_between_constraint_and_encroaching_vertex(va, vb); + if(min_dist < min_distance) { + min_distance = CGAL::to_double(min_dist); + min_va = va; + min_vb = vb; + min_vertex = min_v; + } + he = next(he, mesh); + } while((he = next(he, mesh)) != end); + } + } + if(!options.quiet) { + std::cout << "Min distance between constraint segment and vertex: " << min_distance << '\n' + << " between segment : " + << CGAL::IO::oformat(min_va, CDT::Conforming_Dt::with_point) << " " + << CGAL::IO::oformat(min_vb, CDT::Conforming_Dt::with_point) << '\n' + << " and vertex: : " + << CGAL::IO::oformat(min_vertex, CDT::Conforming_Dt::with_point) << "\n\n"; + } + cdt.check_segment_vertex_distance_or_throw(min_va, min_vb, min_vertex, min_distance, + CDT::Check_distance::NON_SQUARED_DISTANCE); + } + if(!options.quiet) { + std::cout << "[timings] compute distances on " << std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - start_time).count() << " ms\n"; + } int poly_id = 0; CDT_3_try { start_time = std::chrono::high_resolution_clock::now();