optional exact constructions of Steiner points

Constrained_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h

@@ -16,10 +16,20 @@
 
 #include <CGAL/Constrained_triangulation_3/internal/config.h>
 
+#include <CGAL/Algebraic_structure_traits.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Cartesian_converter.h>
+#include <CGAL/Compact_container.h>
 #include <CGAL/Conforming_constrained_Delaunay_triangulation_vertex_data_3.h>
+#include <CGAL/enum.h>
+#include <CGAL/functional.h>
+#include <CGAL/kernel_assertions.h>
+#include <CGAL/Number_types/internal/Exact_type_selector.h>
 #include <CGAL/Real_timer.h>
+#include <CGAL/SMDS_3/io_signature.h>
 #include <CGAL/Triangulation_2/internal/Polyline_constraint_hierarchy_2.h>
 #include <CGAL/Triangulation_segment_traverser_3.h>
+#include <CGAL/unordered_flat_map.h>
 #include <CGAL/unordered_flat_set.h>
 
 #include <CGAL/Mesh_3/io_signature.h>
@@ -30,8 +40,21 @@
 #include <boost/container/small_vector.hpp>
 #include <boost/iterator/function_output_iterator.hpp>
 
-#include <fstream>
+#include <algorithm>
+#include <array>
 #include <bitset>
+#include <fstream>
+#include <functional>
+#include <ios>
+#include <iostream>
+#include <limits>
+#include <optional>
+#include <ostream>
+#include <sstream>
+#include <stack>
+#include <string>
+#include <utility>
+#include <vector>
 
 #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<int>(Flags::use_epeck_for_normals)]; }
   void use_epeck_for_normals(bool b) { flags.set(static_cast<int>(Flags::use_epeck_for_normals), b); }
 
+  bool use_epeck_for_Steiner_points() const { return flags[static_cast<int>(Flags::use_epeck_for_Steiner_points)]; }
+  void use_epeck_for_Steiner_points(bool b) { flags.set(static_cast<int>(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<Point>) 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<typename LambdaComputer, typename MidpointCheck>
+  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<typename T_3::Geom_traits::FT>::Is_exact::value) {
+      if(debug().use_epeck_for_Steiner_points()) {
+        return std::invoke([&]() {
+          using Epeck_ft = internal::Exact_field_selector<double>::Type;
+          using Exact_kernel = Simple_cartesian<Epeck_ft>;
+          Exact_kernel exact_kernel;
+          Cartesian_converter<Exact_kernel, Geom_traits> back_from_exact;
+          Cartesian_converter<Geom_traits, Exact_kernel> 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<typename LambdaComputer>
+  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<LambdaComputer>(lambda_computer),
+        [](auto lambda, const std::optional<Point>&) {
+          // 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 return_orig_result_point = [&](auto lambda_val, Point orig_pa_param,
-          [&](auto lambda, Point orig_pa, Point orig_pb)
+                                          Point orig_pb_param) -> Construct_Steiner_point_return_type {
-              -> Construct_Steiner_point_return_type
+        // Compute projected point with distance-based ratio threshold check
-          {
+        const auto result_point = compute_projected_point_with_threshold(
-            const auto vector_orig_ab = vector_functor(orig_pa, orig_pb);
+            orig_pa_param, orig_pb_param,  // Projection segment
-            const auto inter_point = translate_functor(orig_pa, scaled_vector_functor(vector_orig_ab, lambda));
+            pa, pb,                         // Midpoint segment
-            const auto dist_a_result = CGAL::approximate_sqrt(sq_length_functor(vector_functor(pa, inter_point)));
+            [lambda_val](auto&&, auto&&) { return lambda_val; },  // Lambda computer
-            const auto ratio = dist_a_result / length_ab;
+            [&](auto, const std::optional<Point>& projected_pt) {  // Threshold check based on distance ratio
-            const auto result_point = (ratio < 0.2 || ratio > 0.8)
+              // If projected_pt is not computed yet (std::nullopt), can't check, so return false
-                                          ? midpoint_functor(pa, pb)
+              if(!projected_pt) return false;
-                                          : inter_point;
+
+              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 = compute_projected_point(
-    const auto result_point = (lambda < 0.2 || lambda > 0.8)
+        pa, pb,
-                                  ? midpoint_functor(pa, pb)
+        [&](auto&& kernel, auto&& converter) {
-                                  : translate_functor(pa, scaled_vector_functor(vector_ab, lambda));
+          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

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<typename GetConstructor, typename... Args>
+  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<double>::Type;
+      using Exact_kernel = Simple_cartesian<Epeck_ft>;
+      Exact_kernel exact_kernel;
+      Cartesian_converter<Exact_kernel, Geom_traits> back_from_exact;
+      Cartesian_converter<Geom_traits, Exact_kernel> 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 <typename VertexIterator>
-  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),

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 <typename T>
+  constexpr T&& operator()(T&& t) const noexcept { return static_cast<T&&>(t); }
+};
+
+} // namespace cpp20
+
 } // namespace CGAL
 
 #endif // CGAL_FUNCTIONAL_H