diff --git a/Distance_3/include/CGAL/squared_distance_3_0.h b/Distance_3/include/CGAL/squared_distance_3_0.h index 7b496c5967a..835c2b97b66 100644 --- a/Distance_3/include/CGAL/squared_distance_3_0.h +++ b/Distance_3/include/CGAL/squared_distance_3_0.h @@ -26,6 +26,7 @@ #include #include #include +#include namespace CGAL { @@ -187,6 +188,21 @@ squared_distance(const typename K::Point_3 & pt1, return k.compute_squared_distance_3_object()(pt1, pt2); } +template +void +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) +{ + typedef typename K::RT RT; + RT dot, squared_length; + dot = wdot(normal, diff, k); + squared_length = wdot(normal, normal, k); + num = square(dot); + den = wmult((K*)0, squared_length, diff.hw(), diff.hw()); +} template typename K::FT @@ -196,13 +212,24 @@ squared_distance_to_plane(const typename K::Vector_3 & normal, { typedef typename K::RT RT; typedef typename K::FT FT; - RT dot, squared_length; - dot = wdot(normal, diff, k); - squared_length = wdot(normal, normal, k); - return FT(dot*dot) / - FT(wmult((K*)0, squared_length, diff.hw(), diff.hw())); + RT num, den; + squared_distance_to_plane_RT(normal, diff, num, den, k); + return Rational_traits().make_rational(num, den); } +template +void +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) +{ + typedef typename K::Vector_3 Vector_3; + Vector_3 wcr = wcross(dir, diff, k); + num = wdot(wcr, wcr, k); + den = wmult((K*)0, wdot(dir, dir, k), diff.hw(), diff.hw()); +} template typename K::FT @@ -210,15 +237,13 @@ squared_distance_to_line(const typename K::Vector_3 & dir, const typename K::Vector_3 & diff, const K& k) { - typedef typename K::Vector_3 Vector_3; typedef typename K::RT RT; typedef typename K::FT FT; - Vector_3 wcr = wcross(dir, diff, k); - return FT(wcr*wcr)/FT(wmult( - (K*)0, RT(wdot(dir, dir, k)), diff.hw(), diff.hw())); + RT num, den; + squared_distance_to_line_RT(dir, diff, num, den, k); + return Rational_traits().make_rational(num, den); } - template inline bool diff --git a/Distance_3/include/CGAL/squared_distance_3_1.h b/Distance_3/include/CGAL/squared_distance_3_1.h index e0e63a6fb23..fdc02382dca 100644 --- a/Distance_3/include/CGAL/squared_distance_3_1.h +++ b/Distance_3/include/CGAL/squared_distance_3_1.h @@ -31,17 +31,33 @@ namespace CGAL { namespace internal { template -typename K::FT -squared_distance( +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 construct_vector; Vector_3 dir(line.direction().vector()); Vector_3 diff = construct_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 @@ -56,6 +72,31 @@ squared_distance( } +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) +{ + typename K::Construct_vector_3 construct_vector; + typedef typename K::Vector_3 Vector_3; + typedef typename K::RT RT; + + Vector_3 diff = construct_vector(ray.source(), pt); + const Vector_3 &dir = ray.direction().vector(); + 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( @@ -63,16 +104,25 @@ squared_distance( const typename K::Ray_3 &ray, const K& k) { - typename K::Construct_vector_3 construct_vector; + // 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; - Vector_3 diff = construct_vector(ray.source(), pt); - const Vector_3 &dir = ray.direction().vector(); - if (!is_acute_angle(dir,diff, k) ) - return (typename K::FT)(diff*diff); - return squared_distance_to_line(dir, diff, k); -} + typename K::Construct_vector_3 construct_vector = k.construct_vector_3_object(); + Vector_3 diff = construct_vector(ray.source(), pt); + const Vector_3 &dir = ray.direction().vector(); + if (!is_acute_angle(dir,diff, k) ) + return (typename K::FT)(diff*diff); + + return squared_distance_to_line(dir, diff, k); +} template inline @@ -85,75 +135,77 @@ squared_distance( return squared_distance(pt, ray, k); } - - - template -typename K::FT -squared_distance( +void +squared_distance_RT( const typename K::Point_3 &pt, const typename K::Segment_3 &seg, - const K& k, - const Homogeneous_tag) + typename K::RT& num, + typename K::RT& den, + const K& k) { - typename K::Construct_vector_3 construct_vector; - typedef typename K::Vector_3 Vector_3; - typedef typename K::RT RT; - typedef typename K::FT FT; - // assert that the segment is valid (non zero length). - Vector_3 diff = construct_vector(seg.source(), pt); - Vector_3 segvec = construct_vector(seg.source(), seg.target()); - RT d = wdot(diff,segvec, k); - if (d <= (RT)0) - return (FT(diff*diff)); - RT e = wdot(segvec,segvec, k); - if ( (d * segvec.hw()) > (e * diff.hw())) - return squared_distance(pt, seg.target(), k); + typedef typename K::Vector_3 Vector_3; + typedef typename K::RT RT; - Vector_3 wcr = wcross(segvec, diff, k); - return FT(wcr*wcr)/FT(e * diff.hw() * diff.hw()); + typename K::Construct_vector_3 construct_vector; + + // assert that the segment is valid (non zero length). + const Vector_3 diff_s = construct_vector(seg.source(), pt); + const Vector_3 segvec = construct_vector(seg.source(), seg.target()); + + const RT d = wdot(diff_s, segvec, k); + if (d <= RT(0)) + { + // 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; + } + + 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 = construct_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&) -{ - typename K::Construct_vector_3 construct_vector; - typedef typename K::Vector_3 Vector_3; - typedef typename K::RT RT; - typedef typename K::FT FT; - // assert that the segment is valid (non zero length). - Vector_3 diff = construct_vector(seg.source(), pt); - Vector_3 segvec = construct_vector(seg.source(), seg.target()); - RT d = wdot(diff,segvec, k); - if (d <= (RT)0) - return (FT(diff*diff)); - RT e = wdot(segvec,segvec, k); - if (d > e) - return squared_distance(pt, seg.target(), k); - - 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::Vector_3 Vector_3; + typedef typename K::RT RT; + typedef typename K::FT FT; + typename K::Construct_vector_3 construct_vector = k.construct_vector_3_object(); + + // assert that the segment is valid (non zero length). + Vector_3 diff = construct_vector(seg.source(), pt); + Vector_3 segvec = construct_vector(seg.source(), seg.target()); + RT d = wdot(diff,segvec, k); + if (d <= RT(0)) + return (FT(diff*diff)); + + 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' + Vector_3 wcr = wcross(segvec, diff, k); + return FT(wcr*wcr) / wmult((K*)0, e, diff.hw(), diff.hw()); +} template inline @@ -166,8 +218,6 @@ squared_distance( return squared_distance(pt, seg, k); } - - template typename K::Comparison_result compare_distance_pssC3( @@ -187,7 +237,7 @@ compare_distance_pssC3( Vector_3 diff = construct_vector(seg1.source(), pt); Vector_3 segvec = construct_vector(seg1.source(), seg1.target()); RT d = wdot(diff,segvec, k); - if (d <= (RT)0){ + if (d <= RT(0)){ d1 = (FT(diff*diff)); }else{ RT e = wdot(segvec,segvec, k); @@ -204,7 +254,7 @@ compare_distance_pssC3( Vector_3 diff = construct_vector(seg2.source(), pt); Vector_3 segvec = construct_vector(seg2.source(), seg2.target()); RT d = wdot(diff,segvec, k); - if (d <= (RT)0){ + if (d <= RT(0)){ d2 = (FT(diff*diff)); }else{ RT e = wdot(segvec,segvec, k); @@ -240,7 +290,7 @@ compare_distance_ppsC3( Vector_3 diff = construct_vector(seg.source(), pt); Vector_3 segvec = construct_vector(seg.source(), seg.target()); RT d = wdot(diff,segvec, k); - if (d <= (RT)0){ + if (d <= RT(0)){ d2 = (FT(diff*diff)); }else{ RT e = wdot(segvec,segvec, k); diff --git a/Distance_3/include/CGAL/squared_distance_3_2.h b/Distance_3/include/CGAL/squared_distance_3_2.h index 78af20db08f..c5cd5cf1694 100644 --- a/Distance_3/include/CGAL/squared_distance_3_2.h +++ b/Distance_3/include/CGAL/squared_distance_3_2.h @@ -218,6 +218,91 @@ on_left_of_triangle_edge(const typename K::Point_3 & pt, return result; } +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) +{ + typename K::Construct_vector_3 vector; + typedef typename K::Vector_3 Vector_3; + 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, typename K::Segment_3(t2, t0), num, den, k); + + typename K::RT num2, den2; + squared_distance_RT(pt, typename K::Segment_3(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 + squared_distance_RT(pt, typename K::Segment_3(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, typename K::Segment_3(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, typename K::Segment_3(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, typename K::Segment_3(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( @@ -234,31 +319,37 @@ squared_distance_to_triangle( const Vector_3 oe3 = vector(t0, t2); const Vector_3 normal = wcross(e1, oe3, k); - if(normal != NULL_VECTOR - && on_left_of_triangle_edge(pt, normal, t0, t1, k) - && on_left_of_triangle_edge(pt, normal, t1, t2, k) - && on_left_of_triangle_edge(pt, normal, t2, t0, k)) - { - // the projection of pt is inside the triangle - inside = true; - return squared_distance_to_plane(normal, vector(t0, pt), k); - } - else { - // 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. - typename K::FT d1 = squared_distance(pt, - typename K::Segment_3(t2, t0), - k); - typename K::FT d2 = squared_distance(pt, - typename K::Segment_3(t1, t2), - k); - typename K::FT d3 = squared_distance(pt, - typename K::Segment_3(t0, t1), - 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. + typename K::FT d1 = squared_distance(pt, typename K::Segment_3(t2, t0), k); + typename K::FT d2 = squared_distance(pt, typename K::Segment_3(t1, t2), k); - return (std::min)( (std::min)(d1, d2), d3); - } + // should not be needed since at most 2 edges cover the full triangle in the degenerate case + typename K::FT d3 = squared_distance(pt, typename K::Segment_3(t0, t1), k); + + return (std::min)( (std::min)(d1, d2), d3); + } + + const bool b01 = on_left_of_triangle_edge(pt, normal, t0, t1, k); + if(!b01) { + return internal::squared_distance(pt, typename K::Segment_3(t0, t1), k); + } + + const bool b12 = on_left_of_triangle_edge(pt, normal, t1, t2, k); + if(!b12) { + return internal::squared_distance(pt, typename K::Segment_3(t1, t2), k); + } + + const bool b20 = on_left_of_triangle_edge(pt, normal, t2, t0, k); + if(!b20) { + return internal::squared_distance(pt, typename K::Segment_3(t2, t0), k); + } + + // The projection of pt is inside the triangle + inside = true; + return squared_distance_to_plane(normal, vector(t0, pt), k); } template @@ -268,7 +359,7 @@ squared_distance( const typename K::Triangle_3 & t, const K& k) { - typename K::Construct_vertex_3 vertex; + typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object(); bool inside = false; return squared_distance_to_triangle(pt, vertex(t, 0), diff --git a/Distance_3/package_info/Distance_3/dependencies b/Distance_3/package_info/Distance_3/dependencies index 2372c4ab8bb..d8f0e86a70c 100644 --- a/Distance_3/package_info/Distance_3/dependencies +++ b/Distance_3/package_info/Distance_3/dependencies @@ -2,7 +2,9 @@ Algebraic_foundations Distance_2 Distance_3 Installation +Interval_support Kernel_23 +Modular_arithmetic Number_types Profiling_tools STL_Extension diff --git a/Distance_3/test/Distance_3/test_distance_3.cpp b/Distance_3/test/Distance_3/test_distance_3.cpp index 3a1f72fe786..0632b699765 100644 --- a/Distance_3/test/Distance_3/test_distance_3.cpp +++ b/Distance_3/test/Distance_3/test_distance_3.cpp @@ -112,6 +112,17 @@ struct Test { { std::cout << "Point - Triangle\n"; check_squared_distance (p(0, 1, 2), T(p(0, 0, 0), p( 2, 0, 0), p( 0, 2, 0)), 4); + + T t(p(0,0,0), p(3,0,0), p(3,3,0)); + check_squared_distance (p(-1, -1, 0), t, 2); + check_squared_distance (p(-1, 0, 0), t, 1); + check_squared_distance (p(0, 0, 0), t, 0); + check_squared_distance (p(1, 0, 0), t, 0); + check_squared_distance (p(4, 0, 0), t, 1); + check_squared_distance (p(1, -1, 0), t, 1); + check_squared_distance (p(1, 1, 1), t, 1); + check_squared_distance (p(1, 0, 1), t, 1); + check_squared_distance (p(0, 0, 1), t, 1); } void P_Tet() diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h index f8fa015669a..e35346a9d5a 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h @@ -15,7 +15,7 @@ #include #include - +#include namespace CGAL { template @@ -36,9 +36,15 @@ inline typename K::Boolean do_intersect(const typename K::Sphere_3 &sp, const typename K::Triangle_3 &tr, - const K & /* k */) + const K & k) { - return squared_distance(sp.center(), tr) <= sp.squared_radius(); + typedef typename K::RT RT; + RT num, den; + + CGAL::internal::squared_distance_RT(sp.center(), tr, num, den, k); + return ! (compare_quotients(num, den, + Rational_traits().numerator(sp.squared_radius()), + Rational_traits().denominator(sp.squared_radius())) == LARGER); } template @@ -46,18 +52,9 @@ inline typename K::Boolean do_intersect(const typename K::Triangle_3 &tr, const typename K::Sphere_3 &sp, - const K & /* k */) + const K & k) { - return squared_distance(sp.center(), tr) <= sp.squared_radius(); -} -template -inline -typename K::Boolean -do_intersect(const typename K::Sphere_3 &sp, - const typename K::Line_3 &lin, - const K & /* k */) -{ - return squared_distance(sp.center(), lin) <= sp.squared_radius(); + return do_intersect(sp, tr, k); } @@ -66,9 +63,26 @@ inline typename K::Boolean do_intersect(const typename K::Line_3 &lin, const typename K::Sphere_3 &sp, - const K & /* k */) + const K & k) { - return squared_distance(sp.center(), lin) <= sp.squared_radius(); + typedef typename K::RT RT; + RT num, den; + + CGAL::internal::squared_distance_RT(sp.center(), lin, num, den, k); + return ! (compare_quotients(num, den, + Rational_traits().numerator(sp.squared_radius()), + Rational_traits().denominator(sp.squared_radius())) == LARGER); + +} + +template +inline +typename K::Boolean +do_intersect(const typename K::Sphere_3 &sp, + const typename K::Line_3 &lin, + const K & k) +{ + return do_intersect(lin,sp,k); } @@ -78,9 +92,15 @@ inline typename K::Boolean do_intersect(const typename K::Sphere_3 &sp, const typename K::Ray_3 &ray, - const K & /* k */) + const K & k) { - return squared_distance(sp.center(), ray) <= sp.squared_radius(); + typedef typename K::RT RT; + RT num, den; + + CGAL::internal::squared_distance_RT(sp.center(), ray, num, den, k); + return ! (compare_quotients(num, den, + Rational_traits().numerator(sp.squared_radius()), + Rational_traits().denominator(sp.squared_radius())) == LARGER); } @@ -89,9 +109,9 @@ inline typename K::Boolean do_intersect(const typename K::Ray_3 &ray, const typename K::Sphere_3 &sp, - const K & /* k */) + const K & k) { - return squared_distance(sp.center(), ray) <= sp.squared_radius(); + return do_intersect(sp,ray,k); } template @@ -99,9 +119,15 @@ inline typename K::Boolean do_intersect(const typename K::Sphere_3 &sp, const typename K::Segment_3 &seg, - const K & /* k */) + const K & k) { - return squared_distance(sp.center(), seg) <= sp.squared_radius(); + typedef typename K::RT RT; + RT num, den; + + CGAL::internal::squared_distance_RT(sp.center(), seg, num, den, k); + return ! (compare_quotients(num, den, + Rational_traits().numerator(sp.squared_radius()), + Rational_traits().denominator(sp.squared_radius())) == LARGER); } @@ -110,9 +136,9 @@ inline typename K::Boolean do_intersect(const typename K::Segment_3 &seg, const typename K::Sphere_3 &sp, - const K & /* k */) + const K & k) { - return squared_distance(sp.center(), seg) <= sp.squared_radius(); + return do_intersect(sp,seg,k); } } // namespace internal diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp index b8953c92b8b..d407f6c313d 100644 --- a/Intersections_3/test/Intersections_3/test_intersections_3.cpp +++ b/Intersections_3/test/Intersections_3/test_intersections_3.cpp @@ -161,13 +161,35 @@ struct Test { void P_do_intersect() { - P p(0,0,0), q(1,0,0), r(2,0,0), s(10,10,10); + P p(0,0,0), q(1,0,0), r(2,0,0), s(10,10,10), s2(10,10,0), s3(10,-10,-1), t(3,0,0); Sph sph(p,1); Cub cub(p,r); assert(do_intersect(q,sph)); assert(do_intersect(sph,q)); assert(! do_intersect(s,cub)); assert(! do_intersect(cub,s)); + + assert(do_intersect(sph, Tr(p,q,s))); + assert(do_intersect(sph, Tr(r,q,s))); + assert(! do_intersect(sph, Tr(r,s,t))); + + assert(! do_intersect(sph, Tr(s,s2,s3))); + + assert(do_intersect(sph, S(p,r))); + assert(do_intersect(sph, S(q,r))); + assert(do_intersect(sph, S(s,q))); + assert(! do_intersect(sph, S(s,r))); + + assert(do_intersect(sph, L(p,r))); + assert(do_intersect(sph, L(q,r))); + assert(do_intersect(sph, L(s,q))); + assert(! do_intersect(sph, L(s,r))); + + assert(do_intersect(sph, R(p,r))); + assert(do_intersect(sph, R(q,r))); + assert(do_intersect(sph, R(s,q))); + assert(! do_intersect(sph, R(s,r))); + } @@ -1265,6 +1287,7 @@ struct Test { Bbox_L(); Bbox_R(); Bbox_Tr(); + } };