diff --git a/Distance_3/test/Distance_3/test_distance_3.cpp b/Distance_3/test/Distance_3/test_distance_3.cpp index 31265cfbb01..3d219582b6e 100644 --- a/Distance_3/test/Distance_3/test_distance_3.cpp +++ b/Distance_3/test/Distance_3/test_distance_3.cpp @@ -8,18 +8,15 @@ #include #include -#include - +// #define CGAL_USE_GTE_AS_SANITY_CHECK +#ifdef CGAL_USE_GTE_AS_SANITY_CHECK #include #include - -#include +#endif #include #include -const double epsilon = 1e-14; - struct randomint { randomint() ; @@ -64,16 +61,17 @@ struct Test 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 bool exact; + const double epsilon = 1e-14; int N = 1000; double m = 0, M = 1; public: - Test(CGAL::Random& r, bool exact) : r(r), exact(exact) { } + Test(CGAL::Random& r, const double epsilon) : r(r), epsilon(epsilon) { } private: inline RT to_nt(int d) const { return RT(d); } @@ -99,31 +97,16 @@ private: template bool are_equal(const Type& t1, const Type& t2, const bool verbose = true) { - if(exact) + const FT diff = CGAL::abs(t1 - t2); + if(diff > std::numeric_limits::epsilon() && + diff > epsilon * (CGAL::abs(t1) + CGAL::abs(t2))) { - if(t1 != t2) + if(verbose) { - if(verbose) - { - std::cerr << "Approximate comparison failed: got " << t1 << " but expected " << t2 << std::endl; - std::cerr << "Diff: " << CGAL::abs(t1 - t2) << std::endl; - } - return false; - } - } - else // !exact - { - 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; + 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; @@ -143,6 +126,18 @@ private: std::cout << "result (new) = " << asd << std::endl; } + void do_intersect_check(const P&, const P&) { } + + 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) { @@ -153,14 +148,6 @@ private: } } - template - void do_intersect_check(const O1& o1, const O2& o2, const FT res) - { - // @todo remove is_degenerate() checks when do_intersect can handle degenerate objects - if(!o1.is_degenerate() && !o2.is_degenerate() && CGAL::do_intersect(o1, o2)) - assert(are_equal(res, FT(0))); - } - template void check_squared_distance(const O1& o1, const O2& o2, const FT& expected_result) { @@ -170,20 +157,18 @@ private: assert(are_equal(res_o1o2, expected_result)); assert(are_equal(res_o2o1, expected_result)); - do_intersect_check(o1, o2, res_o1o2); - do_intersect_check(o1, o2, res_o2o1); + do_intersect_check(o1, o2); + do_intersect_check(o1, o2); } template void check_squared_distance_with_bound(const O1& o1, const O2& o2, const FT& ubound) { - std::cout << "ex: " << CGAL::squared_distance(o1, o2) << " vs bound " << ubound << std::endl; - const FT res_o1o2 = CGAL::squared_distance(o1, o2); const FT res_o2o1 = CGAL::squared_distance(o2, o1); - do_intersect_check(o1, o2, res_o1o2); - do_intersect_check(o1, o2, res_o2o1); + do_intersect_check(o1, o2); + do_intersect_check(o1, o2); assert(res_o1o2 <= ubound); assert(res_o2o1 <= ubound); @@ -310,7 +295,7 @@ private: const double lambda_7 = r.get_double(1, 2); const P p7 = p3 + lambda_7 * (p3 - p2); - // [p2;p3] left of [p6;p7] + // [p2;p3] disjoint && left of [p6;p7] check_squared_distance(S{p2, p3}, S{p6, p7}, CGAL::squared_distance(p3, p6)); check_squared_distance(S{p2, p3}, S{p7, p6}, CGAL::squared_distance(p3, p6)); check_squared_distance(S{p3, p2}, S{p6, p7}, CGAL::squared_distance(p3, p6)); @@ -321,7 +306,6 @@ private: const P p8 = p2 + lambda_8 * (p3 - p2); const double lambda_9 = r.get_double(-M, M); const P p9 = p2 + lambda_9 * (p3 - p2); - std::cout << "p2389 " << p2 << " " << p3 << " " << p8 << " " << p9 << std::endl; S s89(p8, p9); S s32(p3, p2); @@ -336,12 +320,15 @@ private: (std::min)(CGAL::squared_distance(p3, p8), CGAL::squared_distance(p3, p9)))); +#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*/) || @@ -352,7 +339,7 @@ private: S s1{p0, p1}, s2{p2, p3}; do_intersect_check(s1, s2); - // GTE ----------------------------------------------------------------------------------------- +#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()}}; @@ -361,6 +348,7 @@ private: 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) @@ -385,8 +373,8 @@ private: upper_bound = sqd; } - // bit ugly, but if constructions are not exact, building tp introduces some error - if(!exact) + // bit ugly, but if constructions are not exact, building `tp` introduces some error + if(epsilon != 0) upper_bound *= (1 + 1e-10); check_squared_distance_with_bound(s01, s23, upper_bound); @@ -397,14 +385,14 @@ private: { // Note : the value is not verified by hand std::cout << "Point - Ray" << std::endl; - check_squared_distance(p( -8, -7, 0), R{p(23, -27, 2), p( -17, 16, 2)}, 86.3685); + check_squared_distance_with_bound(p( -8, -7, 0), R{p(23, -27, 2), p( -17, 16, 2)}, 86.368512613); } void R_R() { // Note : the values are not verified by hand std::cout << "Ray - Ray" << std::endl; - check_squared_distance(R{p( 0, 0, 30), p( 0, 30, 30)}, R{p(100, -100, 0), p( 200, 1, 0)}, 20899.5); + check_squared_distance_with_bound(R{p( 0, 0, 30), p( 0, 30, 30)}, R{p(100, -100, 0), p( 200, 1, 0)}, 20899.504975002); check_squared_distance(R{p( 1, 0, 0), p( 0, 0, 0)}, R{p( 1, 3, 3), p( 0, 0, 3)}, 9); check_squared_distance(R{p( 0, 0, 0), p( 1, 0, 0)}, R{p( 0, 0, 2), p( -1, 0, 2)}, 4); } @@ -413,15 +401,15 @@ private: { // Note : the values are not verified by hand std::cout << "Segment - Ray" << std::endl; - check_squared_distance(S{p( 0, 0, 30), p( 0, 30, 30)}, R{p(100, -100, 0), p( 200, 1, 0)}, 20899.5); + check_squared_distance_with_bound(S{p( 0, 0, 30), p( 0, 30, 30)}, R{p(100, -100, 0), p( 200, 1, 0)}, 20899.504975002); } void R_L() { // Note : the values are not verified by hand std::cout << "Ray - Line" << std::endl; + check_squared_distance_with_bound(R{p( 0, 0, 30), p( 0, 30, 30)}, L{p(100, -100, 0), p( 200, 1, 0)}, 20899.504975002); check_squared_distance(R{p(10, 0, 0), p( 20, 0, 0)}, L{p( 0, 0, 3), p( 0, 3, 3)}, 109); - check_squared_distance(R{p( 0, 0, 30), p( 0, 30, 30)}, L{p(100, -100, 0), p( 200, 1, 0)}, 20899.5); check_squared_distance(R{p( 1, 0, 0), p( 0, 0, 0)}, L{p( 1, 3, 3), p( 0, 0, 3)}, 9); check_squared_distance(R{p( 0, 0, 0), p( 1, 0, 0)}, L{p( 0, 0, 2), p( -1, 0, 2)}, 4); } @@ -484,7 +472,7 @@ private: std::cout << "Plane - Plane" << std::endl; Pl p1(0, 1, 0, 0); typename K::Vector_3 v = -p1.orthogonal_vector(); - v /= CGAL::sqrt(v.squared_length()); + v /= CGAL::approximate_sqrt(v.squared_length()); Pl p2 = Pl(0,-1,0,6); check_squared_distance(p1,p2, 36); check_squared_distance(Pl(-2, 1, 1, 0), Pl(2, 1, 3, 0), 0); @@ -620,32 +608,19 @@ private: T tr1{p0, p1, p2}, tr2{p3, p4, p5}; do_intersect_check(tr1, tr2); - // PQP ----------------------------------------------------------------------------------------- -// CGAL::PQP_Checker PQP_checker; -// auto pqp_res = PQP_checker(p0, p1, p2, p3, p4, p5); -// std::cout << "dist (PQP) : " << pqp_res << std::endl; - - // GTE ----------------------------------------------------------------------------------------- -// gte::DCPQuery, gte::Triangle3 > GTE_TT_checker; -// gte::Triangle3 gte_tr1{{CGAL::to_double(p0.x()), CGAL::to_double(p0.y()), CGAL::to_double(p0.z())}, -// {CGAL::to_double(p1.x()), CGAL::to_double(p1.y()), CGAL::to_double(p1.z())}, -// {CGAL::to_double(p2.x()), CGAL::to_double(p2.y()), CGAL::to_double(p2.z())}}; -// gte::Triangle3 gte_tr2{{CGAL::to_double(p3.x()), CGAL::to_double(p3.y()), CGAL::to_double(p3.z())}, -// {CGAL::to_double(p4.x()), CGAL::to_double(p4.y()), CGAL::to_double(p4.z())}, -// {CGAL::to_double(p5.x()), CGAL::to_double(p5.y()), CGAL::to_double(p5.z())}}; -// auto gte_res = GTE_TT_checker(gte_tr1, gte_tr2); - +#ifdef CGAL_USE_GTE_AS_SANITY_CHECK gte::DCPQuery, 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; -// std::cout << "diff CGAL PQP : " << (pqp_res - CGAL::squared_distance(tr1, tr2)) << std::endl; - if(!are_equal(CGAL::squared_distance(tr1, tr2), gte_res.sqrDistance)) - std::cerr <<"uh-oh" << std::endl; + 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 } } @@ -654,28 +629,29 @@ public: { std::cout << "Kernel: " << typeid(K).name() << std::endl; -// P_P(); -// P_S(); -// P_R(); -// P_L(); -// P_T(); -// P_Pl(); + P_P(); + P_S(); + P_R(); + P_L(); + P_T(); + P_Pl(); + P_Tet(); S_S(); -// S_R(); -// S_L(); -// S_Pl(); + S_R(); + S_L(); + S_Pl(); -// R_R(); -// R_L(); -// R_Pl(); + R_R(); + R_L(); + R_Pl(); -// L_L(); -// L_Pl(); + L_L(); + L_Pl(); T_T(); -// Pl_Pl(); + Pl_Pl(); } }; @@ -689,11 +665,15 @@ int main() CGAL::Random r; std::cout << "random seed = " << r.get_seed() << std::endl; -// Test >(r, false).run(); -// Test >(r, false).run(); + // @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, false).run(); - Test(r, true).run(); + const double epick_eps = 10 * std::numeric_limits::epsilon(); + Test(r, epick_eps).run(); + + Test(r, 0).run(); std::cout << "Done!" << std::endl;