From 8bdc4b4f5d25035083e98ef7fcc9c52387869a97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o=20Valque?= Date: Wed, 19 Feb 2025 11:55:57 +0100 Subject: [PATCH] add test_compare_distance_3.cpp --- .../CGAL/Distance_3/Point_3_Triangle_3.h | 2 +- .../Distance_3/test_compare_distance_3.cpp | 698 ++++++++++++++++++ 2 files changed, 699 insertions(+), 1 deletion(-) create mode 100644 Distance_3/test/Distance_3/test_compare_distance_3.cpp 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 18de39c4024..02002ff2600 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 @@ -293,7 +293,7 @@ squared_distance(const typename K::Triangle_3& t, } template -inline typename K::Comparison_result +typename K::Comparison_result compare_squared_distance_to_triangle(const typename K::Point_3& pt, const typename K::Point_3& t0, const typename K::Point_3& t1, diff --git a/Distance_3/test/Distance_3/test_compare_distance_3.cpp b/Distance_3/test/Distance_3/test_compare_distance_3.cpp new file mode 100644 index 00000000000..dfc2e1c8ae9 --- /dev/null +++ b/Distance_3/test/Distance_3/test_compare_distance_3.cpp @@ -0,0 +1,698 @@ +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +// #define CGAL_USE_GTE_AS_SANITY_CHECK +#ifdef CGAL_USE_GTE_AS_SANITY_CHECK +#include +#include +#endif + +#include +#include + +struct randomint +{ + randomint() ; + int get() const { return sequence[cur]; } + int next() { + cur = (cur + 1) % 11; + return get(); + } + +private: + int sequence[11]; + int cur; +}; + +inline randomint::randomint() +{ + cur = 0; + sequence[0] = 19; + sequence[1] = 5; + sequence[2] = 17; + sequence[3] = 13; + sequence[4] = 29; + sequence[5] = 2; + sequence[6] = 23; + sequence[7] = 31; + sequence[8] = 3; + sequence[9] = 37; + sequence[10] = 11; +} + +randomint ri; + +template +struct Test +{ + typedef typename K::RT RT; + typedef typename K::FT FT; + typedef typename K::Comparison_result Comparison_result; + typedef typename K::Point_3 P; + typedef typename K::Segment_3 S; + typedef typename K::Vector_3 V; + typedef typename K::Ray_3 R; + typedef typename K::Line_3 L; + typedef typename K::Triangle_3 T; + typedef typename K::Plane_3 Pl; + typedef typename K::Tetrahedron_3 Tet; + typedef typename K::Iso_cuboid_3 Cub; + +private: + CGAL::Random& r; + const double epsilon = 1e-14; + int N = 10; + double m = 0, M = 1; + +public: + Test(CGAL::Random& r, const double epsilon) : r(r), epsilon(epsilon) { } + +private: + inline RT to_nt(int d) const { return RT(d); } + + P p(int x, int y, int z) + { + int w = ri.next(); + return P(to_nt(x*w), to_nt(y*w), to_nt(z*w), to_nt(w)); + } + + P random_point() const + { + return P(FT(r.get_double(m, M)), FT(r.get_double(m, M)), FT(r.get_double(m, M))); + } + + Pl pl(int a, int b, int c, int d) + { + int w = ri.next(); + return Pl(to_nt(a*w), to_nt(b*w), to_nt(c*w), to_nt(d*w)); + } + +private: + template + bool are_equal(const Type& t1, const Type& t2, const bool verbose = true) + { + const FT diff = CGAL::abs(t1 - t2); + if(diff > std::numeric_limits::epsilon() && + diff > epsilon * (CGAL::abs(t1) + CGAL::abs(t2))) + { + if(verbose) + { + std::cerr << "Approximate comparison failed (t1|t2): got " << t1 << " but expected " << t2 << std::endl; + std::cerr << "Diff: " << CGAL::abs(t1 - t2) << " vs eps: " << epsilon * (CGAL::abs(t1) + CGAL::abs(t2)) << std::endl; + } + return false; + } + + 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); + + assert(res_o1o2 == expected_result); + assert(res_o2o1 == expected_result); + } + +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); + } + + 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); + } + +// 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); + +// 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); + +// // 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))); + +// // 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 + +// // 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)))); + +// 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}; + +// // 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)); +// } +// } + +// 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 + +// // 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); + +// #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 +// } + +// // 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}; + +// FT upper_bound = CGAL::squared_distance(p0, p2); + +// V p01 = V{p0, p1} / FT(N); +// for(int l=0; l, gte::Triangle3 > GTE_TT_checker; +// gte::Triangle3 gte_tr1{{p0.x(), p0.y(), p0.z()}, {p1.x(), p1.y(), p1.z()}, {p2.x(), p2.y(), p2.z()}}; +// gte::Triangle3 gte_tr2{{p3.x(), p3.y(), p3.z()}, {p4.x(), p4.y(), p4.z()}, {p5.x(), p5.y(), p5.z()}}; +// auto gte_res = GTE_TT_checker(gte_tr1, gte_tr2); + +// std::cout << "dist (CGAL) : " << CGAL::squared_distance(tr1, tr2) << std::endl; +// std::cout << "dist (GTE) : " << gte_res.sqrDistance << std::endl; +// std::cout << "diff CGAL GTE : " << (gte_res.sqrDistance - CGAL::squared_distance(tr1, tr2)) << std::endl; + +// // don't assert on purpose, GTE has slightly (10^-30 different results, even with an exact NT) +// are_equal(CGAL::squared_distance(tr1, tr2), gte_res.sqrDistance); +// #endif +// } +// } + +public: + void run() + { + std::cout << "Kernel: " << typeid(K).name() << std::endl; + + P_P(); + P_S(); + P_R(); + P_L(); + // P_T(); + P_Pl(); + // P_Tet(); + + // S_S(); + // S_R(); + // S_L(); + // S_Pl(); + + // R_R(); + // R_L(); + // R_Pl(); + + L_L(); + // L_Pl(); + + // T_T(); + + // Pl_Pl(); + } +}; + +int main() +{ + std::cout.precision(17); + std::cerr.precision(17); + + std::cout << "3D Distance tests" << std::endl; + + CGAL::Random r; + std::cout << "random seed = " << r.get_seed() << std::endl; + + // @todo Some tests are too difficult for these kernels +// Test >(r, 1e-5).run(); +// Test >(r, 1e-5).run(); +// Test > >(r, 1e-5).run(); + + Test >(r, 0).run(); + + const double epick_eps = 10 * std::numeric_limits::epsilon(); + Test(r, epick_eps).run(); + + Test(r, 0).run(); + + std::cout << "Done!" << std::endl; + + return EXIT_SUCCESS; +}