diff --git a/Distance_3/benchmark/Distance_3/bench_thingi10k_3.cpp b/Distance_3/benchmark/Distance_3/bench_thingi10k_3.cpp index 70fe503d513..f2b17155f89 100644 --- a/Distance_3/benchmark/Distance_3/bench_thingi10k_3.cpp +++ b/Distance_3/benchmark/Distance_3/bench_thingi10k_3.cpp @@ -107,7 +107,7 @@ int main(int argc, char** argv) std::cout << "3D Distance bench" << std::endl; - std::vector input_points; + std::vector::Point_3> input_points; std::vector> input_triangles; if (!CGAL::IO::read_polygon_soup(filename, input_points, input_triangles)) @@ -126,13 +126,24 @@ int main(int argc, char** argv) // Test().run(filename, max*max); // Test().run(filename, 10*max*max/input_points.size()); - Test().run(filename, 100*max/input_points.size()); - Test().run(filename, max/input_points.size()); - Test().run(filename, max/(100*input_points.size())); - Test().run(filename, 100*max/input_points.size()); - Test().run(filename, max/input_points.size()); - Test().run(filename, max/(100*input_points.size())); + double average_length=0; + double min_sq_length=CGAL::squared_distance(input_points[input_triangles[0][0]],input_points[input_triangles[0][1]]); + for(auto &tr: input_triangles){ + for(int i=0; i<3; ++i){ + double l=CGAL::squared_distance(input_points[tr[i]], input_points[tr[(i+1)%3]]); + min_sq_length=(std::min)(min_sq_length, l); + average_length+=std::sqrt(l); + } + } + average_length/=(3*input_triangles.size()); + + // Test().run(filename, 100*max/input_points.size()); + Test().run(filename, average_length*average_length/256); + Test().run(filename, min_sq_length*4); + + Test().run(filename, average_length*average_length/256); + Test().run(filename, min_sq_length*4); std::cout << "Done!" << std::endl; 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 4ddb7ad1702..6331d938d82 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 @@ -46,7 +46,7 @@ compare_squared_distance(const typename K::Plane_3& plane1, const K& k, const typename K::FT& d2) { - return compare(squared_distance(plane1,plane2,k), d2); + return compare(squared_distance(plane1,plane2,k), d2); } } // namespace internal 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 cb2bdd7401d..bb4195cbb75 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 @@ -94,22 +94,7 @@ compare_squared_distance(const typename K::Point_3& pt, const K& k, const typename K::FT &d2) { - typedef typename K::Vector_3 Vector_3; - - 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); - - //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; - - if(!is_acute_angle(dir, diff, k)) - return compare(diff*diff, d2); - - return res_pl; + return compare(squared_distance(pt, ray, k), d2); } template 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 f13fa5ae5de..07f79bfd912 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 @@ -116,6 +116,7 @@ compare_squared_distance(const typename K::Point_3& pt, const K& k, const typename K::FT& d2) { + //Doing an early exit was slower. return compare(squared_distance(pt, seg, k), d2); } diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Tetrahedron_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Tetrahedron_3.h index e4505ec0f30..031a60a6650 100644 --- a/Distance_3/include/CGAL/Distance_3/Point_3_Tetrahedron_3.h +++ b/Distance_3/include/CGAL/Distance_3/Point_3_Tetrahedron_3.h @@ -155,6 +155,10 @@ compare_squared_distance(const typename K::Point_3& pt, typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object(); typename K::Orientation_3 orientation = k.orientation_3_object(); + /* The content of this function is very similar with the one above, the difference is we can exit earlier if + we found a triangle closer than d or plane farther than d since we do not need the exact distance. + (there are also early exits in calling functions) */ + bool on_bounded_side = true; bool inside_or_far_to_the_plane = false; const Point_3& t0 = vertex(tet, 0); 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 759062de81f..0f5d86f7a06 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 @@ -309,6 +309,9 @@ compare_squared_distance_to_triangle(const typename K::Point_3& pt, 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(); + /* The content of this function is very similar with the one above, the difference is we can exit earlier if + we found a segment closer than d or if the point is farther than d to the plane since we do not need the exact distance */ + const Vector_3 e1 = vector(t0, t1); const Vector_3 oe3 = vector(t0, t2); const Vector_3 normal = wcross(e1, oe3, k); diff --git a/Distance_3/include/CGAL/Distance_3/Segment_3_Line_3.h b/Distance_3/include/CGAL/Distance_3/Segment_3_Line_3.h index 6fa52ef1f4e..f5d3336b5d6 100644 --- a/Distance_3/include/CGAL/Distance_3/Segment_3_Line_3.h +++ b/Distance_3/include/CGAL/Distance_3/Segment_3_Line_3.h @@ -62,27 +62,8 @@ squared_distance(const typename K::Segment_3& seg, const RT sdm_ss2l = wdot(perpend2line, start_min_lp, k); const RT sdm_se2l = wdot(perpend2line, end_min_lp, k); - // if(sdm_ss2l < RT(0)) { - // crossing = (sdm_se2l >= RT(0)); - // } else { - // if(sdm_se2l <= RT(0)) { - // crossing = true; - // } else { - // crossing = (sdm_ss2l == RT(0)); - // } - // } crossing = (sdm_ss2l*sdm_se2l) <= RT(0); - // if(crossing) { - // return squared_distance_to_plane(normal, start_min_lp, k); - // } else { - // const RT dm = distance_measure_sub(sdm_ss2l, sdm_se2l, start_min_lp, end_min_lp, k); - // if(dm <= RT(0)) { - // return squared_distance_to_line(linedir, start_min_lp, k); - // } else { - // return squared_distance_to_line(linedir, end_min_lp, k); - // } - // } if(crossing) return squared_distance_to_plane(normal, start_min_lp, k); else @@ -105,55 +86,7 @@ compare_squared_distance(const typename K::Segment_3& seg, const K& k, const typename K::FT& d2) { - // typedef typename K::RT RT; - // typedef typename K::Point_3 Point_3; - // 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(); - - // const Point_3& linepoint = line.point(); - // const Point_3& start = seg.source(); - // const Point_3& end = seg.target(); - - // if(start == end) - // return csq_dist(start, line, d2); - - // const Vector_3 linedir = line.direction().vector(); - // const Vector_3 segdir = seg.direction().vector(); - // const Vector_3 normal = wcross(segdir, linedir, k); - - // if(is_null(normal, k)) - // return compare(squared_distance_to_line(linedir, vector(linepoint,start), k), d2); - - // bool crossing; - - // const Vector_3 perpend2line = wcross(linedir, normal, k); - // const Vector_3 start_min_lp = vector(linepoint, start); - // const Vector_3 end_min_lp = vector(linepoint, end); - // const RT sdm_ss2l = wdot(perpend2line, start_min_lp, k); - // const RT sdm_se2l = wdot(perpend2line, end_min_lp, k); - - // if(sdm_ss2l < RT(0)) { - // crossing = (sdm_se2l >= RT(0)); - // } else { - // if(sdm_se2l <= RT(0)) { - // crossing = true; - // } else { - // crossing = (sdm_ss2l == RT(0)); - // } - // } - - // if(crossing) { - // return squared_distance_to_plane(normal, start_min_lp, k); - // } else { - // const RT dm = distance_measure_sub(sdm_ss2l, sdm_se2l, start_min_lp, end_min_lp, k); - // if(dm <= RT(0)) { - // return squared_distance_to_line(linedir, start_min_lp, k); - // } else { - // return squared_distance_to_line(linedir, end_min_lp, k); - // } - // } + // Perform an early exit was slower return compare(squared_distance(seg, line, k), d2); } 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 310fd80f1e2..2808ad85fa7 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 @@ -206,6 +206,10 @@ compare_squared_distance(const typename K::Segment_3& s1, 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(); + /* The content of this function is very similar with the one above, the difference is we can exit earlier if + the supporting line are farther than d since we do not need the exact distance. */ + + const Point_3& p1 = vertex(s1, 0); const Point_3& q1 = vertex(s1, 1); const Point_3& p2 = vertex(s2, 0); 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 119af808e72..6584e669f34 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 @@ -226,10 +226,14 @@ compare_squared_distance_disjoint(const typename K::Triangle_3& tr1, typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object(); typename K::Compare_squared_distance_3 csq_dist = k.compare_squared_distance_3_object(); - typename K::Comparison_result res(LARGER); + typename K::Comparison_result res(LARGER); + + // The tiangle are supposed to be disjoint + assertion(!do_intersect(tr1, tr2)); for(int i=0; i<3; ++i) { + //Compare the distance between edges for(int j=0; j<3; ++j) { typename K::Comparison_result temp_res_ss=csq_dist(Segment_3(vertex(tr1, i%3), vertex(tr1, (i+1)%3)),Segment_3(vertex(tr2, j%3), vertex(tr2, (j+1)%3)),d2); @@ -238,14 +242,15 @@ compare_squared_distance_disjoint(const typename K::Triangle_3& tr1, res=smaller_of(res, temp_res_ss); } + //Compare the distance between vertices and triangles typename K::Comparison_result temp_res_v_pl= csq_dist(vertex(tr1, i), tr2,d2); if(temp_res_v_pl==SMALLER) - return SMALLER; + return SMALLER; res=smaller_of(res, temp_res_v_pl); temp_res_v_pl= csq_dist(vertex(tr2, i), tr1,d2); if(temp_res_v_pl==SMALLER) - return SMALLER; + return SMALLER; res=smaller_of(res, temp_res_v_pl); } return res; @@ -258,6 +263,7 @@ compare_squared_distance(const typename K::Triangle_3& tr1, const typename K::Triangle_3& tr2, const K& k, const typename K::FT& d2){ + //TODO did something more intelligent (sq_dist and csq_dist does not exist for Segment-Triangle) if(tr1.is_degenerate() || tr2.is_degenerate()) return compare(squared_distance(tr1,tr2, k), d2); if(do_intersect(tr1, tr2)) 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 c7787758a12..1c367d7fc8b 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 @@ -97,8 +97,10 @@ wdot(const typename K::Point_3 &p, return wdot_tag(p, q, r, k, tag); } + static double diff_of_products(const double a, const double b, const double c, const double d) { +// Kahan method, less numerical error #if 1 double w = d * c; double e = std::fma(c, -d, w); @@ -122,10 +124,6 @@ wcross(const typename K::Vector_3 &u, const K&) { typedef typename K::Vector_3 Vector_3; - // return Vector_3( - // u.hy()*v.hz() - u.hz()*v.hy(), - // u.hz()*v.hx() - u.hx()*v.hz(), - // u.hx()*v.hy() - u.hy()*v.hx()); return Vector_3( diff_of_products(u.hy(),v.hz(),u.hz(),v.hy()), diff_of_products(u.hz(),v.hx(),u.hx(),v.hz()), @@ -229,35 +227,6 @@ squared_distance_to_plane(const typename K::Vector_3& normal, return Rational_traits().make_rational(num, den); } -template -void -signed_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 = dot<0?-square(dot):square(dot); - den = wmult((K*)0, squared_length, diff.hw(), diff.hw()); -} - -template -typename K::FT -signed_squared_distance_to_plane(const typename K::Vector_3& normal, - const typename K::Vector_3& diff, - const K& k) -{ - typedef typename K::RT RT; - typedef typename K::FT FT; - RT num, den; - signed_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, diff --git a/Distance_3/include/CGAL/global_functions_distance_3.h b/Distance_3/include/CGAL/global_functions_distance_3.h index b4868ddc930..23775ec8535 100644 --- a/Distance_3/include/CGAL/global_functions_distance_3.h +++ b/Distance_3/include/CGAL/global_functions_distance_3.h @@ -82,10 +82,11 @@ 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(Plane_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(Ray_3) +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(Segment_3) CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(Triangle_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) @@ -96,18 +97,12 @@ CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Point_3, Tetrahedron_3) CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Line_3, Plane_3) -//Not yet modify +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Ray_3, Line_3) +CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Ray_3, Plane_3) CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Segment_3, Line_3) CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Segment_3, Ray_3) CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Segment_3, Plane_3) -// CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Segment_3, Triangle_3) - -CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Ray_3, Line_3) -CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION(Ray_3, Plane_3) -CGAL_COMPUTE_AND_COMPARE_SQUARED_DISTANCE_FUNCTION_SELF(Ray_3) - - } //namespace CGAL diff --git a/Distance_3/test/Distance_3/test_distance_3.cpp b/Distance_3/test/Distance_3/test_distance_3.cpp index 4686a7cae8f..5612ef37772 100644 --- a/Distance_3/test/Distance_3/test_distance_3.cpp +++ b/Distance_3/test/Distance_3/test_distance_3.cpp @@ -266,9 +266,6 @@ private: P p3 = random_point(); P q = random_point(); - // Failed on seed (1740146243) with an inexact kernel without the bound increasing - // The failed is due of the minimum realized by an edge and the computation is slightly different - // between the distance compute between the two triangles check_squared_distance_with_bound(q, Tet(p0, p1, p2, p3), squared_distance(q, T(p0, p1, p2))*(1 + 1e-10)); check_squared_distance_with_bound(q, Tet(p0, p1, p2, p3), squared_distance(q, T(p1, p3, p2))*(1 + 1e-10)); check_squared_distance_with_bound(q, Tet(p0, p1, p2, p3), squared_distance(q, T(p0, p3, p2))*(1 + 1e-10)); diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 7c7de293777..c932ea5acfe 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -937,7 +937,6 @@ namespace CommonKernelFunctors { Needs_FT operator()(const T1& p, const T2& q, const FT& d2) const { - // return CGAL::compare(internal::squared_distance(p, q, K()), d2); return internal::compare_squared_distance(p, q, K(), d2); } @@ -945,10 +944,7 @@ namespace CommonKernelFunctors { Needs_FT operator()(const T1& p, const T2& q, const T3& r, const T4& s) const { - // 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())); } };