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 55abdaa8d0b..fd14bab2aa1 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 @@ -289,6 +289,232 @@ public: return Base::operator()(s,b); } + + + result_type + operator()(const Ray_3 &r, const Bbox_3& b) const + { + CGAL_BRANCH_PROFILER_3("semi-static failures/attempts/calls to : Do_intersect_3", tmp); + + Get_approx get_approx; // Identity functor for all points + // but lazy points. + const Point_3& p = r.source(); + const Point_3& q = r.point(1); + + 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) && + fit_in_double(get_approx(q).z(), qz) ) + { + // std::cerr << "\n" + // << px << " " << py << " " << pz << "\n" + // << qx << " " << qy << " " << qz << "\n" + // << bxmin << " " << bymin << " " << bzmin << "\n" + // << bxmax << " " << bymax << " " << bzmax << "\n"; + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + + + // AF: I copy pasted the code to call the max + + // ----------------------------------- + // treat x coord + // ----------------------------------- + double dmin, dmax, tmin, tmax; + if ( qx >= px ) // this is input and needs no epsilon + { + // different from segment: if(px > bxmax) return false; // segment on the right of bbox + // if(qx < bxmin) return false; // segment on the left of bbox + if (bxmax < px) return false; // different from segment: added + 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 + { + // different from segment: if(qx > bxmax) return false; // segment on the right of bbox + // if(px < bxmin) return false; // segment on the left of bbox + if(px < bxmin) return false; // different from segment: added + 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; + } + } + + // std::cerr << "t1 "; + + 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()(r,b); + } + + const double EPS_1 = 3.55618e-15; + + double error = EPS_1 * m; + + switch(sign_with_error( tmax - dmax, error)) { + case POSITIVE: + // std::cerr << "t2 "; + tmax = 1; + dmax = 1; + break; + case NEGATIVE: + break; + default: + // ambiguity of the comparison tmax > dmin + // let's call the exact predicate + // std::cerr << "\ntest2 NEED EXACT\n"; + return Base::operator()(r,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()(r,b); + } + + error = EPS_1 * m * m; + + // std::cerr << dmin << " " << tmax_ << " " << d_ << " " + // << tmin << " " << dmax << " " << tmin_ << std::endl; + + 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 segment is outside the box, on one + // side or the other. + if(sign1 == ZERO || sign2 == ZERO) { + // std::cerr << "\ntest3 NEED EXACT\n"; + return Base::operator()(r,b); // We are *unsure*: one *may be* + // positive. + } + + // std::cerr << "t3 "; + + // 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_; + // std::cerr << "t4 "; + break; + case NEGATIVE: + break; + default: // uncertainty + // std::cerr << "\ntest4 NEED EXACT\n"; + return Base::operator()(r,b); + } + + // epsilons needed + switch(sign_with_error((d_*tmax) - (dmax*tmax_) , error)) { + case POSITIVE: + tmax = tmax_; + dmax = d_; + // std::cerr << "t5 ";break; + case NEGATIVE: + break; + default: // uncertainty + // std::cerr << "\ntest5 NEED EXACT\n"; + 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) { + // std::cerr << "f6"; + return false; // We are *sure* the segment is outside the box, on one + // side or the other. + } + if(sign1 == ZERO || sign2 == ZERO) { + // std::cerr << "\test6 NEED EXACT\n"; + return Base::operator()(r,b); // We are *unsure*: one *may be* + // negative. + } + // std::cerr << "t6"; + return true; // We are *sure* the two signs are positive. + } + return Base::operator()(r,b); + } + + + // Computes the epsilon for Bbox_3_Segment_3_do_intersect. static double compute_epsilon_bbox_segment_3() {