From 8f551aa20a8f2772f5dcf9484c981e1282bb8b4e Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 6 Nov 2025 12:09:13 +0100 Subject: [PATCH] optional exact constructions of Steiner points --- .../Conforming_Delaunay_triangulation_3.h | 182 +++++++++++++++--- ...ing_constrained_Delaunay_triangulation_3.h | 45 ++++- STL_Extension/include/CGAL/functional.h | 11 ++ 3 files changed, 209 insertions(+), 29 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 5939f7d035c..ac5888d420e 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h @@ -16,10 +16,20 @@ #include +#include +#include +#include +#include #include +#include +#include +#include +#include #include +#include #include #include +#include #include #include @@ -30,8 +40,21 @@ #include #include -#include +#include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #ifndef DOXYGEN_RUNNING @@ -64,6 +87,7 @@ struct Debug_options { debug_polygon_insertion, display_statistics, use_epeck_for_normals, + use_epeck_for_Steiner_points, nb_of_flags }; @@ -122,6 +146,9 @@ struct Debug_options { bool use_epeck_for_normals() const { return flags[static_cast(Flags::use_epeck_for_normals)]; } void use_epeck_for_normals(bool b) { flags.set(static_cast(Flags::use_epeck_for_normals), b); } + bool use_epeck_for_Steiner_points() const { return flags[static_cast(Flags::use_epeck_for_Steiner_points)]; } + void use_epeck_for_Steiner_points(bool b) { flags.set(static_cast(Flags::use_epeck_for_Steiner_points), b); } + double segment_vertex_epsilon() const { return segment_vertex_epsilon_; } void set_segment_vertex_epsilon(double eps) { segment_vertex_epsilon_ = eps; } @@ -256,7 +283,7 @@ public: using Line = typename T_3::Geom_traits::Line_3; using Locate_type = typename T_3::Locate_type; - inline static With_offset_tag with_offset{ {-1} }; + inline static With_offset_tag with_offset{ -1 }; inline static With_point_tag with_point{ {-1} }; inline static With_point_and_info_tag with_point_and_info{ {-1} }; @@ -971,6 +998,95 @@ protected: return vector_of_encroaching_vertices; } + // Helper to compute a projected point with optional exact kernel and custom threshold check + // lambda_computer receives (kernel, converter) and returns the lambda parameter for projection + // use_midpoint_check receives (lambda, std::optional) and returns true if midpoint should be used instead + // - First call with std::nullopt allows checking lambda before computing projection + // - Second call with actual projected point allows distance-based checks + template + Point compute_projected_point_with_threshold(const Point& start_pt, const Point& end_pt, + const Point& midpoint_start, const Point& midpoint_end, + LambdaComputer&& lambda_computer, + MidpointCheck&& use_midpoint_check) const + { + auto& gt = tr().geom_traits(); + auto vector_functor = gt.construct_vector_3_object(); + auto midpoint_functor = gt.construct_midpoint_3_object(); + auto scaled_vector_functor = gt.construct_scaled_vector_3_object(); + auto translate_functor = gt.construct_translated_point_3_object(); + + if constexpr (!Algebraic_structure_traits::Is_exact::value) { + if(debug().use_epeck_for_Steiner_points()) { + return std::invoke([&]() { + using Epeck_ft = internal::Exact_field_selector::Type; + using Exact_kernel = Simple_cartesian; + Exact_kernel exact_kernel; + Cartesian_converter back_from_exact; + Cartesian_converter to_exact; + + auto&& exact_vector = exact_kernel.construct_vector_3_object(); + auto&& exact_midpoint = exact_kernel.construct_midpoint_3_object(); + auto&& exact_scaled_vector = exact_kernel.construct_scaled_vector_3_object(); + auto&& exact_translate = exact_kernel.construct_translated_point_3_object(); + + const auto lambda = lambda_computer(exact_kernel, to_exact); + + // Check threshold before computing projection + if(use_midpoint_check(lambda, std::nullopt)) { + // Only convert midpoint points when needed + const auto midpoint_start_exact = to_exact(midpoint_start); + const auto midpoint_end_exact = to_exact(midpoint_end); + return back_from_exact(exact_midpoint(midpoint_start_exact, midpoint_end_exact)); + } + + // Only convert projection points when needed + const auto start_pt_exact = to_exact(start_pt); + const auto end_pt_exact = to_exact(end_pt); + const auto vector_exact = exact_vector(start_pt_exact, end_pt_exact); + const auto projected_exact = exact_translate(start_pt_exact, exact_scaled_vector(vector_exact, lambda)); + const auto projected_approx = back_from_exact(projected_exact); + + // Second threshold check with actual projected point if needed + if(use_midpoint_check(lambda, projected_approx)) { + const auto midpoint_start_exact = to_exact(midpoint_start); + const auto midpoint_end_exact = to_exact(midpoint_end); + return back_from_exact(exact_midpoint(midpoint_start_exact, midpoint_end_exact)); + } + + return projected_approx; + }); + } + } + return std::invoke([&]() { + const auto lambda = lambda_computer(gt, CGAL::cpp20::identity{}); + + // Check threshold before computing projection + if(use_midpoint_check(lambda, std::nullopt)) { + return midpoint_functor(midpoint_start, midpoint_end); + } + + const auto vector_ab = vector_functor(start_pt, end_pt); + const auto projected_pt = translate_functor(start_pt, scaled_vector_functor(vector_ab, lambda)); + + // Second threshold check with actual projected point if needed + return use_midpoint_check(lambda, projected_pt) ? midpoint_functor(midpoint_start, midpoint_end) : projected_pt; + }); + } + + // Convenience wrapper for simple lambda-based threshold (lambda < 0.2 || lambda > 0.8) + template + Point compute_projected_point(const Point& start_pt, const Point& end_pt, + LambdaComputer&& lambda_computer) const + { + return compute_projected_point_with_threshold( + start_pt, end_pt, start_pt, end_pt, + std::forward(lambda_computer), + [](auto lambda, const std::optional&) { + // Only need lambda for this threshold check + return lambda < 0.2 || lambda > 0.8; + }); + } + Construct_Steiner_point_return_type construct_Steiner_point(Constrained_polyline_id constrained_polyline_id, Subconstraint subconstraint) { @@ -978,10 +1094,7 @@ protected: auto compare_angle_functor = gt.compare_angle_3_object(); auto vector_functor = gt.construct_vector_3_object(); auto midpoint_functor = gt.construct_midpoint_3_object(); - auto scaled_vector_functor = gt.construct_scaled_vector_3_object(); auto sq_length_functor = gt.compute_squared_length_3_object(); - auto sc_product_functor = gt.compute_scalar_product_3_object(); - auto translate_functor = gt.construct_translated_point_3_object(); const Vertex_handle va = subconstraint.first; const Vertex_handle vb = subconstraint.second; @@ -1033,24 +1146,30 @@ protected: #endif // CGAL_CDT_3_DEBUG_CONFORMING const auto vector_orig_ab = vector_functor(orig_pa, orig_pb); const auto length_ab = CGAL::approximate_sqrt(sq_length_functor(vector_ab)); - auto return_orig_result_point = - [&](auto lambda, Point orig_pa, Point orig_pb) - -> Construct_Steiner_point_return_type - { - const auto vector_orig_ab = vector_functor(orig_pa, orig_pb); - const auto inter_point = translate_functor(orig_pa, scaled_vector_functor(vector_orig_ab, lambda)); - const auto dist_a_result = CGAL::approximate_sqrt(sq_length_functor(vector_functor(pa, inter_point))); - const auto ratio = dist_a_result / length_ab; - const auto result_point = (ratio < 0.2 || ratio > 0.8) - ? midpoint_functor(pa, pb) - : inter_point; + auto return_orig_result_point = [&](auto lambda_val, Point orig_pa_param, + Point orig_pb_param) -> Construct_Steiner_point_return_type { + // Compute projected point with distance-based ratio threshold check + const auto result_point = compute_projected_point_with_threshold( + orig_pa_param, orig_pb_param, // Projection segment + pa, pb, // Midpoint segment + [lambda_val](auto&&, auto&&) { return lambda_val; }, // Lambda computer + [&](auto, const std::optional& projected_pt) { // Threshold check based on distance ratio + // If projected_pt is not computed yet (std::nullopt), can't check, so return false + if(!projected_pt) return false; + + const auto dist_a_result = CGAL::approximate_sqrt(sq_length_functor(vector_functor(pa, *projected_pt))); + const auto ratio = dist_a_result / length_ab; +#if CGAL_CDT_3_DEBUG_CONFORMING + std::cerr << " ref ratio = " << ratio << '\n'; +#endif + return ratio < 0.2 || ratio > 0.8; + }); #if CGAL_CDT_3_DEBUG_CONFORMING - std::cerr << " ref ratio = " << ratio << '\n'; - std::cerr << " -> Steiner point: " << result_point << '\n'; + std::cerr << " -> Steiner point: " << result_point << '\n'; #endif // CGAL_CDT_3_DEBUG_CONFORMING - return {result_point, reference_vertex->cell(), reference_vertex}; - }; + return {result_point, reference_vertex->cell(), reference_vertex}; + }; const auto length_orig_ab = CGAL::approximate_sqrt(sq_length_functor(vector_orig_ab)); if(ref_va == orig_va || ref_vb == orig_va) { @@ -1075,13 +1194,26 @@ protected: } } // compute the projection of the reference point - const auto vector_a_ref = vector_functor(pa, reference_point); - const auto lambda = sc_product_functor(vector_a_ref, vector_ab) / sq_length_functor(vector_ab); - const auto result_point = (lambda < 0.2 || lambda > 0.8) - ? midpoint_functor(pa, pb) - : translate_functor(pa, scaled_vector_functor(vector_ab, lambda)); + + const auto result_point = compute_projected_point( + pa, pb, + [&](auto&& kernel, auto&& converter) { + auto&& vec_func = kernel.construct_vector_3_object(); + auto&& sc_prod_func = kernel.compute_scalar_product_3_object(); + auto&& sq_len_func = kernel.compute_squared_length_3_object(); + + const auto pa_converted = converter(pa); + const auto pb_converted = converter(pb); + const auto ref_pt_converted = converter(reference_point); + + const auto vector_ab = vec_func(pa_converted, pb_converted); + const auto vector_a_ref = vec_func(pa_converted, ref_pt_converted); + return sc_prod_func(vector_a_ref, vector_ab) / sq_len_func(vector_ab); + }); #if CGAL_CDT_3_DEBUG_CONFORMING + const auto vector_a_ref = vector_functor(pa, reference_point); + const auto lambda = sc_product_functor(vector_a_ref, vector_ab) / sq_length_functor(vector_ab); std::cerr << " lambda = " << lambda << '\n'; std::cerr << " -> Steiner point: " << result_point << '\n'; #endif // CGAL_CDT_3_DEBUG_CONFORMING 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 56547d232ee..4210af8a1f3 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 @@ -1825,7 +1825,7 @@ public: auto& borders = face_index.has_value() ? this->non_const_face_borders(face_index.value()) - : this->face_data.emplace_back(CDT_2_traits{compute_accumulated_normal(first_it, size)}).borders; + : this->face_data.emplace_back(CDT_2_traits{construct_accumulated_normal(first_it, size)}).borders; auto& border = borders.emplace_back(); border.reserve(std::size(vertex_handles)); const auto polygon_constraint_id = face_index.value_or(this->face_data.size() - 1); @@ -2492,8 +2492,41 @@ private: Next_region(const std::string& what, CDT_2_face_handle fh) : std::logic_error(what), fh_2d(fh) {} }; + template + Point_3 construct_with_optional_exact_kernel(GetConstructor&& get_ctor, Args&&... args) const + { + if(this->debug().use_epeck_for_Steiner_points()) { + // Use exact Epeck kernel for robust computation + using Epeck_ft = internal::Exact_field_selector::Type; + using Exact_kernel = Simple_cartesian; + Exact_kernel exact_kernel; + Cartesian_converter back_from_exact; + Cartesian_converter to_exact; + + auto&& construct = get_ctor(exact_kernel); + return back_from_exact(construct(to_exact(args)...)); + } else { + auto&& construct = get_ctor(tr().geom_traits()); + return construct(args...); + } + } + + Point_3 construct_midpoint(const Point_3& a, const Point_3& b) const + { + return construct_with_optional_exact_kernel( + [](auto&& kernel) { return kernel.construct_midpoint_3_object(); }, + a, b); + } + + Point_3 construct_centroid(const Point_3& p1, const Point_3& p2, const Point_3& p3) const + { + return construct_with_optional_exact_kernel( + [](auto&& kernel) { return kernel.construct_centroid_3_object(); }, + p1, p2, p3); + } + template - Vector_3 compute_accumulated_normal(VertexIterator first_it, std::size_t size) const + Vector_3 construct_accumulated_normal(VertexIterator first_it, std::size_t size) const { const auto before_last_it = std::next(first_it, size - 2); const auto last_it = std::next(before_last_it); @@ -3627,7 +3660,10 @@ private: CDT_2_face_handle fh_2d) { const auto& cdt_2 = non_const_cdt_2; - auto steiner_pt = CGAL::centroid(cdt_2.triangle(fh_2d)); + + auto steiner_pt = construct_centroid(cdt_2.point(fh_2d->vertex(0)), + cdt_2.point(fh_2d->vertex(1)), + cdt_2.point(fh_2d->vertex(2))); if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug().verbose_special_cases()) { std::cerr << cdt_3_format("Trying to insert Steiner (centroid) point {} in non-coplanar face {}.\n", IO::oformat(steiner_pt), IO::oformat(cdt_2.triangle(fh_2d))); @@ -3728,7 +3764,8 @@ private: void insert_mid_point_in_constrained_edge(Vertex_handle va_3d, Vertex_handle vb_3d) { const auto a = this->point(va_3d); const auto b = this->point(vb_3d); - const auto mid = CGAL::midpoint(a, b); + + const auto mid = construct_midpoint(a, b); if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug().Steiner_points()) { std::cerr << cdt_3_format("Inserting Steiner (midpoint) point {} of constrained edge ({:.6} , {:.6})\n", IO::oformat(mid), IO::oformat(va_3d, with_point_and_info), diff --git a/STL_Extension/include/CGAL/functional.h b/STL_Extension/include/CGAL/functional.h index c440c5e71f3..6510235e025 100644 --- a/STL_Extension/include/CGAL/functional.h +++ b/STL_Extension/include/CGAL/functional.h @@ -40,6 +40,17 @@ namespace cpp98 { } // namespace cpp98 +namespace cpp20 { + +/// Replacement for `std::identity` that is added in C++20. +struct identity +{ + template + constexpr T&& operator()(T&& t) const noexcept { return static_cast(t); } +}; + +} // namespace cpp20 + } // namespace CGAL #endif // CGAL_FUNCTIONAL_H