From 68b48a715184140e156029c397feca1a89219a0a Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 14 Mar 2012 15:54:19 +0000 Subject: [PATCH 02/38] Commit current version To be reviewed fully. --- .gitignore | 10 +++ .../Bbox_3_Segment_3_do_intersect.h | 73 +++++++++++++++++-- .../bbox_other_do_intersect_test.cpp | 31 +++++++- 3 files changed, 107 insertions(+), 7 deletions(-) diff --git a/.gitignore b/.gitignore index 451529c1056..d1b7b65716b 100644 --- a/.gitignore +++ b/.gitignore @@ -163,6 +163,16 @@ Installation/auxiliary/gdb/test Installation/cmake/modules/*.tmp Installation/test/Installation/cgal_test Installation/test/Installation/deprecation_warning +Intersections_3/test/Intersections_3/CMakeLists.txt +Intersections_3/test/Intersections_3/bbox_other_do_intersect_test +Intersections_3/test/Intersections_3/cgal_test_with_cmake +Intersections_3/test/Intersections_3/circle_other +Intersections_3/test/Intersections_3/line_line +Intersections_3/test/Intersections_3/segment_segment +Intersections_3/test/Intersections_3/test_intersections_3 +Intersections_3/test/Intersections_3/triangle_3_triangle_3_intersection +Intersections_3/test/Intersections_3/triangle_other +Intersections_3/test/Intersections_3/triangle_other_intersection_test Jet_fitting_3/examples/Jet_fitting_3/*.exe Jet_fitting_3/examples/Jet_fitting_3/*.sln Jet_fitting_3/examples/Jet_fitting_3/*.vcproj diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 43e6052e9ad..294b212fd0b 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -28,6 +28,15 @@ // inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf +// This algorithm intersects the line with the x-, y-, and z-slabs of the +// bounding box, and computes the interval [t1, t2], in the +// parameterization of the line given by the segment (for t=0, that is the +// source of the segment, and for t=1 that is its target), where the line +// intersects the three slabs of the the bounding box. + +// For a segment, the intersection is non-empty iff +// [t1, t2] intersects [0, 1]. + namespace CGAL { namespace internal { @@ -41,6 +50,11 @@ namespace internal { const FT& bxmin, const FT& bymin, const FT& bzmin, const FT& bxmax, const FT& bymax, const FT& bzmax) { + // The following code encode t1 and t2 by: + // t1 = tmin/dmin + // t2 = tmax/dmax + // For the first lines, dmax==dmin and is not explicitly defined. + // ----------------------------------- // treat x coord // ----------------------------------- @@ -58,26 +72,47 @@ namespace internal { dmin = px - qx; } - if ( tmax < FT(0) || tmin > dmin ) + if ( tmax < FT(0) ) // test t2 < 0, for a segment or a ray + return false; + if ( tmin > dmin ) // test t1 > 1, for a segment return false; + // If the query is vertical for x, then check its x-coordinate is in + // the x-slab. + if( (px == qx) && // <=> (dmin == 0) + ( CGAL::sign(tmin) * CGAL::sign(tmax) ) > 0 ) return false; + + // If dmin == 0, at this point, [t1, t2] == ]-inf, +inf[, or t1 or t2 + // is a NaN. But the case with NaNs is treated as if the interval + // [t1, t2] was ]-inf, +inf[. + FT dmax = dmin; + + // set t1=max(t1, 0), for a segment or a ray if ( tmin < FT(0) ) { tmin = FT(0); dmin = FT(1); } + // set t2=min(t2, 1), for a segment if ( tmax > dmax ) { tmax = FT(1); dmax = FT(1); } + CGAL_assertion(dmin >= 0); + CGAL_assertion(dmax >= 0); + // ----------------------------------- // treat y coord // ----------------------------------- + FT d_, tmin_, tmax_; + // Say: + // tymin = tmin_ / d_ + // tymax = tmax_ / d_ if ( qy >= py ) { tmin_ = bymin - py; @@ -91,19 +126,42 @@ namespace internal { d_ = py - qy; } - if ( (dmin*tmax_) < (d_*tmin) || (dmax*tmin_) > (d_*tmax) ) - return false; + // If t1 > tymax || tymin > t2, return false. + if ( dmin > 0 ) { + if( (dmin*tmax_) < (d_*tmin) ) return false; + if( (dmax*tmin_) > (d_*tmax) ) return false; + } - if( (dmin*tmin_) > (d_*tmin) ) + // If the segment is vertical for y, then check its y coordinate is in + // the y-slab. + if( (py == qy) && // <=> dmin == 0 + ( sign(tmin_) * sign(tmax_) ) > 0 ) return false; + + // If tymin > t1, set t1 = tymin. + if( dmin == 0 || (dmin*tmin_) > (d_*tmin) ) { tmin = tmin_; dmin = d_; + if(tmin > dmin) return false; // if t1 > 1, for a segment + if(tmin < FT(0)) { + // set t1=max(t1, 0), for a ray or a segment + tmin = FT(0); + dmin = FT(1); + } } - if( (dmax*tmax_) < (d_*tmax) ) + // If tymax < t2, set t2 = tymax. + if( dmax == 0 || (dmax*tmax_) < (d_*tmax) ) { tmax = tmax_; dmax = d_; + if(tmax < FT(0)) return false; // if t2 < 0, for a segment or a ray + if ( tmax > dmax ) + { + // set t2=min(t2, 1), for a segment + tmax = FT(1); + dmax = FT(1); + } } // ----------------------------------- @@ -122,7 +180,10 @@ namespace internal { d_ = pz - qz; } - return ( (dmin*tmax_) >= (d_*tmin) && (dmax*tmin_) <= (d_*tmax) ); + return ( (dmin*tmax_) >= (d_*tmin) && + (dmax*tmin_) <= (d_*tmax) && + tmin_ < d_ && + tmax_ > FT(0) ); } template diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index ba3c575e73a..d8cd3d7ae64 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -156,7 +156,7 @@ bool test() Point pB(FT(5.), FT(10.), FT(100.)); Point pC(FT(5.), FT(10.), FT(-1.)); Point pE(FT(5.), FT(10.), FT(0.)); - + Segment s12(p1,p2); Segment s13(p1,p3); Segment s14(p1,p4); @@ -184,8 +184,37 @@ bool test() b &= test_aux(sBC,"sBC",bbox,true); b &= test_aux(sCE,"sCE",bbox,false); b &= test_aux(sEC,"sEC",bbox,false); + + CGAL::Bbox_3 bbox_elem_1834(2, 1.81818, 0.166666, + 2.18182, 2.18182, 0.333333); + Point source_1834(2, 2, 1); + Point target_1834(2, 2, 0.75); + Segment segment_query_1834( source_1834, target_1834 ); + Point source_1834a(2, 2, 0.5); + Point target_1834a(2, 2, 0.75); + Segment segment_query_1834a( source_1834a, target_1834a ); + + b &= test_aux(segment_query_1834, + "segment_query_1834", bbox_elem_1834, false); + b &= test_aux(segment_query_1834.opposite(), + "segment_query_1834.opposite()", bbox_elem_1834, false); + b &= test_aux(segment_query_1834a, + "segment_query_1834a", bbox_elem_1834, false); + b &= test_aux(segment_query_1834a.opposite(), + "segment_query_1834a.opposite()", bbox_elem_1834, false); + + CGAL::Bbox_3 bbox_elem_1834b(1.81818, 2, 0, + 2, 2.18182, 0.333333); + Segment segment_query_1834b(Point(2, 2, 1), + Point(2, 2, 0.75)); + b &= test_aux(segment_query_1834b, + "segment_query_1834b", bbox_elem_1834b, false); + b &= test_aux(segment_query_1834b.opposite(), + "segment_query_1834b.opposite()", bbox_elem_1834b, false); + + Ray r12(p1,p2); Ray r13(p1,p3); Ray r14(p1,p4); From 300057e17a9855dc934428d99d738b1b6b870129 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 14 Mar 2012 16:24:41 +0000 Subject: [PATCH 03/38] Test also with Filtered_kernel without static filters And count the percentage of intersection among random objects --- .../bbox_other_do_intersect_test.cpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index d8cd3d7ae64..82550c9f15b 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -99,20 +99,24 @@ void speed(const std::string& name) CGAL::Timer timer; timer.start(); + std::size_t success = 0; while ( timer.time() < 0.1 ) { for ( typename std::vector::iterator it = segment_vector.begin(); it != segment_vector.end() ; ++it ) { - do_intersect(bbox_small, *it); + success += do_intersect(bbox_small, *it); } ++nb_loops; } timer.stop(); + std::cout << std::fixed << std::setprecision(1); std::cout << "\tDo_intersect(bbox, " << name << "): " << (nb_loops*segment_vector.size()) / (timer.time()*1000) - << " computations / ms " << std::endl; + << " computations / ms " + << (success / ((0.+ nb_loops*segment_vector.size()) / 100)) + << "% of intersection" << std::endl; } template @@ -373,6 +377,11 @@ int main() b &= test >(); test_speed >(); + std::cout << std::endl << "Testing with Filtered_kernel > without static filters..." << std::endl ; + typedef CGAL::Filtered_kernel, false> Fk_no_static; + b &= test(); + test_speed(); + std::cout << std::endl << "Testing with Exact_predicates_inexact_constructions_kernel..." << std::endl ; b &= test(); test_speed(); From c8bc77e46d5a0489bb56fd757e31c97535e501e9 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 14 Mar 2012 16:41:03 +0000 Subject: [PATCH 04/38] Submit work in progress. The old code was really wrong! --- .../Bbox_3_Segment_3_do_intersect.h | 37 +++++++++++++++---- 1 file changed, 30 insertions(+), 7 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 294b212fd0b..2763a4eaf86 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -126,19 +126,22 @@ namespace internal { d_ = py - qy; } - // If t1 > tymax || tymin > t2, return false. - if ( dmin > 0 ) { - if( (dmin*tmax_) < (d_*tmin) ) return false; - if( (dmax*tmin_) > (d_*tmax) ) return false; - } + CGAL_assertion(d_ >= 0); // If the segment is vertical for y, then check its y coordinate is in // the y-slab. if( (py == qy) && // <=> dmin == 0 ( sign(tmin_) * sign(tmax_) ) > 0 ) return false; + // If t1 > tymax || tymin > t2, return false. + if ( dmin > 0 && d_ > 0) { + if( (dmin*tmax_) < (d_*tmin) ) return false; + if( (dmax*tmin_) > (d_*tmax) ) return false; + } + // If tymin > t1, set t1 = tymin. - if( dmin == 0 || (dmin*tmin_) > (d_*tmin) ) + if( dmin == 0 || + ( d_ > 0 && (dmin*tmin_) > (d_*tmin) ) ) { tmin = tmin_; dmin = d_; @@ -151,7 +154,8 @@ namespace internal { } // If tymax < t2, set t2 = tymax. - if( dmax == 0 || (dmax*tmax_) < (d_*tmax) ) + if( dmax == 0 || + ( d_ > 0 && (dmax*tmax_) < (d_*tmax) ) ) { tmax = tmax_; dmax = d_; @@ -164,9 +168,16 @@ namespace internal { } } + CGAL_assertion(dmin >= 0); + CGAL_assertion(dmax >= 0); + // ----------------------------------- // treat z coord // ----------------------------------- + + // Say: + // tzmin = tmin_ / d_ + // tzmax = tmax_ / d_ if ( qz >= pz ) { tmin_ = bzmin - pz; @@ -180,6 +191,18 @@ namespace internal { d_ = pz - qz; } + CGAL_assertion(d_ >= 0); + + // If t1 > tzmax, return false. + + // If tzmin > t2, return false. + + // If tzmin > t1, set t1 = tzmin. + + // If tzmax < t2, set t2 = tzmax. + + // Return (t1 <= 1 && t2 >= 0) + return ( (dmin*tmax_) >= (d_*tmin) && (dmax*tmin_) <= (d_*tmax) && tmin_ < d_ && From c0a37addfb8db88b52f7111bfac8fe084d9be63f Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 15 Mar 2012 11:11:47 +0000 Subject: [PATCH 05/38] Also test with Sc. Simple_cartesian is a non-filtering kernel whose FT is not a IEEE 754 type. That number type is interesting to test with. --- .../test/Intersections_3/bbox_other_do_intersect_test.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index 82550c9f15b..c76ab69df9d 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -369,6 +369,10 @@ int main() b &= test >(); test_speed >(); + std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; + b &= test >(); + test_speed >(); + std::cout << std::endl << "Testing with Cartesian..." << std::endl ; b &= test >(); test_speed >(); From 71d1ba01e17969457e89a9d03f52435516240b55 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 15 Mar 2012 13:27:06 +0000 Subject: [PATCH 06/38] This version should be correct. --- .../Bbox_3_Segment_3_do_intersect.h | 28 +++++++++---------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 2763a4eaf86..4f32f8e0bf5 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -32,7 +32,7 @@ // bounding box, and computes the interval [t1, t2], in the // parameterization of the line given by the segment (for t=0, that is the // source of the segment, and for t=1 that is its target), where the line -// intersects the three slabs of the the bounding box. +// intersects the three slabs of the bounding box. // For a segment, the intersection is non-empty iff // [t1, t2] intersects [0, 1]. @@ -108,8 +108,8 @@ namespace internal { // ----------------------------------- // treat y coord // ----------------------------------- - FT d_, tmin_, tmax_; + // Say: // tymin = tmin_ / d_ // tymax = tmax_ / d_ @@ -130,7 +130,7 @@ namespace internal { // If the segment is vertical for y, then check its y coordinate is in // the y-slab. - if( (py == qy) && // <=> dmin == 0 + if( (py == qy) && // <=> d_ == 0 ( sign(tmin_) * sign(tmax_) ) > 0 ) return false; // If t1 > tymax || tymin > t2, return false. @@ -193,20 +193,18 @@ namespace internal { CGAL_assertion(d_ >= 0); - // If t1 > tzmax, return false. + // If the query is vertical for z, then check its z-coordinate is in + // the z-slab. + if( (pz == qz) && // <=> (d_ == 0) + ( CGAL::sign(tmin_) * CGAL::sign(tmax_) ) > 0 ) return false; - // If tzmin > t2, return false. + // If t1 > tzmax || tzmin > t2, return false. + if ( dmin > 0 && d_ > 0) { + if( (dmin*tmax_) < (d_*tmin) ) return false; + if( (dmax*tmin_) > (d_*tmax) ) return false; + } - // If tzmin > t1, set t1 = tzmin. - - // If tzmax < t2, set t2 = tzmax. - - // Return (t1 <= 1 && t2 >= 0) - - return ( (dmin*tmax_) >= (d_*tmin) && - (dmax*tmin_) <= (d_*tmax) && - tmin_ < d_ && - tmax_ > FT(0) ); + return ( d_ == 0 || tmin_ <= d_ && tmax_ >= FT(0) ); // t1 <= 1 && t2 >= 0 } template From 8dbae2494a130b518992cbbb8efa9bb2156b2d93 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 15 Mar 2012 14:48:12 +0000 Subject: [PATCH 07/38] Beginning of an intensive test suite of do_intersect(Bbox_3, Segment_3) --- .../bbox_other_do_intersect_test.cpp | 167 ++++++++++++++++-- 1 file changed, 156 insertions(+), 11 deletions(-) diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index c76ab69df9d..9d14cd2d7bb 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -72,6 +72,132 @@ bool test_aux(const T& t, return (b == expected); } +template +bool test_case(const FT& px, const FT& py, const FT& pz, + const FT& qx, const FT& qy, const FT& qz, + FT bxmin, FT bymin, FT bzmin, + FT bxmax, FT bymax, FT bzmax, + const bool expected, + const bool change_signs = true, + const bool swap_coords = true, + const bool opposite_seg = true, + const bool translate = true, + const bool scale = true) +{ + bool b = true; + if(change_signs) { + b &= test_case( px, py, pz, + qx, qy, qz, + bxmin, bymin, bzmin, + bxmax, bymax, bzmax, expected, false); + b &= test_case(-px, py, pz, + -qx, qy, qz, + -bxmin, bymin, bzmin, + -bxmax, bymax, bzmax, expected, false); + b &= test_case( px, -py, pz, + qx, -qy, qz, + bxmin, -bymin, bzmin, + bxmax, -bymax, bzmax, expected, false); + b &= test_case( px, py, -pz, + qx, qy, -qz, + bxmin, bymin, -bzmin, + bxmax, bymax, -bzmax, expected, false); + b &= test_case(-px, -py, pz, + -qx, -qy, qz, + -bxmin, -bymin, bzmin, + -bxmax, -bymax, bzmax, expected, false); + b &= test_case( px, -py, -pz, + qx, -qy, -qz, + bxmin, -bymin, -bzmin, + bxmax, -bymax, -bzmax, expected, false); + b &= test_case(-px, py, -pz, + -qx, qy, -qz, + -bxmin, bymin, -bzmin, + -bxmax, bymax, -bzmax, expected, false); + b &= test_case(-px, -py, -pz, + -qx, -qy, -qz, + -bxmin, -bymin, -bzmin, + -bxmax, -bymax, -bzmax, expected, false); + } else if(swap_coords) { + // xyz + b &= test_case( px, py, pz, + qx, qy, qz, + bxmin, bymin, bzmin, + bxmax, bymax, bzmax, expected, + false, false); + // xzy + b &= test_case( px, pz, py, + qx, qz, qy, + bxmin, bzmin, bymin, + bxmax, bzmax, bymax, expected, + false, false); + // yxz + b &= test_case( py, px, pz, + qy, qx, qz, + bymin, bxmin, bzmin, + bymax, bxmax, bzmax, expected, + false, false); + // zxy + b &= test_case( pz, px, py, + qz, qx, qy, + bzmin, bxmin, bymin, + bzmax, bxmax, bymax, expected, + false, false); + + // yzx + b &= test_case( py, pz, px, + qy, qz, qx, + bymin, bzmin, bxmin, + bymax, bzmax, bxmax, expected, + false, false); + // zyx + b &= test_case( pz, py, px, + qz, qy, qx, + bzmin, bymin, bxmin, + bzmax, bymax, bxmax, expected, + false, false); + } else if(opposite_seg) { + b &= test_case(px, py, pz, + qx, qy, qz, + bxmin, bymin, bzmin, + bxmax, bymax, bzmax, expected, + false, false, false); + b &= test_case(qx, qy, qz, + px, py, pz, + bxmin, bymin, bzmin, + bxmax, bymax, bzmax, expected, + false, false, false); + } else { + using CGAL::do_intersect; + using CGAL::Bbox_3; + typedef typename K::Point_3 Point_3; + typedef typename K::Segment_3 Segment_3; + if(bxmin > bxmax) std::swap(bxmin, bxmax); + if(bymin > bymax) std::swap(bymin, bymax); + if(bzmin > bzmax) std::swap(bzmin, bzmax); + if(do_intersect(Bbox_3(bxmin, bymin, bzmin, + bxmax, bymax, bzmax), + Segment_3(Point_3(px, py, pz), + Point_3(qx, qy, qz))) != expected) + { + b = false; + CGAL::set_pretty_mode(std::cerr); + std::cerr << "Wrong result for do_intersect(" + << Bbox_3(bxmin, bymin, bzmin, + bxmax, bymax, bzmax) + << ",\n" + << " " + << Segment_3(Point_3(px, py, pz), + Point_3(qx, qy, qz)) + << ")\n" + << " it should have been " << std::boolalpha << expected + << std::endl; + } + } + return b; +} + template void speed(const std::string& name) { @@ -208,6 +334,14 @@ bool test() "segment_query_1834a", bbox_elem_1834, false); b &= test_aux(segment_query_1834a.opposite(), "segment_query_1834a.opposite()", bbox_elem_1834, false); + b &= test_case(2., 2., 1., + 2., 2., 0.75, + 1.81818, 2., 0., + 2., 2.18182, 0.333333, false); + b &= test_case(2., 2., 0.5, + 2., 2., 0.75, + 1.81818, 2., 0., + 2., 2.18182, 0.333333, false); CGAL::Bbox_3 bbox_elem_1834b(1.81818, 2, 0, 2, 2.18182, 0.333333); @@ -217,8 +351,11 @@ bool test() "segment_query_1834b", bbox_elem_1834b, false); b &= test_aux(segment_query_1834b.opposite(), "segment_query_1834b.opposite()", bbox_elem_1834b, false); + b &= test_case(2., 2., 1., + 2., 2., 0.75, + 1.81818, 2., 0., + 2., 2.18182, 0.333333, false); - Ray r12(p1,p2); Ray r13(p1,p3); Ray r14(p1,p4); @@ -338,22 +475,30 @@ bool test() Line line3(seg3); Line line4(seg4); - test_aux(seg2, "seg2", bbox2, false); - test_aux(seg3, "seg3", bbox3, true); - test_aux(seg4, "seg4", bbox4, false); + b &= test_aux(seg2, "seg2", bbox2, false); + b &= test_aux(seg3, "seg3", bbox3, true); + b &= test_aux(seg4, "seg4", bbox4, false); - test_aux(ray2, "ray2", bbox2, false); - test_aux(ray3, "ray3", bbox3, true); - test_aux(ray4, "ray4", bbox4, false); + b &= test_aux(ray2, "ray2", bbox2, false); + b &= test_aux(ray3, "ray3", bbox3, true); + b &= test_aux(ray4, "ray4", bbox4, false); - test_aux(line2, "line2", bbox2, false); - test_aux(line3, "line3", bbox3, true); - test_aux(line4, "line4", bbox4, false); + b &= test_aux(line2, "line2", bbox2, false); + b &= test_aux(line3, "line3", bbox3, true); + b &= test_aux(line4, "line4", bbox4, false); // Use do_intersect(bbox,bbox) CGAL::do_intersect(bbox2,bbox4); - return b; + b &= test_case(1., 1., 0., + 1., 1., 1., + 0., 0., 0., + 1., 1., 1., true); + b &= test_case(1., 1., 0., + 2., 2., 2., + 0., 0., 0., + 1., 1., 1., true); + return b; } int main() From 611a2c439bdb3692832681b5db3b42a6ef115cd5 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 15 Mar 2012 15:09:05 +0000 Subject: [PATCH 08/38] Naively translate and scale the inputs, for test of do_intersect(Bbox, ..) --- .../bbox_other_do_intersect_test.cpp | 63 +++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index 9d14cd2d7bb..98fc720e6c8 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -168,6 +168,69 @@ bool test_case(const FT& px, const FT& py, const FT& pz, bxmin, bymin, bzmin, bxmax, bymax, bzmax, expected, false, false, false); + } else if(translate) { + b &= test_case(px, py, pz, + qx, qy, qz, + bxmin, bymin, bzmin, + bxmax, bymax, bzmax, expected, + false, false, false, false); + FT dx = 10, dy = 20, dz = 30; + b &= test_case(dx + px, dy + py, dz + pz, + dx + qx, dy + qy, dz + qz, + dx + bxmin, dy + bymin, dz + bzmin, + dx + bxmax, dy + bymax, dz + bzmax, expected, + false, false, false, false); + dx = (1 >> 10), dy = dx, dz = dx; + b &= test_case(dx + px, dy + py, dz + pz, + dx + qx, dy + qy, dz + qz, + dx + bxmin, dy + bymin, dz + bzmin, + dx + bxmax, dy + bymax, dz + bzmax, expected, + false, false, false, false); + dx = -(1 >> 10), dy = dx, dz = dx; + b &= test_case(dx + px, dy + py, dz + pz, + dx + qx, dy + qy, dz + qz, + dx + bxmin, dy + bymin, dz + bzmin, + dx + bxmax, dy + bymax, dz + bzmax, expected, + false, false, false, false); + } else if(scale) { + b &= test_case(px, py, pz, + qx, qy, qz, + bxmin, bymin, bzmin, + bxmax, bymax, bzmax, expected, + false, false, false, false, false); + FT delta = 9; + b &= test_case(delta * qx, delta * qy, delta * qz, + delta * px, delta * py, delta * pz, + delta * bxmin, delta * bymin, delta * bzmin, + delta * bxmax, delta * bymax, delta * bzmax, expected, + false, false, false, false, false); + delta = ((unsigned int)-1 / 4); + b &= test_case(delta * qx, delta * qy, delta * qz, + delta * px, delta * py, delta * pz, + delta * bxmin, delta * bymin, delta * bzmin, + delta * bxmax, delta * bymax, delta * bzmax, expected, + false, false, false, false, false); + delta = ((unsigned int)-1 / 4); + delta = 1/delta; + b &= test_case(delta * qx, delta * qy, delta * qz, + delta * px, delta * py, delta * pz, + delta * bxmin, delta * bymin, delta * bzmin, + delta * bxmax, delta * bymax, delta * bzmax, expected, + false, false, false, false, false); + delta = 7; + delta /= 3; + b &= test_case(delta * qx, delta * qy, delta * qz, + delta * px, delta * py, delta * pz, + delta * bxmin, delta * bymin, delta * bzmin, + delta * bxmax, delta * bymax, delta * bzmax, expected, + false, false, false, false, false); + delta = 1; + delta /= 10; + b &= test_case(delta * qx, delta * qy, delta * qz, + delta * px, delta * py, delta * pz, + delta * bxmin, delta * bymin, delta * bzmin, + delta * bxmax, delta * bymax, delta * bzmax, expected, + false, false, false, false, false); } else { using CGAL::do_intersect; using CGAL::Bbox_3; From bd09875feadc29a411b9696cb5f8800a9a5bac43 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 15 Mar 2012 15:45:31 +0000 Subject: [PATCH 09/38] Add a few test cases. One of then use boost::math::nextafter to move around a critical case. That is strange that even non-exact FT can deal with that without filtering. --- .../bbox_other_do_intersect_test.cpp | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index 98fc720e6c8..bb9c5e7177e 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -37,6 +37,7 @@ #include +#include // for nextafter double random_in(const double a, @@ -557,10 +558,56 @@ bool test() 1., 1., 1., 0., 0., 0., 1., 1., 1., true); + b &= test_case(0.5, 0.5, -0.5, + 0.5, 0.5, 0.5, + -0.5, -0.5, -0.5, + 0.5, 0.5, 0.5, true); + float f = 0.5f; + double d = boost::math::nextafter(f, f+1); + double d2 = boost::math::nextafter(f, f-1); + b &= test_case(d, 0.5, -0.5, + d, 0.5, 0.5, + -0.5, -0.5, -0.5, + 0.5, 0.5, 0.5, false, + false, false, false, false, false); + b &= test_case(d2, 0.5, -0.5, + d, 0.5, 0.5, + -0.5, -0.5, -0.5, + 0.5, 0.5, 0.5, true, + false, false, false, false, false); + b &= test_case(1., 1., 0., 2., 2., 2., 0., 0., 0., 1., 1., 1., true); + b &= test_case(1., 1., 1., + 1., 1., 1., + 0., 0., 0., + 1., 1., 1., true); + b &= test_case(0.9, 0.9, 0.9, + 0.9, 0.9, 0.9, + 0., 0., 0., + 1., 1., 1., true); + b &= test_case(0., 0., 0., + 0., 0., 0., + 0., 0., 0., + 1., 1., 1., true); + b &= test_case(0.1, 0., 0.1, + 0.1, 0., 0.1, + 0., 0., 0., + 1., 1., 1., true); + b &= test_case(0.1, -0.1, 0.1, + 0.1, -0.1, 0.1, + 0., 0., 0., + 1., 1., 1., false); + b &= test_case(0.1, 0.1, 0.1, + 0.1, 0.1, 0.1, + 0., 0., 0., + 1., 1., 1., true); + b &= test_case(1., 1., 1.1, + 1., 1., 1.1, + 0., 0., 0., + 1., 1., 1., false); return b; } From 3d36a2e4c359e4e0618d908ee6e57fa714ca4301 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 15 Mar 2012 15:47:36 +0000 Subject: [PATCH 10/38] Display coordinates with full precision --- .../test/Intersections_3/bbox_other_do_intersect_test.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index bb9c5e7177e..9711119b045 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -247,6 +247,7 @@ bool test_case(const FT& px, const FT& py, const FT& pz, { b = false; CGAL::set_pretty_mode(std::cerr); + std::cerr.precision(17); std::cerr << "Wrong result for do_intersect(" << Bbox_3(bxmin, bymin, bzmin, bxmax, bymax, bzmax) From 58809c064a400bfe9f4fb1aac5fef4e7ccfcd266 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 16 Mar 2012 16:10:44 +0000 Subject: [PATCH 11/38] Full test suite of do_intersect(Bbox_3, Segment_3). I have tested with gcov that all branches of the predicates are tested. --- .../bbox_other_do_intersect_test.cpp | 76 +++++++++++++++---- 1 file changed, 60 insertions(+), 16 deletions(-) diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index 9711119b045..0356defd69b 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -322,6 +322,49 @@ void test_speed() speed("line"); } +template +bool intensive_test(bool exact_predicates = true) +{ + bool b = true; + + // Test vertical segments + for(double x = 4.; x <= 7.; x+=1.) + for(double ymin = 0.; ymin <= 7.; ymin+=1.) + for(double ymax = ymin; ymax <= 7.; ymax+=1.) + { + const bool expected = + x >= -5. && x<= 5. && + ymin <= 5. && ymax >= -5; + b &= test_case(x, ymin, 0., + x, ymax, 0., + -5., -5., -5., 5., 5., 5., expected); + } + // Test slanted segments + for(double x = -7.; x <= 6.; x+=1.) + for(double y = -1.; y <= 6.; y+=1.) + { + const bool expected = + x >= -6. && x <= 5. && + y >= -6. && y <= 5. && + y <= x + 10. && y >= x - 10.; + b &= test_case(x, y, 0., + x + 1., y + 1., 0., + -5., -5., -5., 5., 5., 5., expected); + } + for(double x = -9.; x <= 6.; x+=1.) + for(double y = -3.; y <= 6.; y+=1.) + { + const bool expected = + x >= -8. && x <= 5. && + y >= -7. && y <= 5. && + 3 * y <= 2 * x + 25. && 3 * y >= 2 * x - 25.; + b &= test_case(x, y, 0., + x + 3., y + 2., 0., + -5., -5., -5., 5., 5., 5., expected); + } + return b; +} + template bool test() { @@ -612,43 +655,44 @@ bool test() return b; } +template +bool test_kernel(bool exact_predicates = true, K k = K()) +{ + bool b = test() && + intensive_test(exact_predicates); + test_speed >(); + return b; +} + int main() { srand(0); std::cout << std::setprecision(5); std::cout << "Testing with Simple_cartesian..." << std::endl ; - bool b = test >(); - test_speed >(); + bool b = test_kernel >(false); std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; - b &= test >(); - test_speed >(); + b &= test_kernel >(false); std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; - b &= test >(); - test_speed >(); + b &= test_kernel >(false); std::cout << std::endl << "Testing with Cartesian..." << std::endl ; - b &= test >(); - test_speed >(); + b &= test_kernel >(false); std::cout << std::endl << "Testing with Cartesian..." << std::endl ; - b &= test >(); - test_speed >(); + b &= test_kernel >(false); std::cout << std::endl << "Testing with Filtered_kernel > without static filters..." << std::endl ; typedef CGAL::Filtered_kernel, false> Fk_no_static; - b &= test(); - test_speed(); + b &= test_kernel(); std::cout << std::endl << "Testing with Exact_predicates_inexact_constructions_kernel..." << std::endl ; - b &= test(); - test_speed(); + b &= test_kernel(); std::cout << std::endl << "Testing with Exact_predicates_exact_constructions_kernel..." << std::endl ; - b &= test(); - test_speed(); + b &= test_kernel(); if ( b ) return EXIT_SUCCESS; From b5703e9cf8f9b4eed50d3ca7a39d1a29d874c5b1 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 16 Mar 2012 16:14:47 +0000 Subject: [PATCH 12/38] Add a note about a test that is already covered by previous test. --- .../internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 4f32f8e0bf5..1fd921605a3 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -81,6 +81,9 @@ namespace internal { // the x-slab. if( (px == qx) && // <=> (dmin == 0) ( CGAL::sign(tmin) * CGAL::sign(tmax) ) > 0 ) return false; + // Note: for a segment the condition sign(tmin)*sign(tmax) > 0 has + // already been tested by the two previous tests tmax<0 || tmin>dmin + // (with dmin==0). // If dmin == 0, at this point, [t1, t2] == ]-inf, +inf[, or t1 or t2 // is a NaN. But the case with NaNs is treated as if the interval From 1aa69b5a8bde23c92d6ccd5badb42e9b16acab90 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Mon, 19 Mar 2012 12:07:20 +0000 Subject: [PATCH 13/38] Factorize the _aux function for do_intersect(BBox_3, Ray_3|Segment_3) --- .../Bbox_3_Ray_3_do_intersect.h | 3 +- .../Bbox_3_Segment_3_do_intersect.h | 29 +++++++++++-------- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h index bc189388db9..baae9924b8a 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h @@ -24,6 +24,7 @@ #include #include +#include // inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf @@ -135,7 +136,7 @@ namespace internal { const Point_3& source = ray.source(); const Point_3& point_on_ray = ray.second_point(); - return bbox_ray_do_intersect_aux( + return do_intersect_bbox_segment_aux( source.x(), source.y(), source.z(), point_on_ray.x(), point_on_ray.y(), point_on_ray.z(), FT(bbox.xmin()), FT(bbox.ymin()), FT(bbox.zmin()), diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 1fd921605a3..301cd8a3a84 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -41,7 +41,9 @@ namespace CGAL { namespace internal { - template +template inline bool do_intersect_bbox_segment_aux( @@ -72,15 +74,16 @@ namespace internal { dmin = px - qx; } - if ( tmax < FT(0) ) // test t2 < 0, for a segment or a ray + if ( bounded_0 && tmax < FT(0) ) // test t2 < 0, for a segment or a ray return false; - if ( tmin > dmin ) // test t1 > 1, for a segment + if ( bounded_1 && tmin > dmin ) // test t1 > 1, for a segment return false; // If the query is vertical for x, then check its x-coordinate is in // the x-slab. if( (px == qx) && // <=> (dmin == 0) - ( CGAL::sign(tmin) * CGAL::sign(tmax) ) > 0 ) return false; + ( CGAL::sign(tmin) * CGAL::sign(tmax) ) > 0 ) + return false; // Note: for a segment the condition sign(tmin)*sign(tmax) > 0 has // already been tested by the two previous tests tmax<0 || tmin>dmin // (with dmin==0). @@ -92,14 +95,14 @@ namespace internal { FT dmax = dmin; // set t1=max(t1, 0), for a segment or a ray - if ( tmin < FT(0) ) + if ( bounded_0 && tmin < FT(0) ) { tmin = FT(0); dmin = FT(1); } // set t2=min(t2, 1), for a segment - if ( tmax > dmax ) + if ( bounded_1 && tmax > dmax ) { tmax = FT(1); dmax = FT(1); @@ -148,8 +151,8 @@ namespace internal { { tmin = tmin_; dmin = d_; - if(tmin > dmin) return false; // if t1 > 1, for a segment - if(tmin < FT(0)) { + if(bounded_1 && tmin > dmin) return false; // if t1 > 1, for a segment + if(bounded_0 && tmin < FT(0)) { // set t1=max(t1, 0), for a ray or a segment tmin = FT(0); dmin = FT(1); @@ -162,8 +165,8 @@ namespace internal { { tmax = tmax_; dmax = d_; - if(tmax < FT(0)) return false; // if t2 < 0, for a segment or a ray - if ( tmax > dmax ) + if( bounded_0 && tmax < FT(0)) return false; // if t2 < 0, for a segment or a ray + if( bounded_1 && tmax > dmax ) { // set t2=min(t2, 1), for a segment tmax = FT(1); @@ -207,7 +210,9 @@ namespace internal { if( (dmax*tmin_) > (d_*tmax) ) return false; } - return ( d_ == 0 || tmin_ <= d_ && tmax_ >= FT(0) ); // t1 <= 1 && t2 >= 0 + return ( d_ == 0 || + (!bounded_1 || tmin_ <= d_) && + (!bounded_0 || tmax_ >= FT(0)) ); // t1 <= 1 && t2 >= 0 } template @@ -221,7 +226,7 @@ namespace internal { const Point_3& source = segment.source(); const Point_3& target = segment.target(); - return do_intersect_bbox_segment_aux( + return do_intersect_bbox_segment_aux( source.x(), source.y(), source.z(), target.x(), target.y(), target.z(), FT(bbox.xmin()), FT(bbox.ymin()), FT(bbox.zmin()), From b0cfb5bc1f5e81167bf8e247f144b17270ddd123 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Mon, 19 Mar 2012 14:12:05 +0000 Subject: [PATCH 14/38] Fix a stupid copy-paste error --- .../test/Intersections_3/bbox_other_do_intersect_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index 0356defd69b..05faf57a827 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -394,7 +394,7 @@ bool test() Point pB(FT(5.), FT(10.), FT(100.)); Point pC(FT(5.), FT(10.), FT(-1.)); Point pE(FT(5.), FT(10.), FT(0.)); - + Segment s12(p1,p2); Segment s13(p1,p3); Segment s14(p1,p4); @@ -660,7 +660,7 @@ bool test_kernel(bool exact_predicates = true, K k = K()) { bool b = test() && intensive_test(exact_predicates); - test_speed >(); + test_speed(); return b; } From 0358937f0191b4ddaeeffe71a7a84176c2648b9e Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 20 Mar 2012 16:56:13 +0000 Subject: [PATCH 15/38] Code optimized for x-axis --- .../Bbox_3_Ray_3_do_intersect.h | 9 +- .../Bbox_3_Segment_3_do_intersect.h | 127 +++++++++++++----- .../bbox_other_do_intersect_test.cpp | 116 ++++++++-------- 3 files changed, 157 insertions(+), 95 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h index baae9924b8a..6a402e1cd4d 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h @@ -139,8 +139,13 @@ namespace internal { return do_intersect_bbox_segment_aux( source.x(), source.y(), source.z(), point_on_ray.x(), point_on_ray.y(), point_on_ray.z(), - FT(bbox.xmin()), FT(bbox.ymin()), FT(bbox.zmin()), - FT(bbox.xmax()), FT(bbox.ymax()), FT(bbox.zmax()) ); + bbox.xmin(), bbox.ymin(), bbox.zmin(), + bbox.xmax(), bbox.ymax(), bbox.zmax() ); + // const CGAL::cpp0x::array rray = {source.x(), source.y(), source.z(), + // point_on_ray.x(), point_on_ray.y(), point_on_ray.z() }; + // return do_intersect_bbox_segment_aux + // (rray, + // *reinterpret_cast*>(&*bbox.cartesian_begin()) ); } } // namespace internal diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 301cd8a3a84..b370a0ebe8f 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -41,16 +41,17 @@ namespace CGAL { namespace internal { -template + template +// __attribute__ ((noinline)) inline bool do_intersect_bbox_segment_aux( const FT& px, const FT& py, const FT& pz, const FT& qx, const FT& qy, const FT& qz, - const FT& bxmin, const FT& bymin, const FT& bzmin, - const FT& bxmax, const FT& bymax, const FT& bzmax) + const double& bxmin, const double& bymin, const double& bzmin, + const double& bxmax, const double& bymax, const double& bzmax) { // The following code encode t1 and t2 by: // t1 = tmin/dmin @@ -60,29 +61,57 @@ template = px ) { - tmin = bxmin - px; - tmax = bxmax - px; - dmin = qx - px; + if(bounded_0 && px > bxmax) return false; // segment on the right of bbox + if(bounded_1 && qx < bxmin) return false; // segment on the left of bbox + + if(bounded_1 && bxmax > qx) { + tmax = 1; + dmax = 1; + } else { + tmax = bxmax - px; + dmax = qx - px; + } + + if(bounded_0 && bxmin < px) // tmin < 0 means px is in the x-range of bbox + { + tmin = 0; + dmin = 1; + } else { + tmin = bxmin - px; + dmin = qx - px; + } } else { - tmin = px - bxmax; - tmax = px - bxmin; - dmin = px - qx; - } + if(bounded_1 && qx > bxmax) return false; // segment on the right of bbox + if(bounded_0 && px < bxmin) return false; // segment on the left of bbox - if ( bounded_0 && tmax < FT(0) ) // test t2 < 0, for a segment or a ray - return false; - if ( bounded_1 && tmin > dmin ) // test t1 > 1, for a segment - return false; + if(bounded_1 && bxmin < qx) { + tmax = 1; + dmax = 1; + } else { + tmax = px - bxmin; + dmax = px - qx; + } + + if(bounded_0 && px < bxmax) // tmin < 0 means px is in the x-range of bbox + { + tmin = 0; + dmin = 1; + } else { + tmin = px - bxmax; + dmin = px - qx; + } + } // If the query is vertical for x, then check its x-coordinate is in // the x-slab. - if( (px == qx) && // <=> (dmin == 0) - ( CGAL::sign(tmin) * CGAL::sign(tmax) ) > 0 ) + if( (px == qx) && // <=> (dmin == 0) + (! (bounded_0 && bounded_1) ) && + ( CGAL::sign(tmin) * CGAL::sign(tmax) ) > 0 ) return false; // Note: for a segment the condition sign(tmin)*sign(tmax) > 0 has // already been tested by the two previous tests tmax<0 || tmin>dmin @@ -92,22 +121,6 @@ template dmax ) - { - tmax = FT(1); - dmax = FT(1); - } - CGAL_assertion(dmin >= 0); CGAL_assertion(dmax >= 0); @@ -215,6 +228,40 @@ template = FT(0)) ); // t1 <= 1 && t2 >= 0 } + template + inline + bool + do_intersect_bbox_segment_aux(const CGAL::cpp0x::array seg, + const CGAL::cpp0x::array box) + + { + const FT& px = seg[0]; + const FT& py = seg[1]; + const FT& pz = seg[2]; + const FT& qx = seg[3]; + const FT& qy = seg[4]; + const FT& qz = seg[5]; + const double& bxmin = box[0]; + const double& bymin = box[1]; + const double& bzmin = box[2]; + const double& bxmax = box[3]; + const double& bymax = box[4]; + const double& bzmax = box[5]; + // for(int i = 0; i < 3; ++i) { + // const int sign = seg[3+i] > seg[i]; // (qx > px)? + // if(bounded_0 && seg[3*(1-sign) + i] > box[3+i]) return false; // segment on the right of bbox + // if(bounded_1 && seg[3*sign + i] < box[i]) return false; // segment on the left of bbox + // } + + return do_intersect_bbox_segment_aux + (px, py, pz, + qy, qy, qz, + bxmin, bymin, bymax, + bxmax, bymax, bzmax); + } + template bool do_intersect(const typename K::Segment_3& segment, const CGAL::Bbox_3& bbox, @@ -229,8 +276,14 @@ template ( source.x(), source.y(), source.z(), target.x(), target.y(), target.z(), - FT(bbox.xmin()), FT(bbox.ymin()), FT(bbox.zmin()), - FT(bbox.xmax()), FT(bbox.ymax()), FT(bbox.zmax()) ); + bbox.xmin(), bbox.ymin(), bbox.zmin(), + bbox.xmax(), bbox.ymax(), bbox.zmax() ); + + // const CGAL::cpp0x::array seg = {source.x(), source.y(), source.z(), + // target.x(), target.y(), target.z() }; + // return do_intersect_bbox_segment_aux + // ( seg, + // *reinterpret_cast*>(&*bbox.cartesian_begin()) ); } template diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index 05faf57a827..07ecc7de679 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -62,7 +62,7 @@ template bool test_aux(const T& t, const std::string& name, const CGAL::Bbox_3& bbox, - bool expected) + bool expected, bool exact_predicates = false) { bool b = CGAL::do_intersect(t,bbox); @@ -79,7 +79,9 @@ bool test_case(const FT& px, const FT& py, const FT& pz, const FT& qx, const FT& qy, const FT& qz, FT bxmin, FT bymin, FT bzmin, FT bxmax, FT bymax, FT bzmax, - const bool expected, + const bool expected, + const bool exact_k = false, + const bool exactness_issue = false, const bool change_signs = true, const bool swap_coords = true, const bool opposite_seg = true, @@ -91,146 +93,146 @@ bool test_case(const FT& px, const FT& py, const FT& pz, b &= test_case( px, py, pz, qx, qy, qz, bxmin, bymin, bzmin, - bxmax, bymax, bzmax, expected, false); + bxmax, bymax, bzmax, expected, exact_k, exactness_issue, false); b &= test_case(-px, py, pz, -qx, qy, qz, -bxmin, bymin, bzmin, - -bxmax, bymax, bzmax, expected, false); + -bxmax, bymax, bzmax, expected, exact_k, exactness_issue, false); b &= test_case( px, -py, pz, qx, -qy, qz, bxmin, -bymin, bzmin, - bxmax, -bymax, bzmax, expected, false); + bxmax, -bymax, bzmax, expected, exact_k, exactness_issue, false); b &= test_case( px, py, -pz, qx, qy, -qz, bxmin, bymin, -bzmin, - bxmax, bymax, -bzmax, expected, false); + bxmax, bymax, -bzmax, expected, exact_k, exactness_issue, false); b &= test_case(-px, -py, pz, -qx, -qy, qz, -bxmin, -bymin, bzmin, - -bxmax, -bymax, bzmax, expected, false); + -bxmax, -bymax, bzmax, expected, exact_k, exactness_issue, false); b &= test_case( px, -py, -pz, qx, -qy, -qz, bxmin, -bymin, -bzmin, - bxmax, -bymax, -bzmax, expected, false); + bxmax, -bymax, -bzmax, expected, exact_k, exactness_issue, false); b &= test_case(-px, py, -pz, -qx, qy, -qz, -bxmin, bymin, -bzmin, - -bxmax, bymax, -bzmax, expected, false); + -bxmax, bymax, -bzmax, expected, exact_k, exactness_issue, false); b &= test_case(-px, -py, -pz, -qx, -qy, -qz, -bxmin, -bymin, -bzmin, - -bxmax, -bymax, -bzmax, expected, false); + -bxmax, -bymax, -bzmax, expected, exact_k, exactness_issue, false); } else if(swap_coords) { // xyz b &= test_case( px, py, pz, qx, qy, qz, bxmin, bymin, bzmin, - bxmax, bymax, bzmax, expected, + bxmax, bymax, bzmax, expected, exact_k, exactness_issue, false, false); // xzy b &= test_case( px, pz, py, qx, qz, qy, bxmin, bzmin, bymin, - bxmax, bzmax, bymax, expected, + bxmax, bzmax, bymax, expected, exact_k, exactness_issue, false, false); // yxz b &= test_case( py, px, pz, qy, qx, qz, bymin, bxmin, bzmin, - bymax, bxmax, bzmax, expected, + bymax, bxmax, bzmax, expected, exact_k, exactness_issue, false, false); // zxy b &= test_case( pz, px, py, qz, qx, qy, bzmin, bxmin, bymin, - bzmax, bxmax, bymax, expected, + bzmax, bxmax, bymax, expected, exact_k, exactness_issue, false, false); // yzx b &= test_case( py, pz, px, qy, qz, qx, bymin, bzmin, bxmin, - bymax, bzmax, bxmax, expected, + bymax, bzmax, bxmax, expected, exact_k, exactness_issue, false, false); // zyx b &= test_case( pz, py, px, qz, qy, qx, bzmin, bymin, bxmin, - bzmax, bymax, bxmax, expected, + bzmax, bymax, bxmax, expected, exact_k, exactness_issue, false, false); } else if(opposite_seg) { b &= test_case(px, py, pz, qx, qy, qz, bxmin, bymin, bzmin, - bxmax, bymax, bzmax, expected, + bxmax, bymax, bzmax, expected, exact_k, exactness_issue, false, false, false); b &= test_case(qx, qy, qz, px, py, pz, bxmin, bymin, bzmin, - bxmax, bymax, bzmax, expected, + bxmax, bymax, bzmax, expected, exact_k, exactness_issue, false, false, false); } else if(translate) { b &= test_case(px, py, pz, qx, qy, qz, bxmin, bymin, bzmin, - bxmax, bymax, bzmax, expected, + bxmax, bymax, bzmax, expected, exact_k, exactness_issue, false, false, false, false); FT dx = 10, dy = 20, dz = 30; b &= test_case(dx + px, dy + py, dz + pz, dx + qx, dy + qy, dz + qz, dx + bxmin, dy + bymin, dz + bzmin, - dx + bxmax, dy + bymax, dz + bzmax, expected, + dx + bxmax, dy + bymax, dz + bzmax, expected, exact_k, exactness_issue, false, false, false, false); dx = (1 >> 10), dy = dx, dz = dx; b &= test_case(dx + px, dy + py, dz + pz, dx + qx, dy + qy, dz + qz, dx + bxmin, dy + bymin, dz + bzmin, - dx + bxmax, dy + bymax, dz + bzmax, expected, + dx + bxmax, dy + bymax, dz + bzmax, expected, exact_k, exactness_issue, false, false, false, false); dx = -(1 >> 10), dy = dx, dz = dx; b &= test_case(dx + px, dy + py, dz + pz, dx + qx, dy + qy, dz + qz, dx + bxmin, dy + bymin, dz + bzmin, - dx + bxmax, dy + bymax, dz + bzmax, expected, + dx + bxmax, dy + bymax, dz + bzmax, expected, exact_k, exactness_issue, false, false, false, false); } else if(scale) { b &= test_case(px, py, pz, qx, qy, qz, bxmin, bymin, bzmin, - bxmax, bymax, bzmax, expected, + bxmax, bymax, bzmax, expected, exact_k, exactness_issue, false, false, false, false, false); FT delta = 9; b &= test_case(delta * qx, delta * qy, delta * qz, delta * px, delta * py, delta * pz, delta * bxmin, delta * bymin, delta * bzmin, - delta * bxmax, delta * bymax, delta * bzmax, expected, + delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue, false, false, false, false, false); - delta = ((unsigned int)-1 / 4); + delta = (1 << 10); b &= test_case(delta * qx, delta * qy, delta * qz, delta * px, delta * py, delta * pz, delta * bxmin, delta * bymin, delta * bzmin, - delta * bxmax, delta * bymax, delta * bzmax, expected, + delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue, false, false, false, false, false); - delta = ((unsigned int)-1 / 4); + delta = (1 << 10); delta = 1/delta; b &= test_case(delta * qx, delta * qy, delta * qz, delta * px, delta * py, delta * pz, delta * bxmin, delta * bymin, delta * bzmin, - delta * bxmax, delta * bymax, delta * bzmax, expected, + delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue, false, false, false, false, false); delta = 7; - delta /= 3; + delta /= 2; b &= test_case(delta * qx, delta * qy, delta * qz, delta * px, delta * py, delta * pz, delta * bxmin, delta * bymin, delta * bzmin, - delta * bxmax, delta * bymax, delta * bzmax, expected, + delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue, false, false, false, false, false); delta = 1; - delta /= 10; + delta /= 8; b &= test_case(delta * qx, delta * qy, delta * qz, delta * px, delta * py, delta * pz, delta * bxmin, delta * bymin, delta * bzmin, - delta * bxmax, delta * bymax, delta * bzmax, expected, + delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue, false, false, false, false, false); } else { using CGAL::do_intersect; @@ -245,19 +247,21 @@ bool test_case(const FT& px, const FT& py, const FT& pz, Segment_3(Point_3(px, py, pz), Point_3(qx, qy, qz))) != expected) { - b = false; - CGAL::set_pretty_mode(std::cerr); - std::cerr.precision(17); - std::cerr << "Wrong result for do_intersect(" - << Bbox_3(bxmin, bymin, bzmin, - bxmax, bymax, bzmax) - << ",\n" - << " " - << Segment_3(Point_3(px, py, pz), - Point_3(qx, qy, qz)) - << ")\n" - << " it should have been " << std::boolalpha << expected - << std::endl; + if(!exactness_issue || exact_k) { + b = false; + CGAL::set_pretty_mode(std::cerr); + std::cerr.precision(17); + std::cerr << "Wrong result for do_intersect(" + << Bbox_3(bxmin, bymin, bzmin, + bxmax, bymax, bzmax) + << ",\n" + << " " + << Segment_3(Point_3(px, py, pz), + Point_3(qx, qy, qz)) + << ")\n" + << " it should have been " << std::boolalpha << expected + << std::endl; + } } } return b; @@ -291,7 +295,7 @@ void speed(const std::string& name) CGAL::Timer timer; timer.start(); std::size_t success = 0; - while ( timer.time() < 0.1 ) + while ( timer.time() < 5. ) { for ( typename std::vector::iterator it = segment_vector.begin(); it != segment_vector.end() ; ++it ) @@ -337,7 +341,7 @@ bool intensive_test(bool exact_predicates = true) ymin <= 5. && ymax >= -5; b &= test_case(x, ymin, 0., x, ymax, 0., - -5., -5., -5., 5., 5., 5., expected); + -5., -5., -5., 5., 5., 5., expected, exact_predicates); } // Test slanted segments for(double x = -7.; x <= 6.; x+=1.) @@ -349,7 +353,7 @@ bool intensive_test(bool exact_predicates = true) y <= x + 10. && y >= x - 10.; b &= test_case(x, y, 0., x + 1., y + 1., 0., - -5., -5., -5., 5., 5., 5., expected); + -5., -5., -5., 5., 5., 5., expected, exact_predicates); } for(double x = -9.; x <= 6.; x+=1.) for(double y = -3.; y <= 6.; y+=1.) @@ -360,13 +364,13 @@ bool intensive_test(bool exact_predicates = true) 3 * y <= 2 * x + 25. && 3 * y >= 2 * x - 25.; b &= test_case(x, y, 0., x + 3., y + 2., 0., - -5., -5., -5., 5., 5., 5., expected); + -5., -5., -5., 5., 5., 5., expected, exact_predicates); } return b; } template -bool test() +bool test(bool exact_kernel = false) { // types typedef typename K::FT FT; @@ -612,12 +616,12 @@ bool test() b &= test_case(d, 0.5, -0.5, d, 0.5, 0.5, -0.5, -0.5, -0.5, - 0.5, 0.5, 0.5, false, + 0.5, 0.5, 0.5, false, exact_kernel, true, false, false, false, false, false); b &= test_case(d2, 0.5, -0.5, d, 0.5, 0.5, -0.5, -0.5, -0.5, - 0.5, 0.5, 0.5, true, + 0.5, 0.5, 0.5, true, exact_kernel, true, false, false, false, false, false); b &= test_case(1., 1., 0., @@ -658,7 +662,7 @@ bool test() template bool test_kernel(bool exact_predicates = true, K k = K()) { - bool b = test() && + bool b = test(exact_predicates) && intensive_test(exact_predicates); test_speed(); return b; @@ -673,16 +677,16 @@ int main() bool b = test_kernel >(false); std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; - b &= test_kernel >(false); + b &= test_kernel >(true); std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; - b &= test_kernel >(false); + b &= test_kernel >(true); std::cout << std::endl << "Testing with Cartesian..." << std::endl ; b &= test_kernel >(false); std::cout << std::endl << "Testing with Cartesian..." << std::endl ; - b &= test_kernel >(false); + b &= test_kernel >(true); std::cout << std::endl << "Testing with Filtered_kernel > without static filters..." << std::endl ; typedef CGAL::Filtered_kernel, false> Fk_no_static; From e54df14afb74eed702b328da6bd2ba3a92b9c427 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 20 Mar 2012 17:47:55 +0000 Subject: [PATCH 16/38] Code optimized for all x-, y-, and z-axis --- .../Bbox_3_Segment_3_do_intersect.h | 209 +++++++++++------- 1 file changed, 133 insertions(+), 76 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index b370a0ebe8f..03a48ef9aff 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -110,9 +110,10 @@ namespace internal { // If the query is vertical for x, then check its x-coordinate is in // the x-slab. if( (px == qx) && // <=> (dmin == 0) - (! (bounded_0 && bounded_1) ) && - ( CGAL::sign(tmin) * CGAL::sign(tmax) ) > 0 ) - return false; + (! (bounded_0 && bounded_1) ) ) // do not check for a segment + { + if(px > bxmax || px < bxmin) return false; + } // Note: for a segment the condition sign(tmin)*sign(tmax) > 0 has // already been tested by the two previous tests tmax<0 || tmin>dmin // (with dmin==0). @@ -127,105 +128,161 @@ namespace internal { // ----------------------------------- // treat y coord // ----------------------------------- - FT d_, tmin_, tmax_; - - // Say: - // tymin = tmin_ / d_ - // tymax = tmax_ / d_ + FT dymin, tymin, tymax, dymax; if ( qy >= py ) { - tmin_ = bymin - py; - tmax_ = bymax - py; - d_ = qy - py; + if(bounded_0 && py > bymax) return false; // segment on the right of bbox + if(bounded_1 && qy < bymin) return false; // segment on the left of bbox + + if(bounded_1 && bymax > qy) { + tymax = 1; + dymax = 1; + } else { + tymax = bymax - py; + dymax = qy - py; + } + + if(bounded_0 && bymin < py) // tmin < 0 means py is in the y-range of bbox + { + tymin = 0; + dymin = 1; + } else { + tymin = bymin - py; + dymin = qy - py; + } } else { - tmin_ = py - bymax; - tmax_ = py - bymin; - d_ = py - qy; + if(bounded_1 && qy > bymax) return false; // segment on the right of bbox + if(bounded_0 && py < bymin) return false; // segment on the left of bbox + + if(bounded_1 && bymin < qy) { + tymax = 1; + dymax = 1; + } else { + tymax = py - bymin; + dymax = py - qy; + } + + if(bounded_0 && py < bymax) // tmin < 0 means py is in the y-range of bbox + { + tymin = 0; + dymin = 1; + } else { + tymin = py - bymax; + dymin = py - qy; + } } - CGAL_assertion(d_ >= 0); - - // If the segment is vertical for y, then check its y coordinate is in + // If the query is vertical for y, then check its y-coordinate is in // the y-slab. - if( (py == qy) && // <=> d_ == 0 - ( sign(tmin_) * sign(tmax_) ) > 0 ) return false; + if( (py == qy) && // <=> (dmin == 0) + (! (bounded_0 && bounded_1) ) ) // do not check for a segment + { + if(py > bymax || py < bymin) return false; + } + + // If dmin == 0, at this point, [t1, t2] == ]-inf, +inf[, or t1 or t2 + // is a NaN. But the case with NaNs is treated as if the interval + // [t1, t2] was ]-inf, +inf[. + + CGAL_assertion(dymin >= 0); + CGAL_assertion(dymax >= 0); + + // ----------------------------------- + // treat z coord + // ----------------------------------- + FT dzmin, tzmin, tzmax, dzmax; + if ( qz >= pz ) + { + if(bounded_0 && pz > bzmax) return false; // segment on the right of bbox + if(bounded_1 && qz < bzmin) return false; // segment on the left of bbox + + if(bounded_1 && bzmax > qz) { + tzmax = 1; + dzmax = 1; + } else { + tzmax = bzmax - pz; + dzmax = qz - pz; + } + + if(bounded_0 && bzmin < pz) // tmin < 0 means pz is in the z-range of bbox + { + tzmin = 0; + dzmin = 1; + } else { + tzmin = bzmin - pz; + dzmin = qz - pz; + } + } + else + { + if(bounded_1 && qz > bzmax) return false; // segment on the right of bbox + if(bounded_0 && pz < bzmin) return false; // segment on the left of bbox + + if(bounded_1 && bzmin < qz) { + tzmax = 1; + dzmax = 1; + } else { + tzmax = pz - bzmin; + dzmax = pz - qz; + } + + if(bounded_0 && pz < bzmax) // tmin < 0 means pz is in the z-range of bbox + { + tzmin = 0; + dzmin = 1; + } else { + tzmin = pz - bzmax; + dzmin = pz - qz; + } + } + + // If the querz is vertical for z, then check its z-coordinate is in + // the z-slab. + if( (pz == qz) && // <=> (dmin == 0) + (! (bounded_0 && bounded_1) ) ) // do not check for a segment + { + if(pz > bzmax || pz < bzmin) return false; + } + + // If dmin == 0, at this point, [t1, t2] == ]-inf, +inf[, or t1 or t2 + // is a NaN. But the case with NaNs is treated as if the interval + // [t1, t2] was ]-inf, +inf[. + + CGAL_assertion(dzmin >= 0); + CGAL_assertion(dzmax >= 0); // If t1 > tymax || tymin > t2, return false. - if ( dmin > 0 && d_ > 0) { - if( (dmin*tmax_) < (d_*tmin) ) return false; - if( (dmax*tmin_) > (d_*tmax) ) return false; - } + if( dymax > 0 && dmin > 0 && (dmin*tymax) < (dymax*tmin) ) return false; + if( dymin > 0 && dmax > 0 && (dmax*tymin) > (dymin*tmax) ) return false; // If tymin > t1, set t1 = tymin. if( dmin == 0 || - ( d_ > 0 && (dmin*tmin_) > (d_*tmin) ) ) + ( dymin > 0 && (dmin*tymin) > (dymin*tmin) ) ) { - tmin = tmin_; - dmin = d_; - if(bounded_1 && tmin > dmin) return false; // if t1 > 1, for a segment - if(bounded_0 && tmin < FT(0)) { - // set t1=max(t1, 0), for a ray or a segment - tmin = FT(0); - dmin = FT(1); - } + tmin = tymin; + dmin = dymin; } // If tymax < t2, set t2 = tymax. if( dmax == 0 || - ( d_ > 0 && (dmax*tmax_) < (d_*tmax) ) ) + ( dymax > 0 && (dmax*tymax) < (dymax*tmax) ) ) { - tmax = tmax_; - dmax = d_; - if( bounded_0 && tmax < FT(0)) return false; // if t2 < 0, for a segment or a ray - if( bounded_1 && tmax > dmax ) - { - // set t2=min(t2, 1), for a segment - tmax = FT(1); - dmax = FT(1); - } + tmax = tymax; + dmax = dymax; } CGAL_assertion(dmin >= 0); CGAL_assertion(dmax >= 0); - // ----------------------------------- - // treat z coord - // ----------------------------------- - - // Say: - // tzmin = tmin_ / d_ - // tzmax = tmax_ / d_ - if ( qz >= pz ) - { - tmin_ = bzmin - pz; - tmax_ = bzmax - pz; - d_ = qz - pz; - } - else - { - tmin_ = pz - bzmax; - tmax_ = pz - bzmin; - d_ = pz - qz; - } - - CGAL_assertion(d_ >= 0); - - // If the query is vertical for z, then check its z-coordinate is in - // the z-slab. - if( (pz == qz) && // <=> (d_ == 0) - ( CGAL::sign(tmin_) * CGAL::sign(tmax_) ) > 0 ) return false; - // If t1 > tzmax || tzmin > t2, return false. - if ( dmin > 0 && d_ > 0) { - if( (dmin*tmax_) < (d_*tmin) ) return false; - if( (dmax*tmin_) > (d_*tmax) ) return false; - } + if( dmin > 0 && dzmax > 0 && (dmin*tzmax) < (dzmax*tmin) ) return false; + if( dmax > 0 && dzmin > 0 && (dmax*tzmin) > (dzmin*tmax) ) return false; - return ( d_ == 0 || - (!bounded_1 || tmin_ <= d_) && - (!bounded_0 || tmax_ >= FT(0)) ); // t1 <= 1 && t2 >= 0 + return ( dzmin == 0 || + (!bounded_1 || tzmin <= dzmin) && + (!bounded_0 || tzmax >= FT(0)) ); // t1 <= 1 && t2 >= 0 } template Date: Wed, 21 Mar 2012 11:13:19 +0000 Subject: [PATCH 17/38] Less tests of sign of expressions I have managed to transform most of the tests to simple comparison of input coordinates. That will ease the writing of static filters. --- .../Bbox_3_Segment_3_do_intersect.h | 35 +++++++++++-------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 03a48ef9aff..3ce8b3efd61 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -113,10 +113,9 @@ namespace internal { (! (bounded_0 && bounded_1) ) ) // do not check for a segment { if(px > bxmax || px < bxmin) return false; + // Note: for a segment the condition has already been tested by the two + // previous tests tmax<0 || tmin>dmin (with dmin==0). } - // Note: for a segment the condition sign(tmin)*sign(tmax) > 0 has - // already been tested by the two previous tests tmax<0 || tmin>dmin - // (with dmin==0). // If dmin == 0, at this point, [t1, t2] == ]-inf, +inf[, or t1 or t2 // is a NaN. But the case with NaNs is treated as if the interval @@ -124,6 +123,8 @@ namespace internal { CGAL_assertion(dmin >= 0); CGAL_assertion(dmax >= 0); + // CGAL_assertion(!bounded_0 || ( (dmin == 0) == (px == qx && (px == bxmax || px == bxmin)) ) ); + // CGAL_assertion(!bounded_1 || ( (dmax == 0) == (px == qx && (px == bxmax || px == bxmin)) ) ); // ----------------------------------- // treat y coord @@ -188,6 +189,9 @@ namespace internal { CGAL_assertion(dymin >= 0); CGAL_assertion(dymax >= 0); + // CGAL_assertion(!bounded_0 || (dmin == 0) == (!bounded_0 && px == qx && py == qy)); + // CGAL_assertion(!bounded_1 || (dmax == 0) == (!bounded_1 && px == qx && py == qy)); + // ----------------------------------- // treat z coord @@ -238,7 +242,7 @@ namespace internal { } } - // If the querz is vertical for z, then check its z-coordinate is in + // If the query is vertical for z, then check its z-coordinate is in // the z-slab. if( (pz == qz) && // <=> (dmin == 0) (! (bounded_0 && bounded_1) ) ) // do not check for a segment @@ -253,21 +257,23 @@ namespace internal { CGAL_assertion(dzmin >= 0); CGAL_assertion(dzmax >= 0); - // If t1 > tymax || tymin > t2, return false. + // If t1 > tymax/dymax || tymin/dymin > t2, return false. if( dymax > 0 && dmin > 0 && (dmin*tymax) < (dymax*tmin) ) return false; if( dymin > 0 && dmax > 0 && (dmax*tymin) > (dymin*tmax) ) return false; - // If tymin > t1, set t1 = tymin. - if( dmin == 0 || - ( dymin > 0 && (dmin*tymin) > (dymin*tmin) ) ) + // If tymin/dymin > t1, set t1 = tymin/dymin. + if( (px == qx) || // <=> (dmin == 0) + ( (py != qy) && // <=> (dymin > 0) + (dmin*tymin) > (dymin*tmin) ) ) { tmin = tymin; dmin = dymin; } - // If tymax < t2, set t2 = tymax. - if( dmax == 0 || - ( dymax > 0 && (dmax*tymax) < (dymax*tmax) ) ) + // If tymax/dymax < t2, set t2 = tymax/dymax. + if( (px == qx) || // <=> (dmax > 0) + ( (py != qy) && // <=> dymax > 0 + (dmax*tymax) < (dymax*tmax) ) ) { tmax = tymax; dmax = dymax; @@ -275,14 +281,13 @@ namespace internal { CGAL_assertion(dmin >= 0); CGAL_assertion(dmax >= 0); + // CGAL_assertion((dmin == 0) == (px == qx && py == qy)); + // CGAL_assertion((dmax == 0) == (px == qx && py == qy)); // If t1 > tzmax || tzmin > t2, return false. if( dmin > 0 && dzmax > 0 && (dmin*tzmax) < (dzmax*tmin) ) return false; if( dmax > 0 && dzmin > 0 && (dmax*tzmin) > (dzmin*tmax) ) return false; - - return ( dzmin == 0 || - (!bounded_1 || tzmin <= dzmin) && - (!bounded_0 || tzmax >= FT(0)) ); // t1 <= 1 && t2 >= 0 + return true; } template Date: Wed, 21 Mar 2012 11:40:29 +0000 Subject: [PATCH 18/38] Less tests of sign of expressions Followup to previous commit. I have managed to transform most of the tests to simple comparison of input coordinates. That will ease the writing of static filters. Only six determinant signs have to be exactly determined. --- .../Bbox_3_Segment_3_do_intersect.h | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 3ce8b3efd61..70e4520ef53 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -258,13 +258,15 @@ namespace internal { CGAL_assertion(dzmax >= 0); // If t1 > tymax/dymax || tymin/dymin > t2, return false. - if( dymax > 0 && dmin > 0 && (dmin*tymax) < (dymax*tmin) ) return false; - if( dymin > 0 && dmax > 0 && (dmax*tymin) > (dymin*tmax) ) return false; + if( py != qy && px != qx) { // dmin > 0, dymax >0, dmax > 0, dymin > 0 + if( (dmin*tymax) < (dymax*tmin) ) return false; // TEST TO FILTER + if( (dmax*tymin) > (dymin*tmax) ) return false; // TEST TO FILTER + } // If tymin/dymin > t1, set t1 = tymin/dymin. if( (px == qx) || // <=> (dmin == 0) ( (py != qy) && // <=> (dymin > 0) - (dmin*tymin) > (dymin*tmin) ) ) + (dmin*tymin) > (dymin*tmin) ) ) // TEST TO FILTER { tmin = tymin; dmin = dymin; @@ -273,7 +275,7 @@ namespace internal { // If tymax/dymax < t2, set t2 = tymax/dymax. if( (px == qx) || // <=> (dmax > 0) ( (py != qy) && // <=> dymax > 0 - (dmax*tymax) < (dymax*tmax) ) ) + (dmax*tymax) < (dymax*tmax) ) ) // TEST TO FILTER { tmax = tymax; dmax = dymax; @@ -285,8 +287,13 @@ namespace internal { // CGAL_assertion((dmax == 0) == (px == qx && py == qy)); // If t1 > tzmax || tzmin > t2, return false. - if( dmin > 0 && dzmax > 0 && (dmin*tzmax) < (dzmax*tmin) ) return false; - if( dmax > 0 && dzmin > 0 && (dmax*tzmin) > (dzmin*tmax) ) return false; + if( (px != qx || + py != qy ) && + (pz != qz) ) // dmin > 0, dmax > 0, dzmax > 0, dzmin > 0 + { + if( (dmin*tzmax) < (dzmax*tmin) ) return false; // TEST TO FILTER + if( (dmax*tzmin) > (dzmin*tmax) ) return false; // TEST TO FILTER + } return true; } From 0e6dbc4b4051439e8ceec121b43bb09845356cbf Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 21 Mar 2012 11:49:30 +0000 Subject: [PATCH 19/38] Reorganize expressions that are evaluation of the sign of a determinant There are six of them in do_intersect(Bbox_3, Segment_3) --- .../Intersections_3/Bbox_3_Segment_3_do_intersect.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 70e4520ef53..14bdd904d02 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -259,14 +259,14 @@ namespace internal { // If t1 > tymax/dymax || tymin/dymin > t2, return false. if( py != qy && px != qx) { // dmin > 0, dymax >0, dmax > 0, dymin > 0 - if( (dmin*tymax) < (dymax*tmin) ) return false; // TEST TO FILTER - if( (dmax*tymin) > (dymin*tmax) ) return false; // TEST TO FILTER + if( (dymax* tmin) > ( dmin*tymax) ) return false; // TEST TO FILTER + if( ( dmax*tymin) > (dymin* tmax) ) return false; // TEST TO FILTER } // If tymin/dymin > t1, set t1 = tymin/dymin. if( (px == qx) || // <=> (dmin == 0) ( (py != qy) && // <=> (dymin > 0) - (dmin*tymin) > (dymin*tmin) ) ) // TEST TO FILTER + ( dmin*tymin) > (dymin* tmin) ) ) // TEST TO FILTER { tmin = tymin; dmin = dymin; @@ -275,7 +275,7 @@ namespace internal { // If tymax/dymax < t2, set t2 = tymax/dymax. if( (px == qx) || // <=> (dmax > 0) ( (py != qy) && // <=> dymax > 0 - (dmax*tymax) < (dymax*tmax) ) ) // TEST TO FILTER + (dymax* tmax) > ( dmax*tymax) ) ) // TEST TO FILTER { tmax = tymax; dmax = dymax; @@ -291,8 +291,8 @@ namespace internal { py != qy ) && (pz != qz) ) // dmin > 0, dmax > 0, dzmax > 0, dzmin > 0 { - if( (dmin*tzmax) < (dzmax*tmin) ) return false; // TEST TO FILTER - if( (dmax*tzmin) > (dzmin*tmax) ) return false; // TEST TO FILTER + if( (dzmax* tmin) > ( dmin*tzmax) ) return false; // TEST TO FILTER + if( ( dmax*tzmin) > (dzmin* tmax) ) return false; // TEST TO FILTER } return true; } From e940d97a7511586d946d61257c55b62007668e7f Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 21 Mar 2012 13:59:45 +0000 Subject: [PATCH 20/38] Prepare the code factorization with the static filter of Do_intersect_3 --- .../Bbox_3_Segment_3_do_intersect.h | 34 +++++++++++++++---- 1 file changed, 28 insertions(+), 6 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 14bdd904d02..ca5aae2e628 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -41,6 +41,16 @@ namespace CGAL { namespace internal { + template + struct Do_intersect_bbox_segment_aux_is_greater { + typedef bool result_type; + + template + bool operator()(const FT& a, const FT& b) const { + return a > b; + } + }; + template @@ -257,29 +267,39 @@ namespace internal { CGAL_assertion(dzmin >= 0); CGAL_assertion(dzmax >= 0); + typedef Do_intersect_bbox_segment_aux_is_greater Is_greater; + typedef typename Is_greater::result_type Is_greater_value; + Is_greater is_greater; + // If t1 > tymax/dymax || tymin/dymin > t2, return false. if( py != qy && px != qx) { // dmin > 0, dymax >0, dmax > 0, dymin > 0 - if( (dymax* tmin) > ( dmin*tymax) ) return false; // TEST TO FILTER - if( ( dmax*tymin) > (dymin* tmax) ) return false; // TEST TO FILTER + const Is_greater_value b1 = is_greater(dymax* tmin, dmin*tymax); + if(possibly(b1)) return !b1; // if(is_greater) return false; // or uncertain + if( is_greater(dymax* tmin, dmin*tymax) ) return false; + const Is_greater_value b2 = is_greater( dmax*tymin, dymin* tmax); + if(possibly(b2)) return !b2; } + Is_greater_value b; // If tymin/dymin > t1, set t1 = tymin/dymin. if( (px == qx) || // <=> (dmin == 0) ( (py != qy) && // <=> (dymin > 0) - ( dmin*tymin) > (dymin* tmin) ) ) // TEST TO FILTER + certainly(b = is_greater( dmin*tymin, dymin* tmin)) ) ) { tmin = tymin; dmin = dymin; } + if(is_indeterminate(b)) return b; // If tymax/dymax < t2, set t2 = tymax/dymax. if( (px == qx) || // <=> (dmax > 0) ( (py != qy) && // <=> dymax > 0 - (dymax* tmax) > ( dmax*tymax) ) ) // TEST TO FILTER + certainly(b = is_greater(dymax* tmax, dmax*tymax)) ) ) { tmax = tymax; dmax = dymax; } + if(is_indeterminate(b)) return b; CGAL_assertion(dmin >= 0); CGAL_assertion(dmax >= 0); @@ -291,8 +311,10 @@ namespace internal { py != qy ) && (pz != qz) ) // dmin > 0, dmax > 0, dzmax > 0, dzmin > 0 { - if( (dzmax* tmin) > ( dmin*tzmax) ) return false; // TEST TO FILTER - if( ( dmax*tzmin) > (dzmin* tmax) ) return false; // TEST TO FILTER + const Is_greater_value b1 = is_greater(dzmax* tmin, dmin*tzmax); + if(possibly(b1)) return !b1; // if(is_greater) return false; // or uncertain + const Is_greater_value b2 = is_greater( dmax*tzmin, dzmin* tmax); + if(possibly(b2)) return !b2; // if(is_greater) return false; // or uncertain } return true; } From b5d987a4b56ff2700140d9b2531940b08e2b5024 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 21 Mar 2012 14:24:31 +0000 Subject: [PATCH 21/38] Increasing the perf of the filtered predicate When FT is Interval_nt, it is better that the Do_intersect_bbox_segment_aux_is_greater returns a Uncertain instead of a bool. That can delay the conversion of Uncertain to bool, and hence better perf. --- .../Intersections_3/Bbox_3_Segment_3_do_intersect.h | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index ca5aae2e628..58d08b306b8 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -25,6 +25,7 @@ #include #include +#include // inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf @@ -41,12 +42,12 @@ namespace CGAL { namespace internal { - template - struct Do_intersect_bbox_segment_aux_is_greater { - typedef bool result_type; + template + struct Do_intersect_bbox_segment_aux_is_greater + { + typedef typename Same_uncertainty::type result_type; - template - bool operator()(const FT& a, const FT& b) const { + result_type operator()(const FT& a, const FT& b) const { return a > b; } }; @@ -267,7 +268,7 @@ namespace internal { CGAL_assertion(dzmin >= 0); CGAL_assertion(dzmax >= 0); - typedef Do_intersect_bbox_segment_aux_is_greater Is_greater; + typedef Do_intersect_bbox_segment_aux_is_greater Is_greater; typedef typename Is_greater::result_type Is_greater_value; Is_greater is_greater; From 439dd0abb96acbac4942790b1400b42087d393ed Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 21 Mar 2012 14:34:53 +0000 Subject: [PATCH 22/38] Fix typo - Remove a line that was rewritten but not removed. - Add a comment --- .../internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 58d08b306b8..01848a39f4e 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -276,7 +276,6 @@ namespace internal { if( py != qy && px != qx) { // dmin > 0, dymax >0, dmax > 0, dymin > 0 const Is_greater_value b1 = is_greater(dymax* tmin, dmin*tymax); if(possibly(b1)) return !b1; // if(is_greater) return false; // or uncertain - if( is_greater(dymax* tmin, dmin*tymax) ) return false; const Is_greater_value b2 = is_greater( dmax*tymin, dymin* tmax); if(possibly(b2)) return !b2; } @@ -290,7 +289,9 @@ namespace internal { tmin = tymin; dmin = dymin; } - if(is_indeterminate(b)) return b; + if(is_indeterminate(b)) return b; // Note that the default-constructed + // Is_greater_value cannot be + // indeterminate. // If tymax/dymax < t2, set t2 = tymax/dymax. if( (px == qx) || // <=> (dmax > 0) From 672d2dd2929ef7857931be4adb3d2cb744aaf52b Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 21 Mar 2012 16:58:29 +0000 Subject: [PATCH 23/38] Commit work in progress --- .../Bbox_3_Ray_3_do_intersect.h | 2 +- .../Bbox_3_Segment_3_do_intersect.h | 96 +++++++++++++++++-- .../bbox_other_do_intersect_test.cpp | 27 +++--- 3 files changed, 102 insertions(+), 23 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h index 6a402e1cd4d..582fec09c33 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h @@ -136,7 +136,7 @@ namespace internal { const Point_3& source = ray.source(); const Point_3& point_on_ray = ray.second_point(); - return do_intersect_bbox_segment_aux( + return do_intersect_bbox_segment_aux( source.x(), source.y(), source.z(), point_on_ray.x(), point_on_ray.y(), point_on_ray.z(), bbox.xmin(), bbox.ymin(), bbox.zmin(), diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 01848a39f4e..75d4e613117 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -1,5 +1,4 @@ -// Copyright (c) 2008 ETH Zurich (Switzerland) -// Copyright (c) 2008-2009 INRIA Sophia-Antipolis (France) +// Copyright (c) 2012 GeometryFactory Sarl (France) // All rights reserved. // // This file is part of CGAL (www.cgal.org); you can redistribute it and/or @@ -17,7 +16,7 @@ // $Id$ // // -// Author(s) : Camille Wormser, Jane Tournois, Pierre Alliez, Stephane Tayeb +// Author(s) : Laurent Rineau #ifndef CGAL_INTERNAL_INTERSECTIONS_3_BBOX_3_SEGMENT_3_DO_INTERSECT_H @@ -26,6 +25,8 @@ #include #include #include +#include +#include // inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf @@ -47,17 +48,77 @@ namespace internal { { typedef typename Same_uncertainty::type result_type; + void register_new_input_values(const FT& t, const FT& d) {} + void compute_new_error_bound() {} + bool bound_overflow() { return false; } + bool value_might_underflow() { return false; } + + static result_type uncertain() { + return true; + } + result_type operator()(const FT& a, const FT& b) const { return a > b; } - }; + }; // end struct template Do_intersect_bbox_segment_aux_is_greater + + template + class Do_intersect_bbox_segment_aux_is_greater + { + double error; + double tmax; + double dmax; + + static const double EPS = 8.8872057372592798e-16; + static const double OVERF = 1e153; + static const double UNDERF = 1e-146; + + public: + CGAL_static_assertion((boost::is_same::value)); + + Do_intersect_bbox_segment_aux_is_greater() : error(0.), tmax(0.), dmax(0.) {} + + void register_new_input_values(const double& t, const double& d) { + const double at = CGAL::abs(t); + const double ad = CGAL::abs(d); + if(at > tmax) tmax = at; + if(ad > dmax) dmax = ad; + } + + void compute_new_error_bound() { + error = tmax * dmax * EPS; + } + + bool bound_overflow() { + return dmax > OVERF && tmax > OVERF; + } + + bool value_might_underflow() { + return dmax < UNDERF || dmax < UNDERF; + } + + typedef Uncertain result_type; + + static result_type uncertain() { + return result_type::indeterminate(); + } + + result_type operator()(const FT& a, const FT& b) const { + const FT x = a - b; + if(x > error) return true; + else if(x < -error) return false; + else return uncertain(); + } + + }; // end specialization Do_intersect_bbox_segment_aux_is_greater template + bool bounded_1, + bool use_static_filters> // __attribute__ ((noinline)) inline - bool + typename Do_intersect_bbox_segment_aux_is_greater::result_type do_intersect_bbox_segment_aux( const FT& px, const FT& py, const FT& pz, const FT& qx, const FT& qy, const FT& qz, @@ -268,12 +329,22 @@ namespace internal { CGAL_assertion(dzmin >= 0); CGAL_assertion(dzmax >= 0); - typedef Do_intersect_bbox_segment_aux_is_greater Is_greater; + typedef Do_intersect_bbox_segment_aux_is_greater Is_greater; typedef typename Is_greater::result_type Is_greater_value; Is_greater is_greater; + is_greater.register_new_input_values(tmin, dmin); + is_greater.register_new_input_values(tymin, dymin); + is_greater.register_new_input_values(tmax, dmax); + is_greater.register_new_input_values(tymax, dymax); + + is_greater.compute_new_error_bound(); + if(is_greater.bound_overflow() || is_greater.value_might_underflow()) + return Is_greater::uncertain(); + // If t1 > tymax/dymax || tymin/dymin > t2, return false. - if( py != qy && px != qx) { // dmin > 0, dymax >0, dmax > 0, dymin > 0 + if( py != qy && px != qx ) { // dmin > 0, dymax >0, dmax > 0, dymin > 0 const Is_greater_value b1 = is_greater(dymax* tmin, dmin*tymax); if(possibly(b1)) return !b1; // if(is_greater) return false; // or uncertain const Is_greater_value b2 = is_greater( dmax*tymin, dymin* tmax); @@ -313,6 +384,13 @@ namespace internal { py != qy ) && (pz != qz) ) // dmin > 0, dmax > 0, dzmax > 0, dzmin > 0 { + is_greater.register_new_input_values(tzmin, dzmin); + is_greater.register_new_input_values(tzmax, dzmax); + + is_greater.compute_new_error_bound(); + if(is_greater.bound_overflow() || is_greater.value_might_underflow()) + return Is_greater::uncertain(); + const Is_greater_value b1 = is_greater(dzmax* tmin, dmin*tzmax); if(possibly(b1)) return !b1; // if(is_greater) return false; // or uncertain const Is_greater_value b2 = is_greater( dmax*tzmin, dzmin* tmax); @@ -366,7 +444,7 @@ namespace internal { const Point_3& source = segment.source(); const Point_3& target = segment.target(); - return do_intersect_bbox_segment_aux( + return do_intersect_bbox_segment_aux( source.x(), source.y(), source.z(), target.x(), target.y(), target.z(), bbox.xmin(), bbox.ymin(), bbox.zmin(), diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index 07ecc7de679..c32d01aed04 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -673,24 +673,25 @@ int main() srand(0); std::cout << std::setprecision(5); - std::cout << "Testing with Simple_cartesian..." << std::endl ; - bool b = test_kernel >(false); + bool b; + // std::cout << "Testing with Simple_cartesian..." << std::endl ; + // b = test_kernel >(false); - std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; - b &= test_kernel >(true); + // std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; + // b &= test_kernel >(true); - std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; - b &= test_kernel >(true); + // std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; + // b &= test_kernel >(true); - std::cout << std::endl << "Testing with Cartesian..." << std::endl ; - b &= test_kernel >(false); + // std::cout << std::endl << "Testing with Cartesian..." << std::endl ; + // b &= test_kernel >(false); - std::cout << std::endl << "Testing with Cartesian..." << std::endl ; - b &= test_kernel >(true); + // std::cout << std::endl << "Testing with Cartesian..." << std::endl ; + // b &= test_kernel >(true); - std::cout << std::endl << "Testing with Filtered_kernel > without static filters..." << std::endl ; - typedef CGAL::Filtered_kernel, false> Fk_no_static; - b &= test_kernel(); + // std::cout << std::endl << "Testing with Filtered_kernel > without static filters..." << std::endl ; + // typedef CGAL::Filtered_kernel, false> Fk_no_static; + // b &= test_kernel(); std::cout << std::endl << "Testing with Exact_predicates_inexact_constructions_kernel..." << std::endl ; b &= test_kernel(); From 00621279fb7da917910d371c0aa6b5b8c9166656 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 23 Mar 2012 13:02:46 +0000 Subject: [PATCH 24/38] Commit the new version of the static filter. Too slow for the moment. --- .../CGAL/internal/Static_filters/Do_intersect_3.h | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index 544d51929f9..34c378f4a1d 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include @@ -121,6 +122,16 @@ public: { CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + const Uncertain ub = + do_intersect_bbox_segment_aux + (px, py, pz, + qx, qy, qz, + bxmin, bymin, bzmin, + bxmax, bymax, bzmax); + + if(is_indeterminate(ub)) return Base::operator()(s,b); + else return ub.sup(); + // ----------------------------------- // treat x coord // ----------------------------------- From a72bd80380f5ee91e558b054669c0e458147aefd Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 14 Jun 2012 17:04:35 +0000 Subject: [PATCH 25/38] Pass the bbox as argument instead of the six coordinates That increased the perfs! :-) --- .../internal/Static_filters/Do_intersect_3.h | 3 +- .../Bbox_3_Ray_3_do_intersect.h | 3 +- .../Bbox_3_Segment_3_do_intersect.h | 15 ++++++--- .../bbox_other_do_intersect_test.cpp | 31 +++++++++---------- 4 files changed, 27 insertions(+), 25 deletions(-) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index 34c378f4a1d..7e60d752c17 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -126,8 +126,7 @@ public: do_intersect_bbox_segment_aux (px, py, pz, qx, qy, qz, - bxmin, bymin, bzmin, - bxmax, bymax, bzmax); + b); if(is_indeterminate(ub)) return Base::operator()(s,b); else return ub.sup(); diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h index 582fec09c33..25e56bf062d 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h @@ -139,8 +139,7 @@ namespace internal { return do_intersect_bbox_segment_aux( source.x(), source.y(), source.z(), point_on_ray.x(), point_on_ray.y(), point_on_ray.z(), - bbox.xmin(), bbox.ymin(), bbox.zmin(), - bbox.xmax(), bbox.ymax(), bbox.zmax() ); + bbox); // const CGAL::cpp0x::array rray = {source.x(), source.y(), source.z(), // point_on_ray.x(), point_on_ray.y(), point_on_ray.z() }; // return do_intersect_bbox_segment_aux diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 75d4e613117..c84a850dbfb 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -122,9 +122,15 @@ namespace internal { do_intersect_bbox_segment_aux( const FT& px, const FT& py, const FT& pz, const FT& qx, const FT& qy, const FT& qz, - const double& bxmin, const double& bymin, const double& bzmin, - const double& bxmax, const double& bymax, const double& bzmax) + const Bbox_3& bbox) { + const double& bxmin = bbox.xmin(); + const double& bymin = bbox.ymin(); + const double& bzmin = bbox.zmin(); + const double& bxmax = bbox.xmax(); + const double& bymax = bbox.ymax(); + const double& bzmax = bbox.zmax(); + // The following code encode t1 and t2 by: // t1 = tmin/dmin // t2 = tmax/dmax @@ -447,8 +453,9 @@ namespace internal { return do_intersect_bbox_segment_aux( source.x(), source.y(), source.z(), target.x(), target.y(), target.z(), - bbox.xmin(), bbox.ymin(), bbox.zmin(), - bbox.xmax(), bbox.ymax(), bbox.zmax() ); + bbox); + // bbox.xmin(), bbox.ymin(), bbox.zmin(), + // bbox.xmax(), bbox.ymax(), bbox.zmax() ); // const CGAL::cpp0x::array seg = {source.x(), source.y(), source.z(), // target.x(), target.y(), target.z() }; diff --git a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp index c0877fdfe5e..3a0586e41e9 100644 --- a/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp +++ b/Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp @@ -67,7 +67,7 @@ bool test_aux(const T& t, bool b = CGAL::do_intersect(t,bbox); if ( b != expected ) - std::cout << "ERROR: do_intersect(" << name + std::cerr << "ERROR: do_intersect(" << name << ") did not answer the expected result !" << std::endl; return (b == expected); @@ -674,27 +674,24 @@ int main() std::cout << std::setprecision(5); bool b; - // std::cout << "Testing with Simple_cartesian..." << std::endl ; - // b = test_kernel >(false); + std::cout << "Testing with Simple_cartesian..." << std::endl ; + b = test_kernel >(false); - // std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; - // b &= test_kernel >(true); + std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; + b &= test_kernel >(true); - // std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; - // b &= test >(true); - // test_speed >(); + std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; + b &= test_kernel >(true); - // std::cout << std::endl << "Testing with Cartesian..." << std::endl ; - // b &= test >(false); - // test_speed >(); + std::cout << std::endl << "Testing with Cartesian..." << std::endl ; + b &= test_kernel >(false); - // std::cout << std::endl << "Testing with Cartesian..." << std::endl ; - // b &= test_kernel >(true); + std::cout << std::endl << "Testing with Cartesian..." << std::endl ; + b &= test_kernel >(true); - // std::cout << std::endl << "Testing with Filtered_kernel > without static filters..." << std::endl ; - // typedef CGAL::Filtered_kernel, false> Fk_no_static; - // b &= test_kernel(); - // test_speed(); + std::cout << std::endl << "Testing with Filtered_kernel > without static filters..." << std::endl ; + typedef CGAL::Filtered_kernel, false> Fk_no_static; + b &= test_kernel(); std::cout << std::endl << "Testing with Exact_predicates_inexact_constructions_kernel..." << std::endl ; b &= test_kernel(); From 1c6ba1851b47e322fc87855cb1f3968d2c41ee12 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 15 Jun 2012 14:57:36 +0000 Subject: [PATCH 26/38] Factorize code in CGAL/internal/Static_filters/Do_intersect_3.h The static filters of Do_intersect(BBox_3, Segment_3) and Do_intersect(BBox_3, Segment_3) now use do_intersect_bbox_segment_aux(..) from CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h --- .../internal/Static_filters/Do_intersect_3.h | 328 +----------------- 1 file changed, 7 insertions(+), 321 deletions(-) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index 7e60d752c17..60236923265 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -130,179 +130,6 @@ public: if(is_indeterminate(ub)) return Base::operator()(s,b); else return ub.sup(); - - // ----------------------------------- - // treat x coord - // ----------------------------------- - double dmin, dmax, tmin, tmax; - if ( qx >= px ) // this is input and needs no epsilon - { - if(px > bxmax) return false; // segment on the right of bbox - if(qx < bxmin) return false; // segment on the left of bbox - - tmax = bxmax - px; - - dmax = qx - px; - if ( bxmin < px ) // tmin < 0 means px is in the x-range of bbox - { - tmin = 0; - dmin = 1; - } else { - tmin = bxmin - px; - dmin = dmax; - } - } - else - { - if(qx > bxmax) return false; // segment on the right of bbox - if(px < bxmin) return false; // segment on the left of bbox - - tmax = px - bxmin; - - dmax = px - qx; - if ( px < bxmax ) // tmin < 0 means px is in the x-range of bbox - { - tmin = 0; - dmin = 1; - } else { - tmin = px - bxmax; - dmin = dmax; - } - } - - double m = CGAL::abs(tmin), m2; - m2 = CGAL::abs(tmax); if(m2 > m) { m = m2; } - m2 = CGAL::abs(dmin); if(m2 > m) { m = m2; } - - if(m < 7e-294) { - // underflow in the computation of 'error' - return Base::operator()(s,b); - } - - const double EPS_1 = 3.55618e-15; - - double error = EPS_1 * m; - - switch(sign_with_error( tmax - dmax, error)) { - case POSITIVE: - tmax = 1; - dmax = 1; - break; - case NEGATIVE: - break; - default: - // ambiguity of the comparison tmax > dmin - // let's call the exact predicate - return Base::operator()(s,b); - } - - // ----------------------------------- - // treat y coord - // ----------------------------------- - double d_, tmin_, tmax_; - if ( qy >= py ) // this is input and needs no epsilon - { - tmin_ = bymin - py; - tmax_ = bymax - py; - d_ = qy - py; - } - else - { - tmin_ = py - bymax; - tmax_ = py - bymin; - d_ = py - qy; - } - - m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; } - m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; } - m2 = CGAL::abs(d_); if(m2 > m) { m = m2; } - - if(m < 3e-147) { - // underflow in the computation of 'error' - return Base::operator()(s,b); - } - - error = EPS_1 * m * m; - - if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */ - // potential overflow on the computation of 'sign1' and 'sign2' - return Base::operator()(s,b); - } - Sign sign1 = sign_with_error( (d_*tmin) - (dmin*tmax_) , error); - Sign sign2 = sign_with_error( (dmax*tmin_) - (d_*tmax) , error); - - if(sign1 == POSITIVE || sign2 == POSITIVE) - return false; // We are *sure* the segment is outside the box, on one - // side or the other. - if(sign1 == ZERO || sign2 == ZERO) { - return Base::operator()(s,b); // We are *unsure*: one *may be* - // positive. - } - - // Here we are sure the two signs are negative. We can continue with - // the rest of the function... - - // epsilons needed - switch(sign_with_error((dmin*tmin_) - (d_*tmin) , error)) { - case POSITIVE: - tmin = tmin_; - dmin = d_; - break; - case NEGATIVE: - break; - default: // uncertainty - return Base::operator()(s,b); - } - - // epsilons needed - switch(sign_with_error((d_*tmax) - (dmax*tmax_) , error)) { - case POSITIVE: - tmax = tmax_; - dmax = d_; - case NEGATIVE: - break; - default: // uncertainty - return Base::operator()(s,b); - } - - // ----------------------------------- - // treat z coord - // ----------------------------------- - if ( qz >= pz ) // this is input and needs no epsilon - { - tmin_ = bzmin - pz; - tmax_ = bzmax - pz; - d_ = qz - pz; - } - else - { - tmin_ = pz - bzmax; - tmax_ = pz - bzmin; - d_ = pz - qz; - } - - m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; } - m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; } - m2 = CGAL::abs(d_); if(m2 > m) { m = m2; } - - // m may have changed - error = EPS_1 * m * m; - - if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */ - // potential overflow on the computation of 'sign1' and 'sign2' - return Base::operator()(s,b); - } - sign1 = sign_with_error( (dmin*tmax_) - (d_*tmin) , error); - sign2 = sign_with_error( (d_*tmax) - (dmax*tmin_) , error); - if(sign1 == NEGATIVE || sign2 == NEGATIVE) { - return false; // We are *sure* the segment is outside the box, on one - // side or the other. - } - if(sign1 == ZERO || sign2 == ZERO) { - return Base::operator()(s,b); // We are *unsure*: one *may be* - // negative. - } - return true; // We are *sure* the two signs are positive. } return Base::operator()(s,b); } @@ -337,155 +164,14 @@ public: { CGAL_BRANCH_PROFILER_BRANCH_1(tmp); - // ----------------------------------- - // treat x coord - // ----------------------------------- - double dmin, dmax, tmin, tmax; - if ( qx >= px ) // this is input and needs no epsilon - { - if(px > bxmax) return false; // if tmax < 0, ray on the right of bbox - tmax = bxmax - px; - - dmax = qx - px; - if ( bxmin < px ) // tmin < 0 means px is in the x-range of bbox - { - tmin = 0; - dmin = 1; - } else { - if( px == qx ) return false; // if dmin == 0 - tmin = bxmin - px; - dmin = dmax; - } - } - else - { - if(px < bxmin) return false; // if tmax < 0, ray on the left of bbox - tmax = px - bxmin; + const Uncertain ub = + do_intersect_bbox_segment_aux + (px, py, pz, + qx, qy, qz, + b); - dmax = px - qx; - if ( px < bxmax ) // tmin < 0 means px is in the x-range of bbox - { - tmin = 0; - dmin = 1; - } else { - if( px == qx ) return false; // if dmin == 0 - tmin = px - bxmax; - dmin = dmax; - } - } - - // ----------------------------------- - // treat y coord - // ----------------------------------- - double d_, tmin_, tmax_; - if ( qy >= py ) // this is input and needs no epsilon - { - tmin_ = bymin - py; - tmax_ = bymax - py; - d_ = qy - py; - } - else - { - tmin_ = py - bymax; - tmax_ = py - bymin; - d_ = py - qy; - } - - double m = CGAL::abs(tmin), m2; - m2 = CGAL::abs(tmax); if(m2 > m) { m = m2; } - m2 = CGAL::abs(dmin); if(m2 > m) { m = m2; } - m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; } - m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; } - m2 = CGAL::abs(d_); if(m2 > m) { m = m2; } - - if(m < 3e-147) { - // underflow in the computation of 'error' - return Base::operator()(r,b); - } - - const double EPS_1 = 3.55618e-15; - - double error = EPS_1 * m * m; - - if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */ - // potential overflow on the computation of 'sign1' and 'sign2' - return Base::operator()(r,b); - } - Sign sign1 = sign_with_error( (d_*tmin) - (dmin*tmax_) , error); - Sign sign2 = sign_with_error( (dmax*tmin_) - (d_*tmax) , error); - - if(sign1 == POSITIVE || sign2 == POSITIVE) - return false; // We are *sure* the ray is outside the box, on one - // side or the other. - if(sign1 == ZERO || sign2 == ZERO) { - return Base::operator()(r,b); // We are *unsure*: one *may be* - // positive. - } - - // Here we are sure the two signs are negative. We can continue with - // the rest of the function... - - // epsilons needed - switch(sign_with_error((dmin*tmin_) - (d_*tmin) , error)) { - case POSITIVE: - tmin = tmin_; - dmin = d_; - break; - case NEGATIVE: - break; - default: // uncertainty - return Base::operator()(r,b); - } - - // epsilons needed - switch(sign_with_error((d_*tmax) - (dmax*tmax_) , error)) { - case POSITIVE: - tmax = tmax_; - dmax = d_; - case NEGATIVE: - break; - default: // uncertainty - return Base::operator()(r,b); - } - - // ----------------------------------- - // treat z coord - // ----------------------------------- - if ( qz >= pz ) // this is input and needs no epsilon - { - tmin_ = bzmin - pz; - tmax_ = bzmax - pz; - d_ = qz - pz; - } - else - { - tmin_ = pz - bzmax; - tmax_ = pz - bzmin; - d_ = pz - qz; - } - - m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; } - m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; } - m2 = CGAL::abs(d_); if(m2 > m) { m = m2; } - - // m may have changed - error = EPS_1 * m * m; - - if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */ - // potential overflow on the computation of 'sign1' and 'sign2' - return Base::operator()(r,b); - } - sign1 = sign_with_error( (dmin*tmax_) - (d_*tmin) , error); - sign2 = sign_with_error( (d_*tmax) - (dmax*tmin_) , error); - if(sign1 == NEGATIVE || sign2 == NEGATIVE) { - return false; // We are *sure* the ray is outside the box, on one - // side or the other. - } - if(sign1 == ZERO || sign2 == ZERO) { - return Base::operator()(r,b); // We are *unsure*: one *may be* - // negative. - } - return true; // We are *sure* the two signs are positive. + if(is_indeterminate(ub)) return Base::operator()(r,b); + else return ub.sup(); } return Base::operator()(r,b); } From 41ba29e19abf2c4fd9cb90d916c805532059d0a4 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 15 Jun 2012 15:06:34 +0000 Subject: [PATCH 27/38] Factorize code of Do_intersect(Bbox_3, ) Uniform use of do_intersect_bbox_segment_aux(..) with various Boolean template arguments. --- .../internal/Static_filters/Do_intersect_3.h | 16 ++- .../Bbox_3_Line_3_do_intersect.h | 113 ++--------------- .../Bbox_3_Ray_3_do_intersect.h | 114 ++---------------- 3 files changed, 38 insertions(+), 205 deletions(-) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index 60236923265..1961ad0b438 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -28,8 +28,10 @@ #include #include #include + #include -#include +// for CGAL::internal::do_intersect_bbox_segment_aux + #include // inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf @@ -123,7 +125,11 @@ public: CGAL_BRANCH_PROFILER_BRANCH_1(tmp); const Uncertain ub = - do_intersect_bbox_segment_aux + do_intersect_bbox_segment_aux + // do use static filters (px, py, pz, qx, qy, qz, b); @@ -165,7 +171,11 @@ public: CGAL_BRANCH_PROFILER_BRANCH_1(tmp); const Uncertain ub = - do_intersect_bbox_segment_aux + do_intersect_bbox_segment_aux + // do use static filters (px, py, pz, qx, qy, qz, b); diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Line_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Line_3_do_intersect.h index 6175e1c37ef..282bbca3850 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Line_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Line_3_do_intersect.h @@ -26,107 +26,15 @@ #include #include +#include +// for CGAL::internal::do_intersect_bbox_segment_aux + // inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf namespace CGAL { namespace internal { - template - inline - bool - bbox_line_do_intersect_aux(const FT& px, const FT& py, const FT& pz, - const FT& vx, const FT& vy, const FT& vz, - const FT& bxmin, const FT& bymin, const FT& bzmin, - const FT& bxmax, const FT& bymax, const FT& bzmax) - { - // ----------------------------------- - // treat x coord - // ----------------------------------- - FT dmin, tmin, tmax; - if ( vx >= 0 ) - { - tmin = bxmin - px; - tmax = bxmax - px; - dmin = vx; - } - else - { - tmin = px - bxmax; - tmax = px - bxmin; - dmin = -vx; - } - - //if px is not in the x-slab - if ( dmin == FT(0) && (tmin > FT(0) || tmax < FT(0)) ) return false; - - FT dmax = dmin; - - // ----------------------------------- - // treat y coord - // ----------------------------------- - FT d_, tmin_, tmax_; - if ( vy >= 0 ) - { - tmin_ = bymin - py; - tmax_ = bymax - py; - d_ = vy; - } - else - { - tmin_ = py - bymax; - tmax_ = py - bymin; - d_ = -vy; - } - - - - - if ( d_ == FT(0) ){ - //if py is not in the y-slab - if( (tmin_ > FT(0) || tmax_ < FT(0)) ) return false; - } - else - if ( (dmin*tmax_) < (d_*tmin) || (dmax*tmin_) > (d_*tmax) ) - return false; - - if( (dmin*tmin_) > (d_*tmin) ) - { - tmin = tmin_; - dmin = d_; - } - - if( (dmax*tmax_) < (d_*tmax) ) - { - tmax = tmax_; - dmax = d_; - } - - // ----------------------------------- - // treat z coord - // ----------------------------------- - if ( vz >= 0 ) - { - tmin_ = bzmin - pz; - tmax_ = bzmax - pz; - d_ = vz; - } - else - { - tmin_ = pz - bzmax; - tmax_ = pz - bzmin; - d_ = -vz; - } - - //if pz is not in the z-slab - //if ( d_ == FT(0) && (tmin_ > FT(0) || tmax_ < FT(0)) ) return false; - //The previous line is not needed as either dmin or d_ are not 0 - //(otherwise the direction of the line would be null). - // The following is equivalent to the in z-slab test if d_=0. - - return ( (dmin*tmax_) >= (d_*tmin) && (dmax*tmin_) <= (d_*tmax) ); - } - template bool do_intersect(const typename K::Line_3& line, const CGAL::Bbox_3& bbox, @@ -139,11 +47,16 @@ namespace internal { const Point_3& point = line.point(); const Vector_3& v = line.to_vector(); - return bbox_line_do_intersect_aux( - point.x(), point.y(), point.z(), - v.x(), v.y(), v.z(), - FT(bbox.xmin()), FT(bbox.ymin()), FT(bbox.zmin()), - FT(bbox.xmax()), FT(bbox.ymax()), FT(bbox.zmax()) ); + return do_intersect_bbox_segment_aux + // do not use static filters + ( + point.x(), point.y(), point.z(), + v.x(), v.y(), v.z(), + bbox + ); } } // namespace internal diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h index 25e56bf062d..249264aca74 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h @@ -24,7 +24,9 @@ #include #include + #include +// for CGAL::internal::do_intersect_bbox_segment_aux // inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf @@ -32,99 +34,6 @@ namespace CGAL { namespace internal { - template - inline - bool - bbox_ray_do_intersect_aux(const FT& px, const FT& py, const FT& pz, - const FT& qx, const FT& qy, const FT& qz, - const FT& bxmin, const FT& bymin, const FT& bzmin, - const FT& bxmax, const FT& bymax, const FT& bzmax) - { - // (px, py, pz) is the source - // ----------------------------------- - // treat x coord - // ----------------------------------- - FT dmin, tmin, tmax; - if ( qx >= px ) - { - tmin = bxmin - px; - tmax = bxmax - px; - dmin = qx - px; - if ( tmax < FT(0) ) - return false; - } - else - { - tmin = px - bxmax; - tmax = px - bxmin; - dmin = px - qx; - if ( tmax < FT(0) ) - return false; - } - - FT dmax = dmin; - if ( tmin > FT(0) ) - { - if ( dmin == FT(0) ) - return false; - } - else - { - tmin = FT(0); - dmin = FT(1); - } - - // ----------------------------------- - // treat y coord - // ----------------------------------- - FT d_, tmin_, tmax_; - if ( qy >= py ) - { - tmin_ = bymin - py; - tmax_ = bymax - py; - d_ = qy - py; - } - else - { - tmin_ = py - bymax; - tmax_ = py - bymin; - d_ = py - qy; - } - - if ( (dmin*tmax_) < (d_*tmin) || (dmax*tmin_) > (d_*tmax) ) - return false; - - if( (dmin*tmin_) > (d_*tmin) ) - { - tmin = tmin_; - dmin = d_; - } - - if( (dmax*tmax_) < (d_*tmax) ) - { - tmax = tmax_; - dmax = d_; - } - - // ----------------------------------- - // treat z coord - // ----------------------------------- - if ( qz >= pz ) - { - tmin_ = bzmin - pz; - tmax_ = bzmax - pz; - d_ = qz - pz; - } - else - { - tmin_ = pz - bzmax; - tmax_ = pz - bzmin; - d_ = pz - qz; - } - - return ( (dmin*tmax_) >= (d_*tmin) && (dmax*tmin_) <= (d_*tmax) ); - } - template bool do_intersect(const typename K::Ray_3& ray, const CGAL::Bbox_3& bbox, @@ -136,15 +45,16 @@ namespace internal { const Point_3& source = ray.source(); const Point_3& point_on_ray = ray.second_point(); - return do_intersect_bbox_segment_aux( - source.x(), source.y(), source.z(), - point_on_ray.x(), point_on_ray.y(), point_on_ray.z(), - bbox); - // const CGAL::cpp0x::array rray = {source.x(), source.y(), source.z(), - // point_on_ray.x(), point_on_ray.y(), point_on_ray.z() }; - // return do_intersect_bbox_segment_aux - // (rray, - // *reinterpret_cast*>(&*bbox.cartesian_begin()) ); + return do_intersect_bbox_segment_aux + // do not use static filters + ( + source.x(), source.y(), source.z(), + point_on_ray.x(), point_on_ray.y(), point_on_ray.z(), + bbox + ); } } // namespace internal From a3587eabf326e69033e26a78249c17dfb55364f2 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 15 Jun 2012 15:12:04 +0000 Subject: [PATCH 28/38] A version with less branches and more numerical computation --- .../Bbox_3_Segment_3_do_intersect.h | 166 ++++++------------ 1 file changed, 55 insertions(+), 111 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index c84a850dbfb..d00814f912c 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -139,50 +139,31 @@ namespace internal { // ----------------------------------- // treat x coord // ----------------------------------- - FT dmin, tmin, tmax, dmax; + FT dmin, tmin, tmax, dmax, tmin_minus_dmin, tmax_minus_dmax; if ( qx >= px ) { - if(bounded_0 && px > bxmax) return false; // segment on the right of bbox - if(bounded_1 && qx < bxmin) return false; // segment on the left of bbox - - if(bounded_1 && bxmax > qx) { - tmax = 1; - dmax = 1; - } else { - tmax = bxmax - px; - dmax = qx - px; - } - - if(bounded_0 && bxmin < px) // tmin < 0 means px is in the x-range of bbox - { - tmin = 0; - dmin = 1; - } else { - tmin = bxmin - px; - dmin = qx - px; - } + tmin = bxmin - px; + tmax = bxmax - px; + dmax = dmin = qx - px; + tmin_minus_dmin = bxmin - qx; + tmax_minus_dmax = bxmax - qx; } else { - if(bounded_1 && qx > bxmax) return false; // segment on the right of bbox - if(bounded_0 && px < bxmin) return false; // segment on the left of bbox + tmin = px - bxmax; + tmax = px - bxmin; + dmax = dmin = px - qx; + tmin_minus_dmin = qx - bxmax; + tmax_minus_dmax = qx - bxmin; + } - if(bounded_1 && bxmin < qx) { - tmax = 1; - dmax = 1; - } else { - tmax = px - bxmin; - dmax = px - qx; - } + if( bounded_0 && tmax < 0. ) return false; // t2 < 0 + if( bounded_1 && tmin_minus_dmin > 0. ) return false; // t1 > 1 - if(bounded_0 && px < bxmax) // tmin < 0 means px is in the x-range of bbox - { - tmin = 0; - dmin = 1; - } else { - tmin = px - bxmax; - dmin = px - qx; - } + if(bounded_0) tmin = CGAL::max(tmin, FT(0)); // t1 = max(t1,0) + if(bounded_1 && tmax_minus_dmax > 0) { // t2 = min(t2,1); + tmax = FT(1); + dmax = FT(1); } // If the query is vertical for x, then check its x-coordinate is in @@ -207,50 +188,31 @@ namespace internal { // ----------------------------------- // treat y coord // ----------------------------------- - FT dymin, tymin, tymax, dymax; + FT dymin, tymin, tymax, dymax, tymin_minus_dmin, tymax_minus_dmax; if ( qy >= py ) { - if(bounded_0 && py > bymax) return false; // segment on the right of bbox - if(bounded_1 && qy < bymin) return false; // segment on the left of bbox - - if(bounded_1 && bymax > qy) { - tymax = 1; - dymax = 1; - } else { - tymax = bymax - py; - dymax = qy - py; - } - - if(bounded_0 && bymin < py) // tmin < 0 means py is in the y-range of bbox - { - tymin = 0; - dymin = 1; - } else { - tymin = bymin - py; - dymin = qy - py; - } + tymin = bymin - py; + tymax = bymax - py; + dymax = dymin = qy - py; + tymin_minus_dmin = bymin - qy; + tymax_minus_dmax = bymax - qy; } else { - if(bounded_1 && qy > bymax) return false; // segment on the right of bbox - if(bounded_0 && py < bymin) return false; // segment on the left of bbox + tymin = py - bymax; + tymax = py - bymin; + dymax = dymin = py - qy; + tymin_minus_dmin = qy - bymax; + tymax_minus_dmax = qy - bymin; + } - if(bounded_1 && bymin < qy) { - tymax = 1; - dymax = 1; - } else { - tymax = py - bymin; - dymax = py - qy; - } + if( bounded_0 && tymax < 0. ) return false; // t2 < 0 + if( bounded_1 && tymin_minus_dmin > 0. ) return false; // t1 > 1 - if(bounded_0 && py < bymax) // tmin < 0 means py is in the y-range of bbox - { - tymin = 0; - dymin = 1; - } else { - tymin = py - bymax; - dymin = py - qy; - } + if(bounded_0) tymin = CGAL::max(tymin, FT(0)); // t1 = max(t1,0) + if(bounded_1 && tymax_minus_dmax > 0) { // t2 = min(t2,1); + tymax = FT(1); + dymax = FT(1); } // If the query is vertical for y, then check its y-coordinate is in @@ -274,50 +236,32 @@ namespace internal { // ----------------------------------- // treat z coord // ----------------------------------- - FT dzmin, tzmin, tzmax, dzmax; + FT dzmin, tzmin, tzmax, dzmax, tzmax_minus_dmax, tzmin_minus_dmin; + if ( qz >= pz ) { - if(bounded_0 && pz > bzmax) return false; // segment on the right of bbox - if(bounded_1 && qz < bzmin) return false; // segment on the left of bbox - - if(bounded_1 && bzmax > qz) { - tzmax = 1; - dzmax = 1; - } else { - tzmax = bzmax - pz; - dzmax = qz - pz; - } - - if(bounded_0 && bzmin < pz) // tmin < 0 means pz is in the z-range of bbox - { - tzmin = 0; - dzmin = 1; - } else { - tzmin = bzmin - pz; - dzmin = qz - pz; - } + tzmin = bzmin - pz; + tzmax = bzmax - pz; + dzmax = dzmin = qz - pz; + tzmin_minus_dmin = bzmin - qz; + tzmax_minus_dmax = bzmax - qz; } else { - if(bounded_1 && qz > bzmax) return false; // segment on the right of bbox - if(bounded_0 && pz < bzmin) return false; // segment on the left of bbox + tzmin = pz - bzmax; + tzmax = pz - bzmin; + dzmax = dzmin = pz - qz; + tzmin_minus_dmin = qz - bzmax; + tzmax_minus_dmax = qz - bzmin; + } - if(bounded_1 && bzmin < qz) { - tzmax = 1; - dzmax = 1; - } else { - tzmax = pz - bzmin; - dzmax = pz - qz; - } + if( bounded_0 && tzmax < 0. ) return false; // t2 < 0 + if( bounded_1 && tzmin_minus_dmin > 0. ) return false; // t1 > 1 - if(bounded_0 && pz < bzmax) // tmin < 0 means pz is in the z-range of bbox - { - tzmin = 0; - dzmin = 1; - } else { - tzmin = pz - bzmax; - dzmin = pz - qz; - } + if(bounded_0) tzmin = CGAL::max(tzmin, FT(0)); // t1 = max(t1,0) + if(bounded_1 && tzmax_minus_dmax > 0) { // t2 = min(t2,1); + tzmax = FT(1); + dzmax = FT(1); } // If the query is vertical for z, then check its z-coordinate is in From 02754f17f9455e32916d362cf188996c49c4ed1d Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 15 Jun 2012 15:27:17 +0000 Subject: [PATCH 29/38] Revert previous commit for CGAL/internal/Intersections_3/Bbox_3_Line_3_do_intersect.h The code cannot be easily reused because a line does not store two points but the a, b, c, d coefficients. --- .../Bbox_3_Line_3_do_intersect.h | 113 ++++++++++++++++-- 1 file changed, 100 insertions(+), 13 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Line_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Line_3_do_intersect.h index 282bbca3850..6175e1c37ef 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Line_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Line_3_do_intersect.h @@ -26,15 +26,107 @@ #include #include -#include -// for CGAL::internal::do_intersect_bbox_segment_aux - // inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf namespace CGAL { namespace internal { + template + inline + bool + bbox_line_do_intersect_aux(const FT& px, const FT& py, const FT& pz, + const FT& vx, const FT& vy, const FT& vz, + const FT& bxmin, const FT& bymin, const FT& bzmin, + const FT& bxmax, const FT& bymax, const FT& bzmax) + { + // ----------------------------------- + // treat x coord + // ----------------------------------- + FT dmin, tmin, tmax; + if ( vx >= 0 ) + { + tmin = bxmin - px; + tmax = bxmax - px; + dmin = vx; + } + else + { + tmin = px - bxmax; + tmax = px - bxmin; + dmin = -vx; + } + + //if px is not in the x-slab + if ( dmin == FT(0) && (tmin > FT(0) || tmax < FT(0)) ) return false; + + FT dmax = dmin; + + // ----------------------------------- + // treat y coord + // ----------------------------------- + FT d_, tmin_, tmax_; + if ( vy >= 0 ) + { + tmin_ = bymin - py; + tmax_ = bymax - py; + d_ = vy; + } + else + { + tmin_ = py - bymax; + tmax_ = py - bymin; + d_ = -vy; + } + + + + + if ( d_ == FT(0) ){ + //if py is not in the y-slab + if( (tmin_ > FT(0) || tmax_ < FT(0)) ) return false; + } + else + if ( (dmin*tmax_) < (d_*tmin) || (dmax*tmin_) > (d_*tmax) ) + return false; + + if( (dmin*tmin_) > (d_*tmin) ) + { + tmin = tmin_; + dmin = d_; + } + + if( (dmax*tmax_) < (d_*tmax) ) + { + tmax = tmax_; + dmax = d_; + } + + // ----------------------------------- + // treat z coord + // ----------------------------------- + if ( vz >= 0 ) + { + tmin_ = bzmin - pz; + tmax_ = bzmax - pz; + d_ = vz; + } + else + { + tmin_ = pz - bzmax; + tmax_ = pz - bzmin; + d_ = -vz; + } + + //if pz is not in the z-slab + //if ( d_ == FT(0) && (tmin_ > FT(0) || tmax_ < FT(0)) ) return false; + //The previous line is not needed as either dmin or d_ are not 0 + //(otherwise the direction of the line would be null). + // The following is equivalent to the in z-slab test if d_=0. + + return ( (dmin*tmax_) >= (d_*tmin) && (dmax*tmin_) <= (d_*tmax) ); + } + template bool do_intersect(const typename K::Line_3& line, const CGAL::Bbox_3& bbox, @@ -47,16 +139,11 @@ namespace internal { const Point_3& point = line.point(); const Vector_3& v = line.to_vector(); - return do_intersect_bbox_segment_aux - // do not use static filters - ( - point.x(), point.y(), point.z(), - v.x(), v.y(), v.z(), - bbox - ); + return bbox_line_do_intersect_aux( + point.x(), point.y(), point.z(), + v.x(), v.y(), v.z(), + FT(bbox.xmin()), FT(bbox.ymin()), FT(bbox.zmin()), + FT(bbox.xmax()), FT(bbox.ymax()), FT(bbox.zmax()) ); } } // namespace internal From 28c9507c17f3b09cabae905411ebc5ae79a87613 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 15 Jun 2012 15:38:45 +0000 Subject: [PATCH 30/38] Go back to the version of revision 69620. That was faster. --- .../Bbox_3_Segment_3_do_intersect.h | 166 ++++++++++++------ 1 file changed, 111 insertions(+), 55 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index d00814f912c..c84a850dbfb 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -139,31 +139,50 @@ namespace internal { // ----------------------------------- // treat x coord // ----------------------------------- - FT dmin, tmin, tmax, dmax, tmin_minus_dmin, tmax_minus_dmax; + FT dmin, tmin, tmax, dmax; if ( qx >= px ) { - tmin = bxmin - px; - tmax = bxmax - px; - dmax = dmin = qx - px; - tmin_minus_dmin = bxmin - qx; - tmax_minus_dmax = bxmax - qx; + if(bounded_0 && px > bxmax) return false; // segment on the right of bbox + if(bounded_1 && qx < bxmin) return false; // segment on the left of bbox + + if(bounded_1 && bxmax > qx) { + tmax = 1; + dmax = 1; + } else { + tmax = bxmax - px; + dmax = qx - px; + } + + if(bounded_0 && bxmin < px) // tmin < 0 means px is in the x-range of bbox + { + tmin = 0; + dmin = 1; + } else { + tmin = bxmin - px; + dmin = qx - px; + } } else { - tmin = px - bxmax; - tmax = px - bxmin; - dmax = dmin = px - qx; - tmin_minus_dmin = qx - bxmax; - tmax_minus_dmax = qx - bxmin; - } + if(bounded_1 && qx > bxmax) return false; // segment on the right of bbox + if(bounded_0 && px < bxmin) return false; // segment on the left of bbox - if( bounded_0 && tmax < 0. ) return false; // t2 < 0 - if( bounded_1 && tmin_minus_dmin > 0. ) return false; // t1 > 1 + if(bounded_1 && bxmin < qx) { + tmax = 1; + dmax = 1; + } else { + tmax = px - bxmin; + dmax = px - qx; + } - if(bounded_0) tmin = CGAL::max(tmin, FT(0)); // t1 = max(t1,0) - if(bounded_1 && tmax_minus_dmax > 0) { // t2 = min(t2,1); - tmax = FT(1); - dmax = FT(1); + if(bounded_0 && px < bxmax) // tmin < 0 means px is in the x-range of bbox + { + tmin = 0; + dmin = 1; + } else { + tmin = px - bxmax; + dmin = px - qx; + } } // If the query is vertical for x, then check its x-coordinate is in @@ -188,31 +207,50 @@ namespace internal { // ----------------------------------- // treat y coord // ----------------------------------- - FT dymin, tymin, tymax, dymax, tymin_minus_dmin, tymax_minus_dmax; + FT dymin, tymin, tymax, dymax; if ( qy >= py ) { - tymin = bymin - py; - tymax = bymax - py; - dymax = dymin = qy - py; - tymin_minus_dmin = bymin - qy; - tymax_minus_dmax = bymax - qy; + if(bounded_0 && py > bymax) return false; // segment on the right of bbox + if(bounded_1 && qy < bymin) return false; // segment on the left of bbox + + if(bounded_1 && bymax > qy) { + tymax = 1; + dymax = 1; + } else { + tymax = bymax - py; + dymax = qy - py; + } + + if(bounded_0 && bymin < py) // tmin < 0 means py is in the y-range of bbox + { + tymin = 0; + dymin = 1; + } else { + tymin = bymin - py; + dymin = qy - py; + } } else { - tymin = py - bymax; - tymax = py - bymin; - dymax = dymin = py - qy; - tymin_minus_dmin = qy - bymax; - tymax_minus_dmax = qy - bymin; - } + if(bounded_1 && qy > bymax) return false; // segment on the right of bbox + if(bounded_0 && py < bymin) return false; // segment on the left of bbox - if( bounded_0 && tymax < 0. ) return false; // t2 < 0 - if( bounded_1 && tymin_minus_dmin > 0. ) return false; // t1 > 1 + if(bounded_1 && bymin < qy) { + tymax = 1; + dymax = 1; + } else { + tymax = py - bymin; + dymax = py - qy; + } - if(bounded_0) tymin = CGAL::max(tymin, FT(0)); // t1 = max(t1,0) - if(bounded_1 && tymax_minus_dmax > 0) { // t2 = min(t2,1); - tymax = FT(1); - dymax = FT(1); + if(bounded_0 && py < bymax) // tmin < 0 means py is in the y-range of bbox + { + tymin = 0; + dymin = 1; + } else { + tymin = py - bymax; + dymin = py - qy; + } } // If the query is vertical for y, then check its y-coordinate is in @@ -236,32 +274,50 @@ namespace internal { // ----------------------------------- // treat z coord // ----------------------------------- - FT dzmin, tzmin, tzmax, dzmax, tzmax_minus_dmax, tzmin_minus_dmin; - + FT dzmin, tzmin, tzmax, dzmax; if ( qz >= pz ) { - tzmin = bzmin - pz; - tzmax = bzmax - pz; - dzmax = dzmin = qz - pz; - tzmin_minus_dmin = bzmin - qz; - tzmax_minus_dmax = bzmax - qz; + if(bounded_0 && pz > bzmax) return false; // segment on the right of bbox + if(bounded_1 && qz < bzmin) return false; // segment on the left of bbox + + if(bounded_1 && bzmax > qz) { + tzmax = 1; + dzmax = 1; + } else { + tzmax = bzmax - pz; + dzmax = qz - pz; + } + + if(bounded_0 && bzmin < pz) // tmin < 0 means pz is in the z-range of bbox + { + tzmin = 0; + dzmin = 1; + } else { + tzmin = bzmin - pz; + dzmin = qz - pz; + } } else { - tzmin = pz - bzmax; - tzmax = pz - bzmin; - dzmax = dzmin = pz - qz; - tzmin_minus_dmin = qz - bzmax; - tzmax_minus_dmax = qz - bzmin; - } + if(bounded_1 && qz > bzmax) return false; // segment on the right of bbox + if(bounded_0 && pz < bzmin) return false; // segment on the left of bbox - if( bounded_0 && tzmax < 0. ) return false; // t2 < 0 - if( bounded_1 && tzmin_minus_dmin > 0. ) return false; // t1 > 1 + if(bounded_1 && bzmin < qz) { + tzmax = 1; + dzmax = 1; + } else { + tzmax = pz - bzmin; + dzmax = pz - qz; + } - if(bounded_0) tzmin = CGAL::max(tzmin, FT(0)); // t1 = max(t1,0) - if(bounded_1 && tzmax_minus_dmax > 0) { // t2 = min(t2,1); - tzmax = FT(1); - dzmax = FT(1); + if(bounded_0 && pz < bzmax) // tmin < 0 means pz is in the z-range of bbox + { + tzmin = 0; + dzmin = 1; + } else { + tzmin = pz - bzmax; + dzmin = pz - qz; + } } // If the query is vertical for z, then check its z-coordinate is in From 6c240a4ae24d82c8106604c5407ddeeedd44cc96 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 15 Jun 2012 15:49:33 +0000 Subject: [PATCH 31/38] Cleanup: remove unused code. --- .../Bbox_3_Segment_3_do_intersect.h | 49 ------------------- 1 file changed, 49 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index c84a850dbfb..ee7a2205a50 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -116,7 +116,6 @@ namespace internal { bool bounded_0, bool bounded_1, bool use_static_filters> -// __attribute__ ((noinline)) inline typename Do_intersect_bbox_segment_aux_is_greater::result_type do_intersect_bbox_segment_aux( @@ -201,8 +200,6 @@ namespace internal { CGAL_assertion(dmin >= 0); CGAL_assertion(dmax >= 0); - // CGAL_assertion(!bounded_0 || ( (dmin == 0) == (px == qx && (px == bxmax || px == bxmin)) ) ); - // CGAL_assertion(!bounded_1 || ( (dmax == 0) == (px == qx && (px == bxmax || px == bxmin)) ) ); // ----------------------------------- // treat y coord @@ -267,8 +264,6 @@ namespace internal { CGAL_assertion(dymin >= 0); CGAL_assertion(dymax >= 0); - // CGAL_assertion(!bounded_0 || (dmin == 0) == (!bounded_0 && px == qx && py == qy)); - // CGAL_assertion(!bounded_1 || (dmax == 0) == (!bounded_1 && px == qx && py == qy)); // ----------------------------------- @@ -382,8 +377,6 @@ namespace internal { CGAL_assertion(dmin >= 0); CGAL_assertion(dmax >= 0); - // CGAL_assertion((dmin == 0) == (px == qx && py == qy)); - // CGAL_assertion((dmax == 0) == (px == qx && py == qy)); // If t1 > tzmax || tzmin > t2, return false. if( (px != qx || @@ -405,40 +398,6 @@ namespace internal { return true; } - template - inline - bool - do_intersect_bbox_segment_aux(const CGAL::cpp0x::array seg, - const CGAL::cpp0x::array box) - - { - const FT& px = seg[0]; - const FT& py = seg[1]; - const FT& pz = seg[2]; - const FT& qx = seg[3]; - const FT& qy = seg[4]; - const FT& qz = seg[5]; - const double& bxmin = box[0]; - const double& bymin = box[1]; - const double& bzmin = box[2]; - const double& bxmax = box[3]; - const double& bymax = box[4]; - const double& bzmax = box[5]; - // for(int i = 0; i < 3; ++i) { - // const int sign = seg[3+i] > seg[i]; // (qx > px)? - // if(bounded_0 && seg[3*(1-sign) + i] > box[3+i]) return false; // segment on the right of bbox - // if(bounded_1 && seg[3*sign + i] < box[i]) return false; // segment on the left of bbox - // } - - return do_intersect_bbox_segment_aux - (px, py, pz, - qy, qy, qz, - bxmin, bymin, bymax, - bxmax, bymax, bzmax); - } - template bool do_intersect(const typename K::Segment_3& segment, const CGAL::Bbox_3& bbox, @@ -454,14 +413,6 @@ namespace internal { source.x(), source.y(), source.z(), target.x(), target.y(), target.z(), bbox); - // bbox.xmin(), bbox.ymin(), bbox.zmin(), - // bbox.xmax(), bbox.ymax(), bbox.zmax() ); - - // const CGAL::cpp0x::array seg = {source.x(), source.y(), source.z(), - // target.x(), target.y(), target.z() }; - // return do_intersect_bbox_segment_aux - // ( seg, - // *reinterpret_cast*>(&*bbox.cartesian_begin()) ); } template From 090b93f2293b79ea2d8a687044253d5d48517f70 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Mon, 18 Jun 2012 09:41:39 +0000 Subject: [PATCH 32/38] Fix compilation for MSVC --- .../Intersections_3/Bbox_3_Segment_3_do_intersect.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index ee7a2205a50..44e3b9d352b 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -69,10 +69,6 @@ namespace internal { double tmax; double dmax; - static const double EPS = 8.8872057372592798e-16; - static const double OVERF = 1e153; - static const double UNDERF = 1e-146; - public: CGAL_static_assertion((boost::is_same::value)); @@ -86,14 +82,17 @@ namespace internal { } void compute_new_error_bound() { + const double EPS = 8.8872057372592798e-16; error = tmax * dmax * EPS; } bool bound_overflow() { + const double OVERF = 1e153; return dmax > OVERF && tmax > OVERF; } bool value_might_underflow() { + const double UNDERF = 1e-146; return dmax < UNDERF || dmax < UNDERF; } From ad2b26dbeddacca79cb0ea93a02a0b57b54590ee Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Mon, 18 Jun 2012 13:04:53 +0000 Subject: [PATCH 33/38] Fix a warning about unused variables And fix a comment (t=1 instead of t=0). --- .../CGAL/internal/Static_filters/Do_intersect_3.h | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index 1961ad0b438..667d65aff34 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -114,9 +114,6 @@ public: const Point_3& q = s.target(); double px, py, pz, qx, qy, qz; - double bxmin = b.xmin(), bymin = b.ymin(), bzmin = b.zmin(), - bxmax = b.xmax(), bymax = b.ymax(), bzmax = b.zmax(); - if (fit_in_double(get_approx(p).x(), px) && fit_in_double(get_approx(p).y(), py) && fit_in_double(get_approx(p).z(), pz) && fit_in_double(get_approx(q).x(), qx) && fit_in_double(get_approx(q).y(), qy) && @@ -128,7 +125,7 @@ public: do_intersect_bbox_segment_aux // do use static filters (px, py, pz, qx, qy, qz, @@ -160,9 +157,6 @@ public: const Point_3& q = r.second_point(); double px, py, pz, qx, qy, qz; - double bxmin = b.xmin(), bymin = b.ymin(), bzmin = b.zmin(), - bxmax = b.xmax(), bymax = b.ymax(), bzmax = b.zmax(); - if (fit_in_double(get_approx(p).x(), px) && fit_in_double(get_approx(p).y(), py) && fit_in_double(get_approx(p).z(), pz) && fit_in_double(get_approx(q).x(), qx) && fit_in_double(get_approx(q).y(), qy) && @@ -174,7 +168,7 @@ public: do_intersect_bbox_segment_aux // do use static filters (px, py, pz, qx, qy, qz, From 5142588eb3a4641e22d4320877192915d75ed57f Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Mon, 18 Jun 2012 16:05:04 +0000 Subject: [PATCH 34/38] Better performance (~10-15% better) By removing several tests (and use CGAL::max instead), the generated assembly is more efficient. --- .../Bbox_3_Segment_3_do_intersect.h | 66 +++++-------------- 1 file changed, 18 insertions(+), 48 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 44e3b9d352b..992e6b0c17c 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -151,14 +151,8 @@ namespace internal { dmax = qx - px; } - if(bounded_0 && bxmin < px) // tmin < 0 means px is in the x-range of bbox - { - tmin = 0; - dmin = 1; - } else { - tmin = bxmin - px; - dmin = qx - px; - } + tmin = bxmin - px; + dmin = qx - px; } else { @@ -173,16 +167,12 @@ namespace internal { dmax = px - qx; } - if(bounded_0 && px < bxmax) // tmin < 0 means px is in the x-range of bbox - { - tmin = 0; - dmin = 1; - } else { - tmin = px - bxmax; - dmin = px - qx; - } + tmin = px - bxmax; + dmin = px - qx; } + if(bounded_0) tmin = CGAL::max(FT(0), tmin); + // If the query is vertical for x, then check its x-coordinate is in // the x-slab. if( (px == qx) && // <=> (dmin == 0) @@ -217,14 +207,8 @@ namespace internal { dymax = qy - py; } - if(bounded_0 && bymin < py) // tmin < 0 means py is in the y-range of bbox - { - tymin = 0; - dymin = 1; - } else { - tymin = bymin - py; - dymin = qy - py; - } + tymin = bymin - py; + dymin = qy - py; } else { @@ -239,16 +223,12 @@ namespace internal { dymax = py - qy; } - if(bounded_0 && py < bymax) // tmin < 0 means py is in the y-range of bbox - { - tymin = 0; - dymin = 1; - } else { - tymin = py - bymax; - dymin = py - qy; - } + tymin = py - bymax; + dymin = py - qy; } + if(bounded_0) tymin = CGAL::max(FT(0), tymin); + // If the query is vertical for y, then check its y-coordinate is in // the y-slab. if( (py == qy) && // <=> (dmin == 0) @@ -282,14 +262,8 @@ namespace internal { dzmax = qz - pz; } - if(bounded_0 && bzmin < pz) // tmin < 0 means pz is in the z-range of bbox - { - tzmin = 0; - dzmin = 1; - } else { - tzmin = bzmin - pz; - dzmin = qz - pz; - } + tzmin = bzmin - pz; + dzmin = qz - pz; } else { @@ -304,16 +278,12 @@ namespace internal { dzmax = pz - qz; } - if(bounded_0 && pz < bzmax) // tmin < 0 means pz is in the z-range of bbox - { - tzmin = 0; - dzmin = 1; - } else { - tzmin = pz - bzmax; - dzmin = pz - qz; - } + tzmin = pz - bzmax; + dzmin = pz - qz; } + if(bounded_0) tzmin = CGAL::max(FT(0), tzmin); + // If the query is vertical for z, then check its z-coordinate is in // the z-slab. if( (pz == qz) && // <=> (dmin == 0) From ab28ccf7402c7d35196b801be17e4f6282948728 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Mon, 18 Jun 2012 16:05:29 +0000 Subject: [PATCH 35/38] Fix bugs in overflow/underflow handling --- .../internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 992e6b0c17c..25961025f06 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -88,12 +88,12 @@ namespace internal { bool bound_overflow() { const double OVERF = 1e153; - return dmax > OVERF && tmax > OVERF; + return dmax > OVERF || tmax > OVERF; } bool value_might_underflow() { const double UNDERF = 1e-146; - return dmax < UNDERF || dmax < UNDERF; + return dmax < UNDERF || tmax < UNDERF; } typedef Uncertain result_type; From 23ab39f2dff910731065ea37cbfa8d185dcf7ecf Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 19 Jun 2012 10:08:04 +0000 Subject: [PATCH 36/38] Remove calls to CGAL::abs for Ray_3 and Segment_3 As t1 and t2 are, by construction, greater than 0, for rays and segments, no need to call CGAL::abs on t[xyz]min and t[xyz]max. --- .../Bbox_3_Segment_3_do_intersect.h | 47 +++++++++++++++---- 1 file changed, 37 insertions(+), 10 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 25961025f06..60f4794c0b8 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -43,7 +43,7 @@ namespace CGAL { namespace internal { - template + template struct Do_intersect_bbox_segment_aux_is_greater { typedef typename Same_uncertainty::type result_type; @@ -62,8 +62,8 @@ namespace internal { } }; // end struct template Do_intersect_bbox_segment_aux_is_greater - template - class Do_intersect_bbox_segment_aux_is_greater + template + class Do_intersect_bbox_segment_aux_is_greater { double error; double tmax; @@ -75,10 +75,15 @@ namespace internal { Do_intersect_bbox_segment_aux_is_greater() : error(0.), tmax(0.), dmax(0.) {} void register_new_input_values(const double& t, const double& d) { - const double at = CGAL::abs(t); - const double ad = CGAL::abs(d); - if(at > tmax) tmax = at; - if(ad > dmax) dmax = ad; + if(bounded_0) { + if(t > tmax) tmax = t; + if(d > dmax) dmax = d; + } else { + const double at = CGAL::abs(t); + const double ad = CGAL::abs(d); + if(at > tmax) tmax = at; + if(ad > dmax) dmax = ad; + } } void compute_new_error_bound() { @@ -116,7 +121,12 @@ namespace internal { bool bounded_1, bool use_static_filters> inline - typename Do_intersect_bbox_segment_aux_is_greater::result_type + typename Do_intersect_bbox_segment_aux_is_greater + < + FT, + bounded_0, + use_static_filters + >::result_type do_intersect_bbox_segment_aux( const FT& px, const FT& py, const FT& pz, const FT& qx, const FT& qy, const FT& qz, @@ -189,6 +199,11 @@ namespace internal { CGAL_assertion(dmin >= 0); CGAL_assertion(dmax >= 0); + if(bounded_0) { + CGAL_assertion(tmin >= 0); + CGAL_assertion(tmax >= 0); + } + // ----------------------------------- // treat y coord @@ -243,6 +258,10 @@ namespace internal { CGAL_assertion(dymin >= 0); CGAL_assertion(dymax >= 0); + if(bounded_0) { + CGAL_assertion(tymin >= 0); + CGAL_assertion(tymax >= 0); + } // ----------------------------------- @@ -298,9 +317,17 @@ namespace internal { CGAL_assertion(dzmin >= 0); CGAL_assertion(dzmax >= 0); + if(bounded_0) { + CGAL_assertion(tzmin >= 0); + CGAL_assertion(tzmax >= 0); + } - typedef Do_intersect_bbox_segment_aux_is_greater Is_greater; + + typedef Do_intersect_bbox_segment_aux_is_greater + Is_greater; typedef typename Is_greater::result_type Is_greater_value; Is_greater is_greater; From eecf1e94be45869135ad505e2187879d5787b299 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 20 Jun 2012 09:53:56 +0000 Subject: [PATCH 37/38] Fix for MSVC: protect the CGAL::max calls to prevent macro substitution --- .../Intersections_3/Bbox_3_Segment_3_do_intersect.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 60f4794c0b8..03a8fadb79f 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -181,7 +181,7 @@ namespace internal { dmin = px - qx; } - if(bounded_0) tmin = CGAL::max(FT(0), tmin); + if(bounded_0) tmin = (CGAL::max)(FT(0), tmin); // If the query is vertical for x, then check its x-coordinate is in // the x-slab. @@ -242,7 +242,7 @@ namespace internal { dymin = py - qy; } - if(bounded_0) tymin = CGAL::max(FT(0), tymin); + if(bounded_0) tymin = (CGAL::max)(FT(0), tymin); // If the query is vertical for y, then check its y-coordinate is in // the y-slab. @@ -301,7 +301,7 @@ namespace internal { dzmin = pz - qz; } - if(bounded_0) tzmin = CGAL::max(FT(0), tzmin); + if(bounded_0) tzmin = (CGAL::max)(FT(0), tzmin); // If the query is vertical for z, then check its z-coordinate is in // the z-slab. From 19ab08d6e95bdda237c6cb99de3e50cda792babe Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 20 Jun 2012 10:32:46 +0000 Subject: [PATCH 38/38] Fix a warning from clang The warning was: [...]/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h:353:9: warning: variable 'b' is used uninitialized whenever '||' condition is true [-Wsometimes-uninitialized] if( (px == qx) || // <=> (dmin == 0) ^~~~~~~~~~ [...]/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h:408:12: note: in instantiation of function template specialization 'CGAL::internal::do_intersect_bbox_segment_aux' requested here return do_intersect_bbox_segment_aux( ^ [...]/include/CGAL/Kernel/function_objects.h:2052:14: note: in instantiation of function template specialization 'CGAL::internal::do_intersect >' requested here { return internal::do_intersect(t1, t2, K()); } ^ [...]/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h:428:10: note: in instantiation of function template specialization 'CGAL::CommonKernelFunctors::Do_intersect_3 >::operator()>, CGAL::Bbox_3>' requested here return typename K::Do_intersect_3()(segment, bbox); ^ [...]/cmake/platforms/x86-64_Linux-2.6_llvm-clang-with-g++-4.6.2_F16/test/Intersections_3/bbox_other_do_intersect_test.cpp:67:12: note: in instantiation of function template specialization 'CGAL::do_intersect >' requested here bool b = CGAL::do_intersect(t,bbox); ^ [...]/cmake/platforms/x86-64_Linux-2.6_llvm-clang-with-g++-4.6.2_F16/test/Intersections_3/bbox_other_do_intersect_test.cpp:416:12: note: in instantiation of function template specialization 'test_aux> >' requested here bool b = test_aux(s12,"s12",bbox,false); ^ [...]/cmake/platforms/x86-64_Linux-2.6_llvm-clang-with-g++-4.6.2_F16/test/Intersections_3/bbox_other_do_intersect_test.cpp:665:12: note: in instantiation of function template specialization 'test >' requested here bool b = test(exact_predicates) && ^ [...]/cmake/platforms/x86-64_Linux-2.6_llvm-clang-with-g++-4.6.2_F16/test/Intersections_3/bbox_other_do_intersect_test.cpp:678:7: note: in instantiation of function template specialization 'test_kernel >' requested here b = test_kernel >(false); ^ [...]/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h:360:25: note: uninitialized use occurs here if(is_indeterminate(b)) return b; // Note that the default-constructed ^ [...]/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h:353:9: note: remove the '||' if its condition is always false if( (px == qx) || // <=> (dmin == 0) ^~~~~~~~~~~~~ [...]/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h:351:23: note: initialize the variable 'b' to silence this warning Is_greater_value b; ^ = false --- .../internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h index 03a8fadb79f..41f958ddc5b 100644 --- a/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h +++ b/Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h @@ -348,7 +348,7 @@ namespace internal { if(possibly(b2)) return !b2; } - Is_greater_value b; + Is_greater_value b = Is_greater_value(); // If tymin/dymin > t1, set t1 = tymin/dymin. if( (px == qx) || // <=> (dmin == 0) ( (py != qy) && // <=> (dymin > 0)