diff --git a/.gitignore b/.gitignore index 1a5f61c990b..34ebe46c01c 100644 --- a/.gitignore +++ b/.gitignore @@ -168,6 +168,13 @@ Interpolation/demo/Interpolation/cgal_test_with_cmake 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/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..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 @@ -28,7 +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 @@ -111,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) && @@ -121,178 +121,18 @@ 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; // segment on the right of bbox - if(qx < bxmin) return false; // segment on the left of bbox + const Uncertain ub = + do_intersect_bbox_segment_aux + // do use static filters + (px, py, pz, + qx, qy, qz, + b); - 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. + if(is_indeterminate(ub)) return Base::operator()(s,b); + else return ub.sup(); } return Base::operator()(s,b); } @@ -317,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) && @@ -327,155 +164,18 @@ 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 + // do use static filters + (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); } 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..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 @@ -25,105 +25,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_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, @@ -135,11 +45,16 @@ namespace internal { const Point_3& source = ray.source(); const Point_3& point_on_ray = ray.second_point(); - return bbox_ray_do_intersect_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()) ); + 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 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..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 @@ -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 @@ -25,104 +24,374 @@ #include #include +#include +#include +#include // 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 bounding box. + +// For a segment, the intersection is non-empty iff +// [t1, t2] intersects [0, 1]. + namespace CGAL { namespace internal { - template + template + struct Do_intersect_bbox_segment_aux_is_greater + { + 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; + + 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) { + 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() { + 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 || tmax < 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 inline - bool + 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, - const FT& bxmin, const FT& bymin, const FT& bzmin, - const FT& bxmax, const FT& bymax, const FT& 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 + // For the first lines, dmax==dmin and is not explicitly defined. + // ----------------------------------- // treat x coord // ----------------------------------- - FT dmin, tmin, tmax; + FT dmin, tmin, tmax, 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; + } + tmin = bxmin - px; - tmax = bxmax - px; dmin = qx - px; } 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 + + if(bounded_1 && bxmin < qx) { + tmax = 1; + dmax = 1; + } else { + tmax = px - bxmin; + dmax = px - qx; + } + tmin = px - bxmax; - tmax = px - bxmin; dmin = px - qx; } - if ( tmax < FT(0) || tmin > dmin ) - return false; + if(bounded_0) tmin = (CGAL::max)(FT(0), tmin); - FT dmax = dmin; - if ( tmin < FT(0) ) + // 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) ) ) // do not check for a segment { - tmin = FT(0); - dmin = FT(1); + 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). } - if ( tmax > dmax ) - { - tmax = FT(1); - dmax = FT(1); + // 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(dmin >= 0); + CGAL_assertion(dmax >= 0); + if(bounded_0) { + CGAL_assertion(tmin >= 0); + CGAL_assertion(tmax >= 0); } + // ----------------------------------- // treat y coord // ----------------------------------- - FT d_, tmin_, tmax_; + 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; + } + + 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; + } + + tymin = py - bymax; + dymin = py - qy; } - if ( (dmin*tmax_) < (d_*tmin) || (dmax*tmin_) > (d_*tmax) ) - return false; + if(bounded_0) tymin = (CGAL::max)(FT(0), tymin); - if( (dmin*tmin_) > (d_*tmin) ) + // If the query is vertical for y, then check its y-coordinate is in + // the y-slab. + if( (py == qy) && // <=> (dmin == 0) + (! (bounded_0 && bounded_1) ) ) // do not check for a segment { - tmin = tmin_; - dmin = d_; + if(py > bymax || py < bymin) return false; } - if( (dmax*tmax_) < (d_*tmax) ) - { - tmax = tmax_; - dmax = d_; + // 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); + if(bounded_0) { + CGAL_assertion(tymin >= 0); + CGAL_assertion(tymax >= 0); } + // ----------------------------------- // treat z coord // ----------------------------------- + FT dzmin, tzmin, tzmax, dzmax; if ( qz >= pz ) { - tmin_ = bzmin - pz; - tmax_ = bzmax - pz; - d_ = 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; + } + + tzmin = bzmin - pz; + dzmin = qz - pz; } else { - tmin_ = pz - bzmax; - tmax_ = pz - bzmin; - d_ = pz - qz; + 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; + } + + tzmin = pz - bzmax; + dzmin = pz - qz; } - return ( (dmin*tmax_) >= (d_*tmin) && (dmax*tmin_) <= (d_*tmax) ); + 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) + (! (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(bounded_0) { + CGAL_assertion(tzmin >= 0); + CGAL_assertion(tzmax >= 0); + } + + + 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 + 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); + if(possibly(b2)) return !b2; + } + + Is_greater_value b = Is_greater_value(); + // If tymin/dymin > t1, set t1 = tymin/dymin. + if( (px == qx) || // <=> (dmin == 0) + ( (py != qy) && // <=> (dymin > 0) + certainly(b = is_greater( dmin*tymin, dymin* tmin)) ) ) + { + tmin = tymin; + dmin = dymin; + } + 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) + ( (py != qy) && // <=> dymax > 0 + 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); + + // If t1 > tzmax || tzmin > t2, return false. + if( (px != qx || + 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); + if(possibly(b2)) return !b2; // if(is_greater) return false; // or uncertain + } + return true; } template @@ -136,11 +405,10 @@ 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()), - FT(bbox.xmax()), FT(bbox.ymax()), FT(bbox.zmax()) ); + bbox); } 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 936afd1b6ed..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 @@ -37,6 +37,7 @@ #include +#include // for nextafter double random_in(const double a, @@ -61,17 +62,211 @@ 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); 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); } +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 exact_k = false, + const bool exactness_issue = false, + 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, exact_k, exactness_issue, false); + b &= test_case(-px, py, pz, + -qx, qy, qz, + -bxmin, bymin, bzmin, + -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, exact_k, exactness_issue, false); + b &= test_case( px, py, -pz, + qx, qy, -qz, + bxmin, bymin, -bzmin, + 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, exact_k, exactness_issue, false); + b &= test_case( px, -py, -pz, + qx, -qy, -qz, + bxmin, -bymin, -bzmin, + 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, exact_k, exactness_issue, false); + b &= test_case(-px, -py, -pz, + -qx, -qy, -qz, + -bxmin, -bymin, -bzmin, + -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, exact_k, exactness_issue, + false, false); + // xzy + b &= test_case( px, pz, py, + qx, qz, qy, + bxmin, bzmin, bymin, + 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, exact_k, exactness_issue, + false, false); + // zxy + b &= test_case( pz, px, py, + qz, qx, qy, + bzmin, bxmin, bymin, + 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, exact_k, exactness_issue, + false, false); + // zyx + b &= test_case( pz, py, px, + qz, qy, qx, + bzmin, bymin, bxmin, + 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, exact_k, exactness_issue, + false, false, false); + b &= test_case(qx, qy, qz, + px, py, pz, + bxmin, bymin, bzmin, + 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, 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, 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, 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, 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, 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, exact_k, exactness_issue, + false, false, false, false, false); + 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, exact_k, exactness_issue, + false, false, false, false, false); + 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, exact_k, exactness_issue, + false, false, false, false, false); + delta = 7; + 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, exact_k, exactness_issue, + false, false, false, false, false); + delta = 1; + 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, exact_k, exactness_issue, + false, false, 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) + { + 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; +} + template void speed(const std::string& name) { @@ -100,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 ) @@ -132,7 +327,50 @@ void test_speed() } template -bool test() +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, exact_predicates); + } + // 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, exact_predicates); + } + 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, exact_predicates); + } + return b; +} + +template +bool test(bool exact_kernel = false) { // types typedef typename K::FT FT; @@ -188,8 +426,48 @@ 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); + 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); + 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); + 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); @@ -309,22 +587,85 @@ 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(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, 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, exact_kernel, 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; +} + +template +bool test_kernel(bool exact_predicates = true, K k = K()) +{ + bool b = test(exact_predicates) && + intensive_test(exact_predicates); + test_speed(); + return b; } int main() @@ -332,38 +673,31 @@ int main() srand(0); std::cout << std::setprecision(5); + bool b; std::cout << "Testing with Simple_cartesian..." << std::endl ; - bool 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 >(true); std::cout << std::endl << "Testing with Simple_cartesian..." << std::endl ; - b &= test >(); - test_speed >(); + b &= test_kernel >(true); 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 >(true); 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;