From 0358937f0191b4ddaeeffe71a7a84176c2648b9e Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 20 Mar 2012 16:56:13 +0000 Subject: [PATCH] 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;