From a4a00d969d5e1a8f7ab896b93a76aab157071a40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Wed, 23 Jun 2021 22:30:23 +0200 Subject: [PATCH] Reapply changes from PR #5680 (merge conflict resolution) Outside of the merge commit for clarity. --- .../include/CGAL/Distance_3/Point_3_Line_3.h | 29 +++-- .../include/CGAL/Distance_3/Point_3_Ray_3.h | 38 ++++++- .../CGAL/Distance_3/Point_3_Segment_3.h | 103 +++++++++--------- .../CGAL/Distance_3/Point_3_Triangle_3.h | 96 +++++++++++++++- .../internal/squared_distance_utils_3.h | 17 +-- 5 files changed, 213 insertions(+), 70 deletions(-) diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Line_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Line_3.h index 24f7593d51d..0d1527c48a7 100644 --- a/Distance_3/include/CGAL/Distance_3/Point_3_Line_3.h +++ b/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 -typename K::FT -squared_distance(const typename K::Point_3& pt, - const typename K::Line_3& line, - const K& k) +void +squared_distance_RT(const typename K::Point_3 &pt, + const typename K::Line_3 &line, + 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 +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().make_rational(num, den); } template diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Ray_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Ray_3.h index 91c8965d26f..7eff899b0eb 100644 --- a/Distance_3/include/CGAL/Distance_3/Point_3_Ray_3.h +++ b/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 +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 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().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; diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Segment_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Segment_3.h index 34ee169c8e1..554c0a58d2d 100644 --- a/Distance_3/include/CGAL/Distance_3/Point_3_Segment_3.h +++ b/Distance_3/include/CGAL/Distance_3/Point_3_Segment_3.h @@ -27,74 +27,75 @@ namespace CGAL { namespace internal { template -typename K::FT -squared_distance(const typename K::Point_3& pt, - const typename K::Segment_3& seg, - const K& k, - const Homogeneous_tag&) +void +squared_distance_RT(const typename K::Point_3& pt, + const typename K::Segment_3& seg, + typename K::RT& num, + 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); - Vector_3 segvec = vector(seg.source(), seg.target()); + const Vector_3 diff_s = vector(seg.source(), pt); + 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); - if((d * segvec.hw()) > (e * diff.hw())) - return sq_dist(pt, seg.target()); + { + // this is squared_distance(pt, seg.source()) + num = wdot(diff_s, diff_s, k); + den = wmult((K*)0, RT(1), diff_s.hw(), diff_s.hw()); + return; + } - Vector_3 wcr = wcross(segvec, diff, k); - return FT(wcr*wcr) / FT(e * diff.hw() * diff.hw()); + const RT e = wdot(segvec, segvec, k); + 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 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 -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; - Tag tag; - return squared_distance(pt, seg, k, tag); + 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(); + + // 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 diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Triangle_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Triangle_3.h index 8859831c6f0..a6f2e7c8e94 100644 --- a/Distance_3/include/CGAL/Distance_3/Point_3_Triangle_3.h +++ b/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 +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 +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 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); diff --git a/Distance_3/include/CGAL/Distance_3/internal/squared_distance_utils_3.h b/Distance_3/include/CGAL/Distance_3/internal/squared_distance_utils_3.h index b9669d86eb7..7dbf82e246e 100644 --- a/Distance_3/include/CGAL/Distance_3/internal/squared_distance_utils_3.h +++ b/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 void -squared_distance_to_plane_RT(const typename K::Vector_3 & normal, - const typename K::Vector_3 & diff, +squared_distance_to_plane_RT(const typename K::Vector_3& normal, + 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 typename K::FT -squared_distance_to_plane(const typename K::Vector_3 & normal, - const typename K::Vector_3 & diff, +squared_distance_to_plane(const typename K::Vector_3& normal, + 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 void -squared_distance_to_line_RT(const typename K::Vector_3 & dir, - const typename K::Vector_3 & diff, +squared_distance_to_line_RT(const typename K::Vector_3& dir, + 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 typename K::FT -squared_distance_to_line(const typename K::Vector_3 & dir, - const typename K::Vector_3 & diff, +squared_distance_to_line(const typename K::Vector_3& dir, + const typename K::Vector_3& diff, const K& k) { typedef typename K::RT RT;