diff --git a/Distance_3/include/CGAL/Distance_3/Line_3_Line_3.h b/Distance_3/include/CGAL/Distance_3/Line_3_Line_3.h index 85f181af18a..a27c6b7aaa6 100644 --- a/Distance_3/include/CGAL/Distance_3/Line_3_Line_3.h +++ b/Distance_3/include/CGAL/Distance_3/Line_3_Line_3.h @@ -57,26 +57,6 @@ compare_squared_distance(const typename K::Line_3& line1, } // namespace internal -template -inline -typename K::FT -squared_distance(const Line_3& line1, - const Line_3& line2) -{ - return K().compute_squared_distance_3_object()(line1, line2); -} - -template -inline -typename K::Comparison_result -compare_squared_distance(const Line_3& line1, - const Line_3& line2, - const typename K::FT& d2) -{ - return K().compare_squared_distance_3_object()(line1,line2,d2); -} - - } // namespace CGAL #endif // CGAL_DISTANCE_3_LINE_3_LINE_3_H diff --git a/Distance_3/include/CGAL/Distance_3/Line_3_Plane_3.h b/Distance_3/include/CGAL/Distance_3/Line_3_Plane_3.h index 6bb4d125c59..df4a14ef8ed 100644 --- a/Distance_3/include/CGAL/Distance_3/Line_3_Plane_3.h +++ b/Distance_3/include/CGAL/Distance_3/Line_3_Plane_3.h @@ -61,26 +61,28 @@ squared_distance(const typename K::Plane_3& pl, return squared_distance(l, pl, k); } +template +inline typename K::Comparison_result +compare_squared_distance(const typename K::Line_3& l, + const typename K::Plane_3& pl, + const K& k, + const typename K::FT& d2) +{ + return compare(squared_distance(l, pl, k), d2); +} + +template +inline typename K::Comparison_result +compare_squared_distance(const typename K::Plane_3& pl, + const typename K::Line_3& l, + const K& k, + const typename K::FT& d2) +{ + return compare_squared_distance(l, pl, k, d2); +} + } // namespace internal -template -inline -typename K::FT -squared_distance(const Line_3& line, - const Plane_3& plane) -{ - return K().compute_squared_distance_3_object()(line, plane); -} - -template -inline -typename K::FT -squared_distance(const Plane_3& plane, - const Line_3& line) -{ - return K().compute_squared_distance_3_object()(plane, line); -} - } // namespace CGAL #endif // CGAL_DISTANCE_3_LINE_3_PLANE_3_H diff --git a/Distance_3/include/CGAL/Distance_3/Plane_3_Plane_3.h b/Distance_3/include/CGAL/Distance_3/Plane_3_Plane_3.h index d18fb366b7f..4ddb7ad1702 100644 --- a/Distance_3/include/CGAL/Distance_3/Plane_3_Plane_3.h +++ b/Distance_3/include/CGAL/Distance_3/Plane_3_Plane_3.h @@ -39,17 +39,18 @@ squared_distance(const typename K::Plane_3& plane1, return sq_dist(plane1.point(), plane2); } -} // namespace internal - template -inline -typename K::FT -squared_distance(const Plane_3& plane1, - const Plane_3& plane2) +inline typename K::Comparison_result +compare_squared_distance(const typename K::Plane_3& plane1, + const typename K::Plane_3& plane2, + const K& k, + const typename K::FT& d2) { - return K().compute_squared_distance_3_object()(plane1, plane2); + return compare(squared_distance(plane1,plane2,k), d2); } +} // namespace internal + } // namespace CGAL #endif // CGAL_DISTANCE_3_PLANE_3_PLANE_3_H 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 e202cb3bdd3..b265894c910 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 @@ -85,49 +85,11 @@ compare_squared_distance(const typename K::Line_3& line, const K& k, const typename K::FT& d2) { - return compare(squared_distance(pt, line, k), d2); + return compare_squared_distance(pt, line, k, d2); } } // namespace internal -template -inline -typename K::FT -squared_distance(const Point_3& pt, - const Line_3& line) -{ - return K().compute_squared_distance_3_object()(pt, line); -} - -template -inline -typename K::FT -squared_distance(const Line_3& line, - const Point_3& pt) -{ - return K().compute_squared_distance_3_object()(line, pt); -} - -template -inline -typename K::Comparison_result -compare_squared_distance(const Point_3& pt, - const Line_3& line, - const typename K::FT& d2) -{ - return K().compare_squared_distance_3_object()(pt, line, d2); -} - -template -inline -typename K::Comparison_result -compare_squared_distance(const Line_3& line, - const Point_3& pt, - const typename K::FT& d2) -{ - return K().compare_squared_distance_3_object()(line, pt, d2); -} - } // namespace CGAL #endif // CGAL_DISTANCE_3_POINT_3_LINE_3_H diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Plane_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Plane_3.h index 41107dfac3e..921b13a446f 100644 --- a/Distance_3/include/CGAL/Distance_3/Point_3_Plane_3.h +++ b/Distance_3/include/CGAL/Distance_3/Point_3_Plane_3.h @@ -68,47 +68,6 @@ compare_squared_distance(const typename K::Plane_3& plane, return compare_squared_distance(pt, plane, k, d2); } -} // namespace internal - -template -inline -typename K::FT -squared_distance(const Point_3& pt, - const Plane_3& plane) -{ - return K().compute_squared_distance_3_object()(pt, plane); -} - -template -inline -typename K::FT -squared_distance(const Plane_3& plane, - const Point_3& pt) -{ - return K().compute_squared_distance_3_object()(plane, pt); -} - -template -inline -typename K::Comparison_result -compare_squared_distance(const Point_3& pt, - const Plane_3& pl, - const typename K::FT& d2) -{ - return K().compare_squared_distance_3_object()(pt, pl, d2); -} - -template -inline -typename K::Comparison_result -compare_squared_distance(const Plane_3& pl, - const Point_3& pt, - const typename K::FT& d2) -{ - return K().compare_squared_distance_3_object()(pl, pt, d2); -} - - -} // namespace CGAL +} } // namespace CGAL::internal #endif // CGAL_DISTANCE_3_POINT_3_PLANE_3_H diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Point_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Point_3.h index ec5102e8bee..4d7de336e31 100644 --- a/Distance_3/include/CGAL/Distance_3/Point_3_Point_3.h +++ b/Distance_3/include/CGAL/Distance_3/Point_3_Point_3.h @@ -45,15 +45,6 @@ compare_squared_distance(const typename K::Point_3& pt1, } // namespace internal -template -inline -typename K::FT -squared_distance(const Point_3& pt1, - const Point_3& pt2) -{ - return internal::squared_distance(pt1, pt2, K()); -} - } // namespace CGAL #endif // CGAL_DISTANCE_3_POINT_3_POINT_3_H 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 71332a693b5..cb2bdd7401d 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 @@ -101,6 +101,7 @@ compare_squared_distance(const typename K::Point_3& pt, const Vector_3 dir = ray.direction().vector(); const Vector_3 diff = vector(ray.source(), pt); + //Compare first the distance to the line, if larger we can exit early const typename K::Comparison_result res_pl = compare_squared_distance_to_line(dir, diff, k, d2); if(res_pl==LARGER) return LARGER; @@ -123,45 +124,6 @@ compare_squared_distance(const typename K::Ray_3& ray, } // namespace internal -template -inline -typename K::FT -squared_distance(const Point_3& pt, - const Ray_3& ray) -{ - return K().compute_squared_distance_3_object()(pt, ray); -} - -template -inline -typename K::FT -squared_distance(const Ray_3& ray, - const Point_3& pt) -{ - return K().compute_squared_distance_3_object()(ray, pt); -} - -template -inline -typename K::Comparison_result -compare_squared_distance(const Point_3& pt, - const Ray_3& ray, - const typename K::FT& d2) -{ - return K().compare_squared_distance_3_object()(pt, ray, d2); -} - -template -inline -typename K::Comparison_result -compare_squared_distance(const Ray_3& ray, - const Point_3& pt, - const typename K::FT& d2) -{ - return K().compare_squared_distance_3_object()(ray, pt, d2); -} - - } // namespace CGAL #endif // CGAL_DISTANCE_3_POINT_3_RAY_3_H 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 5188a13622a..1a94bee31f5 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 @@ -79,6 +79,7 @@ squared_distance(const typename K::Point_3& pt, 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). const Vector_3 diff = vector(seg.source(), pt); @@ -90,7 +91,7 @@ squared_distance(const typename K::Point_3& pt, 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); + return sq_dist(pt, seg.target()); // This is an expanded call to squared_distance_to_line() to avoid recomputing 'e' const Vector_3 wcr = wcross(segvec, diff, k); @@ -113,37 +114,35 @@ typename K::Comparison_result compare_squared_distance(const typename K::Point_3& pt, const typename K::Segment_3& seg, const K& k, - const typename K::FT &d2) + const typename K::FT& d2) { 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::Compare_squared_distance_3 csq_dist = k.compare_squared_distance_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()); - //We first compare the distance of the point and the line + //Compare first the distance to the line, if larger we can exit early const typename K::Comparison_result res_pl= compare_squared_distance_to_line(segvec, diff, k, d2); - - //If greater than d2, we early exit if(res_pl==LARGER) return LARGER; - //If distance is realized by the source + //If the distance is realized by the source const RT d = wdot(diff, segvec, k); if(d <= RT(0)) return compare(FT(diff*diff), d2); - //If distance is realized by the target + //If the distance is realized by the target const RT e = wdot(segvec, segvec, k); if(wmult((K*)0, d, segvec.hw()) > wmult((K*)0, e, diff.hw())) - return compare_squared_distance(pt, seg.target(), k, d2); + return csq_dist(pt, seg.target(), d2); - //If distance is realized by the interior, it is equal to the one of the line + //If the distance is realized by the interior return res_pl; } @@ -160,44 +159,6 @@ compare_squared_distance(const typename K::Segment_3& seg, } // namespace internal -template -inline -typename K::FT -squared_distance(const Point_3& pt, - const Segment_3& seg) -{ - return K().compute_squared_distance_3_object()(pt, seg); -} - -template -inline -typename K::FT -squared_distance(const Segment_3& seg, - const Point_3& pt) -{ - return K().compute_squared_distance_3_object()(seg, pt); -} - -template -inline -typename K::Comparison_result -compare_squared_distance(const Point_3& pt, - const Segment_3& seg, - const typename K::FT &d2) -{ - return K().compare_squared_distance_3_object()(pt, seg, d2); -} - -template -inline -typename K::Comparison_result -compare_squared_distance(const Segment_3& seg, - const Point_3& pt, - const typename K::FT &d2) -{ - return K().compare_squared_distance_3_object()(seg, pt, d2); -} - } // namespace CGAL #endif // CGAL_DISTANCE_3_POINT_3_SEGMENT_3_H diff --git a/Distance_3/include/CGAL/Distance_3/Segment_3_Segment_3.h b/Distance_3/include/CGAL/Distance_3/Segment_3_Segment_3.h index e993eb94ba8..5c9562c53bb 100644 --- a/Distance_3/include/CGAL/Distance_3/Segment_3_Segment_3.h +++ b/Distance_3/include/CGAL/Distance_3/Segment_3_Segment_3.h @@ -140,97 +140,6 @@ squared_distance(const typename K::Segment_3& s1, return res; } -// Using Lumelsky, 'On Fast Computation of Distance Between Line Segments' 1984 -template -typename K::Comparison_result -compare_squared_distance(const typename K::Segment_3& s1, - const typename K::Segment_3& s2, - const K& k, - const typename K::FT& d2) -{ - typedef typename K::FT FT; - typedef typename K::Point_3 Point_3; - typedef typename K::Vector_3 Vector_3; - - typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object(); - typename K::Construct_vector_3 cv = k.construct_vector_3_object(); - typename K::Compute_scalar_product_3 sp = k.compute_scalar_product_3_object(); - typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object(); - - const Point_3& p1 = vertex(s1, 0); - const Point_3& q1 = vertex(s1, 1); - const Point_3& p2 = vertex(s2, 0); - const Point_3& q2 = vertex(s2, 1); - const Vector_3 v1 = cv(p1, q1), v2 = cv(p2, q2); - const Vector_3 p1p2 = cv(p1, p2); - - //If degenerate segment, compare distance between the point and the other segment - if(p1 == q1) - if(p2 == q2) - return compare_square_distance(p1,p2,k,d2); - else - return compare_square_distance(p1,s2,k,d2); - else if(p2 == q2) - return compare_square_distance(s1,p2,k,d2); - - //Compare distance between the lines - const Vector_3 normal = wcross(v1, v2, k); - const Vector_3 diff = p1p2; - - //If the lines are at distance more than d, the segment do not be at distance less than d - typename K::Comparison_result res_ll=compare_squared_distance(v1.supporting_line(), v2.supporting_line(), k, d2); - if(res_ll==LARGER) - return LARGER; //Early exit - - //Compute distance between the segments - //Since we already compare the distance between the lines, we care about only the cases where - //the distance are realized by vertices aka res.x or res.y are 0 or 1 - - const FT a = sp(v1, v1); - const FT b = - sp(v1, v2); - const FT c = - b; - const FT d = - sp(v2, v2); - const FT e = sp(v1, p1p2); - const FT f = sp(v2, p1p2); - - CGAL_assertion(a > 0 && d < 0); - const FT det = a*d - b*c; - const FT res_x; - if(det == 0) - res_x = 0; - else - res_x = boost::algorithm::clamp((e*d - b*f) / det, 0, 1); - - FT xc = res_x*c; - // res.y = (f - xc) / d, by definition, but building it up more efficiently - if(f > xc) // y < 0 <=> f - xc / d < 0 <=> f - xc > 0 (since d is -||v2||) - { - // res_y = 0; - res_x = boost::algorithm::clamp(e/a, 0, 1); // (e + y*c) / a - return compare(sq_dist(p1+res_x*v1,p2), d2); - } - else // y >= 0 - { - FT n = f - xc; // delay the division by d - if(n < d) // y > 1 <=> n / d > 1 <=> n < d (once again, important to note that d is negative!) - { - // res_y = 1; - res_x = boost::algorithm::clamp((e + c) / a, 0, 1); // (e + y*c) / a - return compare(sq_dist(p1+res_x*v1,p2), d2); - } - else if(res_x==0 || res_x==1) // 0 <= y <= 1 - { - FT res_y = n / d; - return compare(sq_dist(p1 + res_x*v1, p2 + res_y*v2), d2); - } - else - { - //Already computed by distance line line - return res_ll; - } - } -} - } // namespace internal } // namespace Distance_3 @@ -279,26 +188,97 @@ squared_distance(const typename K::Segment_3& seg1, return res.squared_distance; } +template +typename K::Comparison_result +compare_squared_distance(const typename K::Segment_3& s1, + const typename K::Segment_3& s2, + const K& k, + const typename K::FT& d2) +{ + typedef typename K::FT FT; + typedef typename K::Point_3 Point_3; + typedef typename K::Vector_3 Vector_3; + + typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object(); + typename K::Construct_vector_3 cv = k.construct_vector_3_object(); + typename K::Compute_scalar_product_3 sp = k.compute_scalar_product_3_object(); + typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object(); + typename K::Compare_squared_distance_3 csq_dist = k.compare_squared_distance_3_object(); + + const Point_3& p1 = vertex(s1, 0); + const Point_3& q1 = vertex(s1, 1); + const Point_3& p2 = vertex(s2, 0); + const Point_3& q2 = vertex(s2, 1); + const Vector_3 v1 = cv(p1, q1), v2 = cv(p2, q2); + const Vector_3 p1p2 = cv(p1, p2); + + // If degenerate segment, compare the distance between the point and the other segment + if(p1 == q1) + if(p2 == q2) + return csq_dist(p1,p2,d2); + else + return csq_dist(p1,s2,d2); + else if(p2 == q2) + return csq_dist(s1,p2,d2); + + const Vector_3 normal = wcross(v1, v2, k); + const Vector_3 diff = p1p2; + + // Compare first the distance between the lines, if larger we can exit early + typename K::Comparison_result res_ll=csq_dist(s1.supporting_line(), s2.supporting_line(), d2); + if(is_certain(res_ll) && res_ll==LARGER) + return LARGER; + + // Compute the distance between the segments + // TODO some factorization with above function? only the last case is different + + const FT a = sp(v1, v1); + const FT b = - sp(v1, v2); + const FT c = - b; + const FT d = - sp(v2, v2); + const FT e = sp(v1, p1p2); + const FT f = sp(v2, p1p2); + + CGAL_assertion(a > 0 && d < 0); + const FT det = a*d - b*c; + FT res_x; + if(det == 0) + res_x = 0; + else + res_x = boost::algorithm::clamp((e*d - b*f) / det, 0, 1); + + FT xc = res_x*c; + // res.y = (f - xc) / d, by definition, but building it up more efficiently + if(f > xc) // y < 0 <=> f - xc / d < 0 <=> f - xc > 0 (since d is -||v2||) + { + // res_y = 0; + res_x = boost::algorithm::clamp(e/a, 0, 1); // (e + y*c) / a + return compare(sq_dist(p1+res_x*v1,p2), d2); + } + else // y >= 0 + { + FT n = f - xc; // delay the division by d + if(n < d) // y > 1 <=> n / d > 1 <=> n < d (once again, important to note that d is negative!) + { + // res_y = 1; + res_x = boost::algorithm::clamp((e + c) / a, 0, 1); // (e + y*c) / a + return compare(sq_dist(p1+res_x*v1,p2), d2); + } + else if(res_x==0 || res_x==1) // 0 <= y <= 1 + { + FT res_y = n / d; + return compare(sq_dist(p1 + res_x*v1, p2 + res_y*v2), d2); + } + else + { + //Already computed by distance line line + return res_ll; + } + } +} + } // namespace internal -template -inline -typename K::FT -squared_distance_parallel(const Segment_3& seg1, - const Segment_3& seg2) -{ - return internal::squared_distance_parallel(seg1, seg2, K()); -} - -template -inline -typename K::FT -squared_distance(const Segment_3& seg1, - const Segment_3& seg2) -{ - return K().compute_squared_distance_3_object()(seg1, seg2); -} - } // namespace CGAL #endif // CGAL_DISTANCE_3_SEGMENT_3_SEGMENT_3_H diff --git a/Distance_3/include/CGAL/Distance_3/Triangle_3_Triangle_3.h b/Distance_3/include/CGAL/Distance_3/Triangle_3_Triangle_3.h index f5e66c37d21..2d0106bf9e9 100644 --- a/Distance_3/include/CGAL/Distance_3/Triangle_3_Triangle_3.h +++ b/Distance_3/include/CGAL/Distance_3/Triangle_3_Triangle_3.h @@ -217,8 +217,8 @@ template typename K::Comparison_result compare_squared_distance(const typename K::Triangle_3& tr1, const typename K::Triangle_3& tr2, - const typename K::FT& d, - const K& k) + const K& k, + const typename K::FT& d2) { typedef typename K::FT FT; @@ -237,15 +237,15 @@ compare_squared_distance(const typename K::Triangle_3& tr1, { for(int j=0; j<3; ++j) { - if(CGAL::possibly(CGAL::compare_squared_distance(Line_3(vertex(tr1, i%3), vertex(tr1, (i+1)%3)),Line_3(vertex(tr2, j%3), vertex(tr2, (j+1)%3)),d)!=CGAL::LARGER)) + if(CGAL::possibly(CGAL::compare_squared_distance(Line_3(vertex(tr1, i%3), vertex(tr1, (i+1)%3)),Line_3(vertex(tr2, j%3), vertex(tr2, (j+1)%3)),d2)!=CGAL::LARGER)) { - typename K::comparison_result res_seg_seg= CGAL::compare_squared_distance(Segment_3(vertex(tr1, i%3), vertex(tr1, (i+1)%3)),Segment_3(vertex(tr2, j%3), vertex(tr2, (j+1)%3)),d); + typename K::comparison_result res_seg_seg= CGAL::compare_squared_distance(Segment_3(vertex(tr1, i%3), vertex(tr1, (i+1)%3)),Segment_3(vertex(tr2, j%3), vertex(tr2, (j+1)%3)),d2); if(is_certain(res_seg_seg) && res_seg_seg==CGAL::SMALLER) return CGAL::SMALLER; } } - typename K::comparison_result res_v_pl= CGAL::compare_squared_distance(vertex(tr1, i), tr2.supporting_plane(),d); + typename K::comparison_result res_v_pl= CGAL::compare_squared_distance(vertex(tr1, i), tr2.supporting_plane(),d2); if(tr_contains_proj_p(tr2, vertex(tr1, i))) if(res_v_pl==CGAL::SMALLER) return CGAL::SMALLER; diff --git a/Distance_3/include/CGAL/global_functions_distance_3.h b/Distance_3/include/CGAL/global_functions_distance_3.h new file mode 100644 index 00000000000..9dd35ef5c55 --- /dev/null +++ b/Distance_3/include/CGAL/global_functions_distance_3.h @@ -0,0 +1,99 @@ +// Copyright (c) 2025 +// GeometryFactory (France), +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Léo Valque + +#ifndef CGAL_KERNEL_GLOBAL_FUNCTIONS_DISTANCE_3_H +#define CGAL_KERNEL_GLOBAL_FUNCTIONS_DISTANCE_3_H + +// Distance functions calling the kernel functor. + +#define CGAL_SQUARED_DISTANCE_FUNCTION(A, B) \ +template \ +inline \ +typename K::FT \ +squared_distance(const A& a, const B& b) \ +{ \ + return K().compute_squared_distance_3_object()(a, b); \ +} \ +template \ +inline \ +typename K::FT \ +squared_distance(const B& a, const A& b) \ +{ \ + return K().compute_squared_distance_3_object()(b, a); \ +} + +#define CGAL_SQUARED_DISTANCE_FUNCTION_SELF(A) \ +template \ +inline \ +typename K::FT \ +squared_distance(const A& a, const A& b) \ +{ \ + return K().compute_squared_distance_3_object()(a, b); \ +} + +#define CGAL_COMPARE_SQUARED_DISTANCE_FUNCTION(A, B) \ +template \ +inline \ +typename K::Comparison_result \ +compare_squared_distance(const A& a, \ + const B& b, \ + const typename K::FT& d2) \ +{ \ + return K().compare_squared_distance_3_object()(a, b, d2); \ +} \ +template \ +inline \ +typename K::Comparison_result \ +compare_squared_distance(const B& b, \ + const A& a, \ + const typename K::FT& d2) \ +{ \ + return K().compare_squared_distance_3_object()(b, a, d2); \ +} + +#define CGAL_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(A) \ +template \ +inline \ +typename K::Comparison_result \ +compare_squared_distance(const A& a, \ + const A& b, \ + const typename K::FT& d2) \ +{ \ + return K().compare_squared_distance_3_object()(a, b, d2); \ +} + +#define CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(A) \ +CGAL_SQUARED_DISTANCE_FUNCTION_SELF(A) \ +CGAL_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(A) + +#define CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(A, B) \ +CGAL_SQUARED_DISTANCE_FUNCTION(A, B) \ +CGAL_COMPARE_SQUARED_DISTANCE_FUNCTION(A, B) + +namespace CGAL { + +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(Point_3) +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(Segment_3) +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(Line_3) +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(Plane_3) + +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Point_3, Line_3) +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Point_3, Ray_3) +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Point_3, Segment_3) + +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Point_3, Plane_3) + +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Line_3, Plane_3) + +} //namespace CGAL + +#endif // CGAL_KERNEL_GLOBAL_FUNCTIONS_DISTANCE_3_H \ No newline at end of file diff --git a/Distance_3/include/CGAL/squared_distance_3.h b/Distance_3/include/CGAL/squared_distance_3.h index 1cf8feb0bfc..f2e51c4b51f 100644 --- a/Distance_3/include/CGAL/squared_distance_3.h +++ b/Distance_3/include/CGAL/squared_distance_3.h @@ -45,4 +45,6 @@ #include +#include + #endif // CGAL_DISTANCE_3_H diff --git a/Distance_3/test/Distance_3/test_compare_distance_3.cpp b/Distance_3/test/Distance_3/test_compare_distance_3.cpp index dfc2e1c8ae9..29380eaf32a 100644 --- a/Distance_3/test/Distance_3/test_compare_distance_3.cpp +++ b/Distance_3/test/Distance_3/test_compare_distance_3.cpp @@ -115,76 +115,112 @@ private: return true; } - template - void check_compare_squared_distance(const O1& o1, const O2& o2, const FT& d, const Comparison_result& expected_result) - { - const FT res_o1o2 = CGAL::compare_squared_distance(o1, o2, d); - const FT res_o2o1 = CGAL::compare_squared_distance(o2, o1, d); + void do_intersect_check(const P&, const P&) { } - assert(res_o1o2 == expected_result); - assert(res_o2o1 == expected_result); + template + void do_intersect_check(const P& p, const O2& o2) + { + if(!o2.is_degenerate() && CGAL::do_intersect(p, o2)) + { + assert(are_equal(CGAL::squared_distance(p, o2), FT(0))); + assert(are_equal(CGAL::squared_distance(o2, p), FT(0))); + } + } + + template + void do_intersect_check(const O1& o1, const O2& o2) + { + if(!o1.is_degenerate() && !o2.is_degenerate() && CGAL::do_intersect(o1, o2)) + { + assert(are_equal(CGAL::squared_distance(o1, o2), FT(0))); + assert(are_equal(CGAL::squared_distance(o2, o1), FT(0))); + } + } + + template + void check_compare_squared_distance(const O1& o1, const O2& o2, const FT& expected_result) + { + // const FT res_o1o2 = K().compare_squared_distance_3_object()(o1, o2, expected_result); + // const FT res_o2o1 = K().compare_squared_distance_3_object()(o2, o1, expected_result); + + const bool res_e_o1o2 = (CGAL::compare_squared_distance(o1, o2, expected_result) == CGAL::EQUAL); + const bool res_e_o2o1 = (CGAL::compare_squared_distance(o2, o1, expected_result) == CGAL::EQUAL); + const bool res_s_o1o2 = (CGAL::compare_squared_distance(o1, o2, expected_result+1) == CGAL::SMALLER); + const bool res_s_o2o1 = (CGAL::compare_squared_distance(o2, o1, expected_result+1) == CGAL::SMALLER); + const bool res_l_o1o2 = (CGAL::compare_squared_distance(o1, o2, expected_result-1) == CGAL::LARGER); + const bool res_l_o2o1 = (CGAL::compare_squared_distance(o2, o1, expected_result-1) == CGAL::LARGER); + + assert(res_e_o1o2); + assert(res_e_o2o1); + assert(res_s_o1o2); + assert(res_s_o2o1); + assert(res_l_o1o2); + assert(res_l_o2o1); + } + + template + void check_compare_squared_distance_with_bound(const O1& o1, const O2& o2, const FT& upper_bound) + { + // const FT res_o1o2 = K().compare_squared_distance_3_object()(o1, o2, expected_result); + // const FT res_o2o1 = K().compare_squared_distance_3_object()(o2, o1, expected_result); + + const bool res_o1o2 = (CGAL::compare_squared_distance(o1, o2, upper_bound) == CGAL::SMALLER); + const bool res_o2o1 = (CGAL::compare_squared_distance(o2, o1, upper_bound) == CGAL::SMALLER); + + assert(res_o1o2); + assert(res_o2o1); } private: void P_P() { std::cout << "Point - Point" << std::endl; - check_compare_squared_distance(p(0, 0, 0), p(0, 0, 0), 0, CGAL::EQUAL); - check_compare_squared_distance(p(0, 0, 0), p(0, 0, 0), 1, CGAL::SMALLER); - check_compare_squared_distance(p(3, 5, 7), p(0, 0, 0), 83, CGAL::EQUAL); - check_compare_squared_distance(p(3, 5, 7), p(0, 0, 0), 80, CGAL::LARGER); + check_compare_squared_distance(p(0, 0, 0), p(0, 0, 0), 0); + check_compare_squared_distance(p(3, 5, 7), p(0, 0, 0), 83); } void P_S() { std::cout << "Point - Segment" << std::endl; - check_compare_squared_distance(p(0, 1, 2), S{p(-3, 0, 0), p( 2, 0, 0)}, 4, CGAL::LARGER); - check_compare_squared_distance(p(0, 1, 2), S{p(-3, 0, 0), p( 2, 0, 0)}, 5, CGAL::EQUAL); - check_compare_squared_distance(p(0, 1, 2), S{p(-3, 0, 0), p( 2, 0, 0)}, 6, CGAL::SMALLER); - - check_compare_squared_distance(p(0, 1, 2), S{p( 3, 0, 0), p( 2, 0, 0)}, 8, CGAL::LARGER); - check_compare_squared_distance(p(0, 1, 2), S{p( 3, 0, 0), p( 2, 0, 0)}, 9, CGAL::EQUAL); - check_compare_squared_distance(p(0, 1, 2), S{p( 3, 0, 0), p( 2, 0, 0)}, 10, CGAL::SMALLER); - check_compare_squared_distance(p(0, 1, 2), S{p( 2, 0, 0), p( 3, 0, 0)}, 10, CGAL::SMALLER); - - check_compare_squared_distance(p(6, 1, 2), S{p( 2, 0, 0), p( 3, 0, 0)}, 13, CGAL::LARGER); - check_compare_squared_distance(p(6, 1, 2), S{p( 2, 0, 0), p( 3, 0, 0)}, 14, CGAL::EQUAL); - check_compare_squared_distance(p(6, 1, 2), S{p( 2, 0, 0), p( 3, 0, 0)}, 15, CGAL::SMALLER); + check_compare_squared_distance(p(0, 1, 2), S{p(-3, 0, 0), p( 2, 0, 0)}, 5); + check_compare_squared_distance(p(0, 1, 2), S{p( 3, 0, 0), p( 2, 0, 0)}, 9); + check_compare_squared_distance(p(0, 1, 2), S{p( 2, 0, 0), p( 3, 0, 0)}, 9); + check_compare_squared_distance(p(6, 1, 2), S{p( 2, 0, 0), p( 3, 0, 0)}, 14); } // void P_T() // { // std::cout << "Point - Triangle" << std::endl; -// check_squared_distance(p(0, 1, 2), T{p(0, 0, 0), p(2, 0, 0), p(0, 2, 0)}, 4); +// check_compare_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); +// check_compare_squared_distance (p(-1, -1, 0), t, 2); +// check_compare_squared_distance (p(-1, 0, 0), t, 1); +// check_compare_squared_distance (p(0, 0, 0), t, 0); +// check_compare_squared_distance (p(1, 0, 0), t, 0); +// check_compare_squared_distance (p(4, 0, 0), t, 1); +// check_compare_squared_distance (p(1, -1, 0), t, 1); +// check_compare_squared_distance (p(1, 1, 1), t, 1); +// check_compare_squared_distance (p(1, 0, 1), t, 1); +// check_compare_squared_distance (p(0, 0, 1), t, 1); // // Degenerate -// check_squared_distance (p(1, 2, 3), T(p(4,3,2), p(4,3,2), p(4,3,2)), squared_distance(p(1, 2, 3), p(4,3,2))); -// check_squared_distance (p(1, 2, 3), T(p(4,3,2), p(10,12,3), p(4,3,2)), squared_distance(p(1, 2, 3), p(4,3,2))); -// check_squared_distance (p(0, 0, 0), T(p(4,3,2), p(4,-3,-2), p(4,3,2)), squared_distance(p(0, 0, 0), p(4,0,0))); +// check_compare_squared_distance (p(1, 2, 3), T(p(4,3,2), p(4,3,2), p(4,3,2)), squared_distance(p(1, 2, 3), p(4,3,2))); +// check_compare_squared_distance (p(1, 2, 3), T(p(4,3,2), p(10,12,3), p(4,3,2)), squared_distance(p(1, 2, 3), p(4,3,2))); +// check_compare_squared_distance (p(0, 0, 0), T(p(4,3,2), p(4,-3,-2), p(4,3,2)), squared_distance(p(0, 0, 0), p(4,0,0))); // // On the triangle -// check_squared_distance (p(7, 1, -5), T(p(2,9,8), p(-4,-3,-5), p(7, 1, -5)), 0); // vertex -// check_squared_distance (p(7, 1, -5), T(p(14,2,-10), p(-7,-1,5), p(8, 3, -1)), 0); // edge -// check_squared_distance (p(1, 4, -3), T(p(0,-8,-3), p(-5,14,-3), p(10, 1, -3)), 0); // face +// check_compare_squared_distance (p(7, 1, -5), T(p(2,9,8), p(-4,-3,-5), p(7, 1, -5)), 0); // vertex +// check_compare_squared_distance (p(7, 1, -5), T(p(14,2,-10), p(-7,-1,5), p(8, 3, -1)), 0); // edge +// check_compare_squared_distance (p(1, 4, -3), T(p(0,-8,-3), p(-5,14,-3), p(10, 1, -3)), 0); // face // // General -// check_squared_distance (p(-15, 1, 0), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), 25); -// check_squared_distance (p(-5, 0, 0), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), squared_distance(p(-5, 0, 0), S(p(-10, 1, 0), p(0,0,0)))); -// check_squared_distance (p(0, -3, 0), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), 9); -// check_squared_distance (p(3, -3, 0), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), squared_distance(p(3, -3, 0), S(p(0,0,0), p(10,0,0)))); -// check_squared_distance (p(16, 1, 1), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), 38); -// check_squared_distance (p(5, 5, 2), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), squared_distance(p(5, 5, 2), S(p(10,0,0), p(-10,1,0)))); +// check_compare_squared_distance (p(-15, 1, 0), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), 25); +// check_compare_squared_distance (p(-5, 0, 0), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), squared_distance(p(-5, 0, 0), S(p(-10, 1, 0), p(0,0,0)))); +// check_compare_squared_distance (p(0, -3, 0), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), 9); +// check_compare_squared_distance (p(3, -3, 0), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), squared_distance(p(3, -3, 0), S(p(0,0,0), p(10,0,0)))); +// check_compare_squared_distance (p(16, 1, 1), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), 38); +// check_compare_squared_distance (p(5, 5, 2), T(p(-10, 1, 0), p(0,0,0), p(10,0,0)), squared_distance(p(5, 5, 2), S(p(10,0,0), p(-10,1,0)))); // for(int i=0; i1) in the code -// for(int j=-2;j<4; j+=2) -// { -// for(int k=-3;k<3; k+=2) -// { -// P p2{j, k, z}; -// P p3{j, k+2, z}; + // translations of (0, -1, z) -- (0, 1, z) to have all variations of x&y (<0, [0;1]; >1) in the code + for(int j=-2;j<4; j+=2) + { + for(int k=-3;k<3; k+=2) + { + P p2{j, k, z}; + P p3{j, k+2, z}; -// // to explicit the expected distances -// if(j == -2 && k == -3) -// check_squared_distance(S{p0, p1}, S{p2, p3}, CGAL::squared_distance(p3, p0)); -// else if(j == -2 && k == -1) -// check_squared_distance(S{p0, p1}, S{p2, p3}, 1); -// else if(j == -2 && k == 1) -// check_squared_distance(S{p0, p1}, S{p2, p3}, CGAL::squared_distance(p2, p0)); -// else if(j == 0 && k == -3) -// check_squared_distance(S{p0, p1}, S{p2, p3}, 1); -// else if(j == 0 && k == -1) -// check_squared_distance(S{p0, p1}, S{p2, p3}, 0); -// else if(j == 0 && k == 1) -// check_squared_distance(S{p0, p1}, S{p2, p3}, 1); -// else if(j == 2 && k == -3) -// check_squared_distance(S{p0, p1}, S{p2, p3}, CGAL::squared_distance(p3, p1)); -// else if(j == 2 && k == -1) -// check_squared_distance(S{p0, p1}, S{p2, p3}, 1); -// else if(j == 2 && k == 1) -// check_squared_distance(S{p0, p1}, S{p2, p3}, CGAL::squared_distance(p2, p1)); -// } -// } + // to explicit the expected distances + if(j == -2 && k == -3) + check_compare_squared_distance(S{p0, p1}, S{p2, p3}, CGAL::squared_distance(p3, p0)); + else if(j == -2 && k == -1) + check_compare_squared_distance(S{p0, p1}, S{p2, p3}, 1); + else if(j == -2 && k == 1) + check_compare_squared_distance(S{p0, p1}, S{p2, p3}, CGAL::squared_distance(p2, p0)); + else if(j == 0 && k == -3) + check_compare_squared_distance(S{p0, p1}, S{p2, p3}, 1); + else if(j == 0 && k == -1) + check_compare_squared_distance(S{p0, p1}, S{p2, p3}, 0); + else if(j == 0 && k == 1) + check_compare_squared_distance(S{p0, p1}, S{p2, p3}, 1); + else if(j == 2 && k == -3) + check_compare_squared_distance(S{p0, p1}, S{p2, p3}, CGAL::squared_distance(p3, p1)); + else if(j == 2 && k == -1) + check_compare_squared_distance(S{p0, p1}, S{p2, p3}, 1); + else if(j == 2 && k == 1) + check_compare_squared_distance(S{p0, p1}, S{p2, p3}, CGAL::squared_distance(p2, p1)); + } + } -// for(int i=0; i, gte::Segment<3, FT> > GTE_SS_checker; -// gte::Segment<3, FT> gte_s1{{p8.x(), p8.y(), p8.z()}, {p9.x(), p9.y(), p9.z()}}; -// gte::Segment<3, FT> gte_s2{{p3.x(), p3.y(), p3.z()}, {p2.x(), p2.y(), p2.z()}}; -// auto gte_res = GTE_SS_checker(gte_s1, gte_s2); -// std::cout << "-------" << std::endl; -// std::cout << "old: " << CGAL::internal::squared_distance_old(s89, s32, K()) << std::endl; -// std::cout << "dist (GTE) : " << gte_res.sqrDistance << std::endl; -// #endif +#ifdef CGAL_USE_GTE_AS_SANITY_CHECK + gte::DCPQuery, gte::Segment<3, FT> > GTE_SS_checker; + gte::Segment<3, FT> gte_s1{{p8.x(), p8.y(), p8.z()}, {p9.x(), p9.y(), p9.z()}}; + gte::Segment<3, FT> gte_s2{{p3.x(), p3.y(), p3.z()}, {p2.x(), p2.y(), p2.z()}}; + auto gte_res = GTE_SS_checker(gte_s1, gte_s2); + std::cout << "-------" << std::endl; + std::cout << "old: " << CGAL::internal::squared_distance_old(s89, s32, K()) << std::endl; + std::cout << "dist (GTE) : " << gte_res.sqrDistance << std::endl; +#endif -// // Because do_intersect() with constructions is likely to return 'false' even for overlaps -// assert(are_equal(CGAL::squared_distance(s89, s32), result, false /*verbose*/) || -// are_equal(CGAL::squared_distance(s32, s89), FT(0))); -// } + // Because do_intersect() with constructions is likely to return 'false' even for overlaps + assert(are_equal(CGAL::squared_distance(s89, s32), result, false /*verbose*/) || + are_equal(CGAL::squared_distance(s32, s89), FT(0))); + } -// // completely generic -// S s1{p0, p1}, s2{p2, p3}; -// do_intersect_check(s1, s2); + // completely generic + S s1{p0, p1}, s2{p2, p3}; + do_intersect_check(s1, s2); -// #ifdef CGAL_USE_GTE_AS_SANITY_CHECK -// gte::DCPQuery, gte::Segment<3, FT> > GTE_SS_checker; -// gte::Segment<3, FT> gte_s1{{p0.x(), p0.y(), p0.z()}, {p1.x(), p1.y(), p1.z()}}; -// gte::Segment<3, FT> gte_s2{{p2.x(), p2.y(), p2.z()}, {p3.x(), p3.y(), p3.z()}}; -// auto gte_res = GTE_SS_checker(gte_s1, gte_s2); +#ifdef CGAL_USE_GTE_AS_SANITY_CHECK + gte::DCPQuery, gte::Segment<3, FT> > GTE_SS_checker; + gte::Segment<3, FT> gte_s1{{p0.x(), p0.y(), p0.z()}, {p1.x(), p1.y(), p1.z()}}; + gte::Segment<3, FT> gte_s2{{p2.x(), p2.y(), p2.z()}, {p3.x(), p3.y(), p3.z()}}; + auto gte_res = GTE_SS_checker(gte_s1, gte_s2); -// std::cout << "dist (CGAL) : " << CGAL::squared_distance(s1, s2) << std::endl; -// std::cout << "dist (GTE) : " << gte_res.sqrDistance << std::endl; -// assert(are_equal(CGAL::squared_distance(s1, s2), gte_res.sqrDistance)); -// #endif -// } + std::cout << "dist (CGAL) : " << CGAL::squared_distance(s1, s2) << std::endl; + std::cout << "dist (GTE) : " << gte_res.sqrDistance << std::endl; + assert(are_equal(CGAL::squared_distance(s1, s2), gte_res.sqrDistance)); +#endif + } -// // a few brute force checks: sample a segments and use squared_distance(P3, S3) -// for(int i=0; i<10; ++i) -// { -// P p0 = random_point(); -// P p1 = random_point(); -// P p2 = random_point(); -// P p3 = random_point(); + // a few brute force checks: sample a segments and use squared_distance(P3, S3) + for(int i=0; i<10; ++i) + { + P p0 = random_point(); + P p1 = random_point(); + P p2 = random_point(); + P p3 = random_point(); -// S s01{p0, p1}; -// S s23{p2, p3}; + S s01{p0, p1}; + S s23{p2, p3}; -// FT upper_bound = CGAL::squared_distance(p0, p2); + FT upper_bound = CGAL::squared_distance(p0, p2); -// V p01 = V{p0, p1} / FT(N); -// for(int l=0; l >(r, 1e-5).run(); // Test > >(r, 1e-5).run(); - Test >(r, 0).run(); + // Test >(r, 0).run(); const double epick_eps = 10 * std::numeric_limits::epsilon(); Test(r, epick_eps).run(); diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index f88a4ea4c2b..7c7de293777 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -945,8 +945,10 @@ namespace CommonKernelFunctors { Needs_FT operator()(const T1& p, const T2& q, const T3& r, const T4& s) const { - return CGAL::compare(internal::squared_distance(p, q, K()), - internal::squared_distance(r, s, K())); + // return internal::compare_squared_distance(p, q, K(), internal::squared_distance(p, q, K())); + return internal::compare_squared_distance(p, q, K(), internal::squared_distance(r, s, K())); + // return CGAL::compare(internal::squared_distance(p, q, K()), + // internal::squared_distance(r, s, K())); } }; diff --git a/Kernel_23/include/CGAL/Kernel/global_functions_3.h b/Kernel_23/include/CGAL/Kernel/global_functions_3.h index 6fb9b975647..b660dd84f05 100644 --- a/Kernel_23/include/CGAL/Kernel/global_functions_3.h +++ b/Kernel_23/include/CGAL/Kernel/global_functions_3.h @@ -412,16 +412,6 @@ compare_slope(const Point_3 &p, return internal::compare_slope(p, q, r, s, K()); } -template < class K > -inline -typename K::Comparison_result -compare_squared_distance(const Point_3 &p, - const Point_3 &q, - const typename K::FT &d2) -{ - return internal::compare_squared_distance(p, q, d2, K()); -} - template < class K > inline typename K::Comparison_result