Reapply changes from PR #5680 (merge conflict resolution)

Outside of the merge commit for clarity.

Distance_3/include/CGAL/Distance_3/Point_3_Line_3.h

@@ -12,7 +12,7 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
 //
 //
-// Author(s)     : Geert-Jan Giezeman, Andreas Fabri
+// Author(s)     : Geert-Jan Giezeman, Andreas Fabri, Mael Rouxel-Labbé
 
 #ifndef CGAL_DISTANCE_3_POINT_3_LINE_3_H
 #define CGAL_DISTANCE_3_POINT_3_LINE_3_H
@@ -26,19 +26,34 @@ namespace CGAL {
 namespace internal {
 
 template <class K>
-typename K::FT
+void
-squared_distance(const typename K::Point_3& pt,
+squared_distance_RT(const typename K::Point_3 &pt,
-                 const typename K::Line_3& line,
+                    const typename K::Line_3 &line,
-                 const K& k)
+                    typename K::RT& num,
+                    typename K::RT& den,
+                    const K& k)
 {
   typedef typename K::Vector_3 Vector_3;
 
   typename K::Construct_vector_3 vector = k.construct_vector_3_object();
 
-  const Vector_3 dir(line.direction().vector());
+  const Vector_3& dir = line.direction().vector();
   const Vector_3 diff = vector(line.point(), pt);
 
-  return internal::squared_distance_to_line(dir, diff, k);
+  return internal::squared_distance_to_line_RT(dir, diff, num, den, k);
+}
+
+template <class K>
+typename K::FT
+squared_distance(const typename K::Point_3& pt,
+                 const typename K::Line_3& line,
+                 const K& k)
+{
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  RT num, den;
+  squared_distance_RT(pt, line, num, den, k);
+  return Rational_traits<FT>().make_rational(num, den);
 }
 
 template <class K>

Distance_3/include/CGAL/Distance_3/Point_3_Ray_3.h

@@ -12,7 +12,7 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
 //
 //
-// Author(s)     : Geert-Jan Giezeman, Andreas Fabri
+// Author(s)     : Geert-Jan Giezeman, Andreas Fabri, Mael Rouxel-Labbé
 
 #ifndef CGAL_DISTANCE_3_POINT_3_RAY_3_H
 #define CGAL_DISTANCE_3_POINT_3_RAY_3_H
@@ -25,18 +25,52 @@
 namespace CGAL {
 namespace internal {
 
+template <class K>
+void
+squared_distance_RT(const typename K::Point_3 &pt,
+                    const typename K::Ray_3 &ray,
+                    typename K::RT& num,
+                    typename K::RT& den,
+                    const K& k)
+{
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::RT RT;
+
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+
+  const Vector_3 dir = ray.direction().vector();
+  const Vector_3 diff = vector(ray.source(), pt);
+
+  if(!is_acute_angle(dir, diff, k))
+  {
+    num = wdot(diff, diff, k);
+    den = wmult((K*)0, RT(1), diff.hw(), diff.hw());
+    return;
+  }
+
+  squared_distance_to_line_RT(dir, diff, num, den, k);
+}
+
 template <class K>
 typename K::FT
 squared_distance(const typename K::Point_3& pt,
                  const typename K::Ray_3& ray,
                  const K& k)
 {
+  // This duplicates code from the _RT functions, but it is a slowdown to do something like:
+  //
+  //   RT num, den;
+  //   squared_distance_RT(pt, ray, num, den, k);
+  //   return Rational_traits<FT>().make_rational(num, den);
+  //
+  // See https://github.com/CGAL/cgal/pull/5680
+
   typedef typename K::Vector_3 Vector_3;
 
   typename K::Construct_vector_3 vector = k.construct_vector_3_object();
 
-  const Vector_3 diff = vector(ray.source(), pt);
   const Vector_3 dir = ray.direction().vector();
+  const Vector_3 diff = vector(ray.source(), pt);
 
   if(!is_acute_angle(dir, diff, k))
     return diff*diff;

Distance_3/include/CGAL/Distance_3/Point_3_Segment_3.h

@@ -27,74 +27,75 @@ namespace CGAL {
 namespace internal {
 
 template <class K>
-typename K::FT
+void
-squared_distance(const typename K::Point_3& pt,
+squared_distance_RT(const typename K::Point_3& pt,
-                 const typename K::Segment_3& seg,
+                    const typename K::Segment_3& seg,
-                 const K& k,
+                    typename K::RT& num,
-                 const Homogeneous_tag&)
+                    typename K::RT& den,
+                    const K& k)
 {
-  typedef typename K::Vector_3 Vector_3;
   typedef typename K::RT RT;
-  typedef typename K::FT FT;
+  typedef typename K::Vector_3 Vector_3;
 
   typename K::Construct_vector_3 vector = k.construct_vector_3_object();
-  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
 
   // assert that the segment is valid (non zero length).
-  Vector_3 diff = vector(seg.source(), pt);
+  const Vector_3 diff_s = vector(seg.source(), pt);
-  Vector_3 segvec = vector(seg.source(), seg.target());
+  const Vector_3 segvec = vector(seg.source(), seg.target());
 
-  RT d = wdot(diff,segvec, k);
+  const RT d = wdot(diff_s, segvec, k);
   if(d <= RT(0))
-    return diff*diff;
+  {
-  RT e = wdot(segvec,segvec, k);
+    // this is squared_distance(pt, seg.source())
-  if((d * segvec.hw()) > (e * diff.hw()))
+    num = wdot(diff_s, diff_s, k);
-    return sq_dist(pt, seg.target());
+    den = wmult((K*)0, RT(1), diff_s.hw(), diff_s.hw());
+    return;
+  }
 
-  Vector_3 wcr = wcross(segvec, diff, k);
+  const RT e = wdot(segvec, segvec, k);
-  return FT(wcr*wcr) / FT(e * diff.hw() * diff.hw());
+  if(wmult((K*)0, d, segvec.hw()) > wmult((K*)0, e, diff_s.hw()))
+  {
+    // this is squared_distance(pt, seg.target())
+    const Vector_3 diff_t = vector(seg.target(), pt);
+    num = wdot(diff_t, diff_t, k);
+    den = wmult((K*)0, RT(1), diff_t.hw(), diff_t.hw());
+    return;
+  }
+
+  // This is an expanded call to squared_distance_to_line_RT() to avoid recomputing 'e'
+  const Vector_3 wcr = wcross(segvec, diff_s, k);
+  num = wdot(wcr, wcr, k);
+  den = wmult((K*)0, e, diff_s.hw(), diff_s.hw());
 }
 
 template <class K>
 typename K::FT
-squared_distance(const typename K::Point_3& pt,
-                 const typename K::Segment_3& seg,
-                 const K& k,
-                 const Cartesian_tag&)
-{
-  typedef typename K::Vector_3 Vector_3;
-  typedef typename K::RT RT;
-  typedef typename K::FT FT;
-
-  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
-  typename K::Compute_squared_distance_3 sqd = k.compute_squared_distance_3_object();
-
-  // assert that the segment is valid (non zero length).
-  Vector_3 diff = vector(seg.source(), pt);
-  Vector_3 segvec = vector(seg.source(), seg.target());
-
-  RT d = wdot(diff,segvec, k);
-  if(d <= RT(0))
-    return diff*diff;
-
-  RT e = wdot(segvec,segvec, k);
-  if(d > e)
-    return sqd(pt, seg.target());
-
-  Vector_3 wcr = wcross(segvec, diff, k);
-  return FT(wcr*wcr)/e;
-}
-
-template <class K>
-inline
-typename K::FT
 squared_distance(const typename K::Point_3& pt,
                  const typename K::Segment_3& seg,
                  const K& k)
 {
-  typedef typename K::Kernel_tag Tag;
+  typedef typename K::RT RT;
-  Tag tag;
+  typedef typename K::FT FT;
-  return squared_distance(pt, seg, k, tag);
+  typedef typename K::Vector_3 Vector_3;
+
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+
+  // assert that the segment is valid (non zero length).
+  const Vector_3 diff = vector(seg.source(), pt);
+  const Vector_3 segvec = vector(seg.source(), seg.target());
+
+  const RT d = wdot(diff, segvec, k);
+  if(d <= RT(0))
+    return (FT(diff*diff));
+
+  const RT e = wdot(segvec, segvec, k);
+  if(wmult((K*)0, d, segvec.hw()) > wmult((K*)0, e, diff.hw()))
+    return squared_distance(pt, seg.target(), k);
+
+  // This is an expanded call to squared_distance_to_line() to avoid recomputing 'e'
+  const Vector_3 wcr = wcross(segvec, diff, k);
+
+  return FT(wcr*wcr) / wmult((K*)0, e, diff.hw(), diff.hw());
 }
 
 template <class K>

Distance_3/include/CGAL/Distance_3/Point_3_Triangle_3.h

@@ -46,6 +46,99 @@ on_left_of_triangle_edge(const typename K::Point_3& pt,
   return (wdot(wcross(edge, normal, k), diff, k) <= RT(0));
 }
 
+template <class K>
+inline void
+squared_distance_to_triangle_RT(const typename K::Point_3& pt,
+                                const typename K::Point_3& t0,
+                                const typename K::Point_3& t1,
+                                const typename K::Point_3& t2,
+                                bool& inside,
+                                typename K::RT& num,
+                                typename K::RT& den,
+                                const K& k)
+{
+  typedef typename K::Vector_3 Vector_3;
+
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+  typename K::Construct_segment_3 segment = k.construct_segment_3_object();
+
+  const Vector_3 e1 = vector(t0, t1);
+  const Vector_3 oe3 = vector(t0, t2);
+  const Vector_3 normal = wcross(e1, oe3, k);
+
+  if(normal == NULL_VECTOR)
+  {
+    // The case normal==NULL_VECTOR covers the case when the triangle
+    // is colinear, or even more degenerate. In that case, we can
+    // simply take also the distance to the three segments.
+    squared_distance_RT(pt, segment(t2, t0), num, den, k);
+
+    typename K::RT num2, den2;
+    squared_distance_RT(pt, segment(t1, t2), num2, den2, k);
+    if(compare_quotients(num2,den2, num,den) == SMALLER)
+    {
+      num = num2;
+      den = den2;
+    }
+
+    // Should not be needed since at most 2 edges cover the full triangle in the degenerate case,
+    // but leaving it for robustness
+    squared_distance_RT(pt, segment(t0, t1), num2, den2, k);
+    if(compare_quotients(num2,den2, num,den) == SMALLER)
+    {
+      num = num2;
+      den = den2;
+    }
+
+    return;
+  }
+
+  const bool b01 = on_left_of_triangle_edge(pt, normal, t0, t1, k);
+  if(!b01)
+  {
+    squared_distance_RT(pt, segment(t0, t1), num, den, k);
+    return;
+  }
+
+  const bool b12 = on_left_of_triangle_edge(pt, normal, t1, t2, k);
+  if(!b12)
+  {
+    squared_distance_RT(pt, segment(t1, t2), num, den, k);
+    return;
+  }
+
+  const bool b20 = on_left_of_triangle_edge(pt, normal, t2, t0, k);
+  if(!b20)
+  {
+    squared_distance_RT(pt, segment(t2, t0), num, den, k);
+    return;
+  }
+
+  // The projection of pt is inside the triangle
+  inside = true;
+  squared_distance_to_plane_RT(normal, vector(t0, pt), num, den, k);
+}
+
+template <class K>
+void
+squared_distance_RT(const typename K::Point_3& pt,
+                    const typename K::Triangle_3& t,
+                    typename K::RT& num,
+                    typename K::RT& den,
+                    const K& k)
+{
+  typename K::Construct_vertex_3 vertex;
+  bool inside = false;
+  squared_distance_to_triangle_RT(pt,
+                                  vertex(t, 0),
+                                  vertex(t, 1),
+                                  vertex(t, 2),
+                                  inside,
+                                  num,
+                                  den,
+                                  k);
+}
+
 template <class K>
 inline typename K::FT
 squared_distance_to_triangle(const typename K::Point_3& pt,
@@ -74,12 +167,11 @@ squared_distance_to_triangle(const typename K::Point_3& pt,
     // Note that in the degenerate case, at most 2 edges cover the full triangle,
     // and only two distances could be used, but leaving 3 for the case of
     // inexact constructions as it might improve the accuracy.
-
     typename K::FT d1 = sq_dist(pt, segment(t2, t0));
     typename K::FT d2 = sq_dist(pt, segment(t1, t2));
     typename K::FT d3 = sq_dist(pt, segment(t0, t1));
 
-    return (std::min)( (std::min)(d1, d2), d3);
+    return (std::min)((std::min)(d1, d2), d3);
   }
 
   const bool b01 = on_left_of_triangle_edge(pt, normal, t0, t1, k);

Distance_3/include/CGAL/Distance_3/internal/squared_distance_utils_3.h

@@ -167,8 +167,8 @@ is_obtuse_angle(const typename K::Point_3 &p,
 }
 template <class K>
 void
-squared_distance_to_plane_RT(const typename K::Vector_3 & normal,
+squared_distance_to_plane_RT(const typename K::Vector_3& normal,
-                             const typename K::Vector_3 & diff,
+                             const typename K::Vector_3& diff,
                              typename K::RT& num,
                              typename K::RT& den,
                              const K& k)
@@ -183,8 +183,8 @@ squared_distance_to_plane_RT(const typename K::Vector_3 & normal,
 
 template <class K>
 typename K::FT
-squared_distance_to_plane(const typename K::Vector_3 & normal,
+squared_distance_to_plane(const typename K::Vector_3& normal,
-                          const typename K::Vector_3 & diff,
+                          const typename K::Vector_3& diff,
                           const K& k)
 {
   typedef typename K::RT RT;
@@ -196,8 +196,8 @@ squared_distance_to_plane(const typename K::Vector_3 & normal,
 
 template <class K>
 void
-squared_distance_to_line_RT(const typename K::Vector_3 & dir,
+squared_distance_to_line_RT(const typename K::Vector_3& dir,
-                            const typename K::Vector_3 & diff,
+                            const typename K::Vector_3& diff,
                             typename K::RT& num,
                             typename K::RT& den,
                             const K& k)
@@ -207,10 +207,11 @@ squared_distance_to_line_RT(const typename K::Vector_3 & dir,
   num = wdot(wcr, wcr, k);
   den = wmult((K*)0, wdot(dir, dir, k), diff.hw(), diff.hw());
 }
+
 template <class K>
 typename K::FT
-squared_distance_to_line(const typename K::Vector_3 & dir,
+squared_distance_to_line(const typename K::Vector_3& dir,
-                         const typename K::Vector_3 & diff,
+                         const typename K::Vector_3& diff,
                          const K& k)
 {
   typedef typename K::RT RT;