Add some real tests for Seg_3-Seg_3 and Tri_3-Tri_3

This commit is contained in:
Mael Rouxel-Labbé 2021-03-12 11:54:19 +01:00
parent 29d8b296ed
commit 5278c3044f
1 changed files with 74 additions and 94 deletions

View File

@ -8,18 +8,15 @@
#include <CGAL/Random.h>
#include <CGAL/Timer.h>
#include <CGAL/tri_tri_distance.h>
// #define CGAL_USE_GTE_AS_SANITY_CHECK
#ifdef CGAL_USE_GTE_AS_SANITY_CHECK
#include <Mathematics/DistTriangle3Triangle3.h>
#include <Mathematics/DistSegmentSegment.h>
#include <PQP.h>
#endif
#include <cassert>
#include <iostream>
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 <typename Type>
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<FT>::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<FT>::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 <typename O2>
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 <typename O1, typename O2>
void do_intersect_check(const O1& o1, const O2& o2)
{
@ -153,14 +148,6 @@ private:
}
}
template <typename O1, typename O2>
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 <typename O1, typename O2>
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 <typename O1, typename O2>
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<FT, gte::Segment<3, FT>, 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<FT, gte::Segment<3, FT>, 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<K> PQP_checker;
// auto pqp_res = PQP_checker(p0, p1, p2, p3, p4, p5);
// std::cout << "dist (PQP) : " << pqp_res << std::endl;
// GTE -----------------------------------------------------------------------------------------
// gte::DCPQuery<double, gte::Triangle3<double>, gte::Triangle3<double> > GTE_TT_checker;
// gte::Triangle3<double> 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<double> 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<FT, gte::Triangle3<FT>, gte::Triangle3<FT> > GTE_TT_checker;
gte::Triangle3<FT> gte_tr1{{p0.x(), p0.y(), p0.z()}, {p1.x(), p1.y(), p1.z()}, {p2.x(), p2.y(), p2.z()}};
gte::Triangle3<FT> 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<CGAL::Simple_cartesian<double> >(r, false).run();
// Test<CGAL::Simple_homogeneous<double> >(r, false).run();
// @todo Some tests are too difficult for these kernels
// Test<CGAL::Simple_cartesian<double> >(r, 1e-5).run();
// Test<CGAL::Simple_homogeneous<double> >(r, 1e-5).run();
// Test<CGAL::Simple_cartesian<CGAL::Interval_nt<true> > >(r, 1e-5).run();
// Test<CGAL::Exact_predicates_inexact_constructions_kernel>(r, false).run();
Test<CGAL::Exact_predicates_exact_constructions_kernel>(r, true).run();
const double epick_eps = 10 * std::numeric_limits<double>::epsilon();
Test<CGAL::Exact_predicates_inexact_constructions_kernel>(r, epick_eps).run();
Test<CGAL::Exact_predicates_exact_constructions_kernel>(r, 0).run();
std::cout << "Done!" << std::endl;