diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Do_intersect_3.h index 31f71e86b58b..44a3ae056690 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Do_intersect_3.h @@ -43,6 +43,7 @@ class Do_intersect_3 typedef typename K_base::Ray_3 Ray_3; typedef typename K_base::Segment_3 Segment_3; typedef typename K_base::Triangle_3 Triangle_3; + typedef typename K_base::Tetrahedron_3 Tetrahedron_3; typedef typename K_base::Sphere_3 Sphere_3; typedef typename K_base::Do_intersect_3 Base; @@ -126,7 +127,46 @@ class Do_intersect_3 return Base::operator()(s,b); } + result_type + operator()(const Bbox_3& b, const Tetrahedron_3 &t) const + { + return this->operator()(t, b); + } + + result_type + operator()(const Tetrahedron_3 &t, const Bbox_3& b) const + { + CGAL_BRANCH_PROFILER_3(std::string("semi-static failures/attempts/calls to : ") + + std::string(CGAL_PRETTY_FUNCTION), tmp); + Get_approx get_approx; + double px, py, pz; + + for(int i = 0; i < 4; ++i) + { + const Point_3& p = t[i]; + 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) ) + { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + + if( (px >= b.xmin() && px <= b.xmax()) && + (py >= b.ymin() && py <= b.ymax()) && + (pz >= b.zmin() && pz <= b.zmax()) ) + { + return true; + } + + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + } + else + { + return Base::operator()(t,b); + } + } + + return Base::operator()(t,b); + } result_type operator()(const Bbox_3& b, const Ray_3 &r) const @@ -169,6 +209,396 @@ class Do_intersect_3 return Base::operator()(r,b); } + + result_type + operator()(const Bbox_3& b, const Triangle_3 &t) const + { + return this->operator()(t, b); + } + + Uncertain sign_of_minor(double px, double py, double qx, double qy, double rx, double ry) const + { + CGAL_BRANCH_PROFILER_3("certain / uncertain / calls to : sign_of_minor", tmp); + + double qx_px = (qx - px); + double ry_py = (ry - py); + double rx_px = (rx - px); + double qy_py = (qy - py); + Sign int_tmp_result; + double double_tmp_result; + double eps; + double_tmp_result = ((qx_px * ry_py) - (rx_px * qy_py)); + double max1 = CGAL::abs(qx_px); + if( (max1 < CGAL::abs(rx_px)) ) + { + max1 = CGAL::abs(rx_px); + } + double max2 = CGAL::abs(ry_py); + if( (max2 < CGAL::abs(qy_py)) ) + { + max2 = CGAL::abs(qy_py); + } + double lower_bound_1; + double upper_bound_1; + lower_bound_1 = max1; + upper_bound_1 = max1; + if( (max2 < lower_bound_1) ) + { + lower_bound_1 = max2; + } + else + { + if( (max2 > upper_bound_1) ) + { + upper_bound_1 = max2; + } + } + if( (lower_bound_1 < 5.00368081960964746551e-147) ) + { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + return Uncertain::indeterminate(); + } + else + { + if( (upper_bound_1 > 1.67597599124282389316e+153) ) + { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + return Uncertain::indeterminate(); + } + eps = (8.88720573725927976811e-16 * (max1 * max2)); + if( (double_tmp_result > eps) ) + { + int_tmp_result = POSITIVE; + } + else + { + if( (double_tmp_result < -eps) ) + { + int_tmp_result = NEGATIVE; + } + else + { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + return Uncertain::indeterminate(); + } + } + } + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return int_tmp_result; + }; + + bool get_cross_product_sign(const std::array< std::array, 3 >& pts, std::array& signs) const + { + CGAL_BRANCH_PROFILER_3("determinate / indeterminate / calls to : get_cross_product_sign", tmp); + Uncertain s = sign_of_minor(pts[0][1], pts[0][2], pts[1][1], pts[1][2], pts[2][1], pts[2][2]); + if (is_indeterminate(s)){ + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + return false; + } + signs[0]=make_certain(s); + s = sign_of_minor(pts[0][2], pts[0][0], pts[1][2], pts[1][0], pts[2][2], pts[2][0]); + if (is_indeterminate(s)){ + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + return false; + } + signs[1]=make_certain(s); + s = sign_of_minor(pts[0][0], pts[0][1], pts[1][0], pts[1][1], pts[2][0], pts[2][1]); + if (is_indeterminate(s)){ + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + return false; + } + signs[2]=make_certain(s); + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return true; + } + + // adaptation of the do-intersect code of Plane_3/Bbox_3 when we don't know the plane equation exactly + bool do_intersect_supporting_plane_bbox(const Triangle_3& t, const std::array< std::array, 3>& pts, const Bbox_3& bbox) const + { + // copy of the static filter of Orientation_3 skipping the fit_in_double + using Uncertain as result type + auto statically_filtered_orientation_3 = + [](const std::array< std::array, 3>& pts, double x, double y, double z) + -> Uncertain + { + CGAL_BRANCH_PROFILER(std::string("Plane-BBox_3 semi-static Orientation_3 calls to/failures : ") + + std::string(CGAL_PRETTY_FUNCTION), tmp); + + double pqx = pts[1][0] - pts[0][0]; + double pqy = pts[1][1] - pts[0][1]; + double pqz = pts[1][2] - pts[0][2]; + double prx = pts[2][0] - pts[0][0]; + double pry = pts[2][1] - pts[0][1]; + double prz = pts[2][2] - pts[0][2]; + double psx = x - pts[0][0]; + double psy = y - pts[0][1]; + double psz = z - pts[0][2]; + + // CGAL::abs uses fabs on platforms where it is faster than (a<0)?-a:a + // Then semi-static filter. + + double maxx = CGAL::abs(pqx); + double maxy = CGAL::abs(pqy); + double maxz = CGAL::abs(pqz); + + double aprx = CGAL::abs(prx); + double apsx = CGAL::abs(psx); + + double apry = CGAL::abs(pry); + double apsy = CGAL::abs(psy); + + double aprz = CGAL::abs(prz); + double apsz = CGAL::abs(psz); +#ifdef CGAL_USE_SSE2_MAX + CGAL::Max mmax; + + maxx = mmax(maxx, aprx, apsx); + maxy = mmax(maxy, apry, apsy); + maxz = mmax(maxz, aprz, apsz); +#else + if (maxx < aprx) maxx = aprx; + if (maxx < apsx) maxx = apsx; + if (maxy < apry) maxy = apry; + if (maxy < apsy) maxy = apsy; + if (maxz < aprz) maxz = aprz; + if (maxz < apsz) maxz = apsz; +#endif + double det = CGAL::determinant(pqx, pqy, pqz, + prx, pry, prz, + psx, psy, psz); + + double eps = 5.1107127829973299e-15 * maxx * maxy * maxz; + +#ifdef CGAL_USE_SSE2_MAX +#if 0 + CGAL::Min mmin; + double tmp = mmin(maxx, maxy, maxz); + maxz = mmax(maxx, maxy, maxz); + maxx = tmp; +#else + sse2minmax(maxx,maxy,maxz); + // maxy can contain ANY element +#endif +#else + // Sort maxx < maxy < maxz. + if (maxx > maxz) + std::swap(maxx, maxz); + if (maxy > maxz) + std::swap(maxy, maxz); + else if (maxy < maxx) + std::swap(maxx, maxy); +#endif + // Protect against underflow in the computation of eps. + if (maxx < 1e-97) /* cbrt(min_double/eps) */ { + if (maxx == 0) + return ZERO; + } + // Protect against overflow in the computation of det. + else if (maxz < 1e102) /* cbrt(max_double [hadamard]/4) */ { + + if (det > eps) return POSITIVE; + if (det < -eps) return NEGATIVE; + } + + CGAL_BRANCH_PROFILER_BRANCH(tmp); + return Uncertain::indeterminate(); + }; + + auto orientation = + [statically_filtered_orientation_3] + (const Triangle_3& t, const std::array< std::array, 3>& pts, double x, double y, double z) + -> Orientation + { + Uncertain res = statically_filtered_orientation_3(pts, x, y, z); + if (!is_indeterminate(res)) return make_certain(res); + typename K_base::Orientation_3 orient = K_base().orientation_3_object(); // skip the static filter and directly call the base + return orient(t[0], t[1], t[2], Point_3(x,y,z)); + }; + + std::array signs; + bool OK = get_cross_product_sign(pts, signs); + + if (OK) + { + // extract extreme directions (copy of the code of get_min_max from Bbox_3_Plane_3_do_intersect.h) + std::array p_min, p_max; + + if(signs[0] == POSITIVE) { + if(signs[1] == POSITIVE) { + if(signs[2] == POSITIVE) { + p_min = {bbox.xmin(), bbox.ymin(),bbox.zmin()}; + p_max = {bbox.xmax(), bbox.ymax(),bbox.zmax()}; + } else { + p_min = {bbox.xmin(), bbox.ymin(),bbox.zmax()}; + p_max = {bbox.xmax(), bbox.ymax(),bbox.zmin()}; + } + } else { + if(signs[2]==POSITIVE) { + p_min = {bbox.xmin(), bbox.ymax(),bbox.zmin()}; + p_max = {bbox.xmax(), bbox.ymin(),bbox.zmax()}; + } else { + p_min = {bbox.xmin(), bbox.ymax(),bbox.zmax()}; + p_max = {bbox.xmax(), bbox.ymin(),bbox.zmin()}; + } + } + } + else{ + if(signs[1]==POSITIVE) { + if(signs[2]==POSITIVE) { + p_min = {bbox.xmax(), bbox.ymin(),bbox.zmin()}; + p_max = {bbox.xmin(), bbox.ymax(),bbox.zmax()}; + } + else{ + p_min = {bbox.xmax(), bbox.ymin(),bbox.zmax()}; + p_max = {bbox.xmin(), bbox.ymax(),bbox.zmin()}; + } + } + else { + if(signs[2]==POSITIVE) { + p_min = {bbox.xmax(), bbox.ymax(),bbox.zmin()}; + p_max = {bbox.xmin(), bbox.ymin(),bbox.zmax()}; + } + else { + p_min = {bbox.xmax(), bbox.ymax(),bbox.zmax()}; + p_max = {bbox.xmin(), bbox.ymin(),bbox.zmin()}; + } + } + } + + return ! (orientation(t, pts, p_max[0], p_max[1], p_max[2]) == ON_NEGATIVE_SIDE || + orientation(t, pts, p_min[0], p_min[1], p_min[2]) == ON_POSITIVE_SIDE); + } + + CGAL::Orientation side = orientation(t, pts, bbox.xmin(), bbox.ymin(),bbox.zmin()); + if(side == COPLANAR) return true; + CGAL::Oriented_side s = orientation(t, pts, bbox.xmax(), bbox.ymax(),bbox.zmax()); + if(s != side) return true; + s = orientation(t, pts, bbox.xmin(), bbox.ymin(),bbox.zmax()); + if(s != side) return true; + s = orientation(t, pts, bbox.xmax(), bbox.ymax(),bbox.zmin()); + if(s != side) return true; + s = orientation(t, pts, bbox.xmin(), bbox.ymax(),bbox.zmin()); + if(s != side) return true; + s = orientation(t, pts, bbox.xmax(), bbox.ymin(),bbox.zmax()); + if(s != side) return true; + s = orientation(t, pts, bbox.xmin(), bbox.ymax(),bbox.zmax()); + if(s != side) return true; + s = orientation(t, pts, bbox.xmax(), bbox.ymin(),bbox.zmin()); + if(s != side) return true; + return false; + }; + + result_type + operator()(const Triangle_3 &t, const Bbox_3& b) const + { + CGAL_BRANCH_PROFILER_3(std::string("semi-static failures/attempts/calls to : ") + + std::string(CGAL_PRETTY_FUNCTION), tmp); + + // check if at least one triangle point is inside the bbox + Get_approx get_approx; + std::array< std::array, 3> pts; + + for (int i=0; i<3; ++i) + { + const Point_3& p = t[i]; + if (fit_in_double(get_approx(p).x(), pts[i][0]) && fit_in_double(get_approx(p).y(), pts[i][1]) && + fit_in_double(get_approx(p).z(), pts[i][2]) ) + { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + + if( (pts[i][0] >= b.xmin() && pts[i][0] <= b.xmax()) && + (pts[i][1] >= b.ymin() && pts[i][1] <= b.ymax()) && + (pts[i][2] >= b.zmin() && pts[i][2] <= b.zmax()) ) + { + return true; + } + + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + } + else + { + return Base::operator()(t,b); + } + } + + // copy of the regular code with do_axis_intersect_aux_impl statically filtered + auto do_axis_intersect_aux_impl = [](double alpha, double beta, double c_alpha, double c_beta) -> Uncertain + { + CGAL_BRANCH_PROFILER_3("certain / uncertain / calls to : axis_inter", tmp2); + + Sign int_tmp_result; + double double_tmp_result; + double eps; + double_tmp_result = ((-c_alpha * alpha) + (c_beta * beta)); + double max1 = CGAL::abs(c_alpha); + if( (max1 < CGAL::abs(c_beta)) ) + { + max1 = CGAL::abs(c_beta); + } + double max2 = CGAL::abs(alpha); + if( (max2 < CGAL::abs(beta)) ) + { + max2 = CGAL::abs(beta); + } + double lower_bound_1; + double upper_bound_1; + lower_bound_1 = max1; + upper_bound_1 = max1; + if( (max2 < lower_bound_1) ) + { + lower_bound_1 = max2; + } + else + { + if( (max2 > upper_bound_1) ) + { + upper_bound_1 = max2; + } + } + if( (lower_bound_1 < 5.00368081960964746551e-147) ) + { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp2); + return Uncertain::indeterminate(); + } + else + { + if( (upper_bound_1 > 1.67597599124282389316e+153) ) + { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp2); + return Uncertain::indeterminate(); + } + eps = (8.88720573725927976811e-16 * (max1 * max2)); + if( (double_tmp_result > eps) ) + { + int_tmp_result = POSITIVE; + } + else + { + if( (double_tmp_result < -eps) ) + { + int_tmp_result = NEGATIVE; + } + else + { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp2); + return Uncertain::indeterminate(); + } + } + } + + CGAL_BRANCH_PROFILER_BRANCH_2(tmp2); + return int_tmp_result; + }; + + if ( !do_intersect_supporting_plane_bbox(t, pts, b) ) + return false; + + Uncertain res = Intersections::internal::do_intersect_bbox_or_iso_cuboid_impl(pts, b, do_axis_intersect_aux_impl); + if ( !is_indeterminate(res) ) + return make_certain(res); + + return Base::operator()(t,b); + } + result_type operator()(const Bbox_3& b, const Sphere_3 &s) const { diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Triangle_3_do_intersect.h b/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Triangle_3_do_intersect.h index ea709e2c9c30..8eca9e6a5172 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Triangle_3_do_intersect.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Triangle_3_do_intersect.h @@ -16,12 +16,15 @@ // Fast Triangle-Cuboid intersection test, following Tomas Akenine-Moeller description. // The code looks slightly different from his code because we avoid the translation at -// a minimal cost (and we use C++;). - -#include +// a minimal cost (and we use C++ ;). +#include #include + #include +#include + +#include namespace CGAL { namespace Intersections { @@ -77,93 +80,101 @@ bool do_bbox_intersect(const typename K::Triangle_3& triangle, // if you do not know it, or if it does not exist, // use get_min_max without the AXE template parameter // available in _plane_is_cuboid_do_intersect.h -template +template inline -void get_min_max(const typename K::FT& px, - const typename K::FT& py, - const typename K::FT& pz, +void get_min_max(const FT& px, + const FT& py, + const FT& pz, const Box3& c, - typename K::Point_3& p_min, - typename K::Point_3& p_max) + std::array& p_min, + std::array& p_max) { - CGAL_USE(px); - CGAL_USE(py); - CGAL_USE(pz); - if(AXE == 0 || px > 0) { if(AXE == 1 || py > 0) { if(AXE == 2 || pz > 0) { - p_min = typename K::Point_3(c.xmin(), c.ymin(),c.zmin()); - p_max = typename K::Point_3(c.xmax(), c.ymax(),c.zmax()); + p_min = CGAL::make_array(c.xmin(), c.ymin(),c.zmin()); + p_max = CGAL::make_array(c.xmax(), c.ymax(),c.zmax()); } else { - p_min = typename K::Point_3(c.xmin(), c.ymin(),c.zmax()); - p_max = typename K::Point_3(c.xmax(), c.ymax(),c.zmin()); + p_min = CGAL::make_array(c.xmin(), c.ymin(),c.zmax()); + p_max = CGAL::make_array(c.xmax(), c.ymax(),c.zmin()); } } else { if(AXE == 2 || pz > 0) { - p_min = typename K::Point_3(c.xmin(), c.ymax(),c.zmin()); - p_max = typename K::Point_3(c.xmax(), c.ymin(),c.zmax()); + p_min = CGAL::make_array(c.xmin(), c.ymax(),c.zmin()); + p_max = CGAL::make_array(c.xmax(), c.ymin(),c.zmax()); } else { - p_min = typename K::Point_3(c.xmin(), c.ymax(),c.zmax()); - p_max = typename K::Point_3(c.xmax(), c.ymin(),c.zmin()); + p_min = CGAL::make_array(c.xmin(), c.ymax(),c.zmax()); + p_max = CGAL::make_array(c.xmax(), c.ymin(),c.zmin()); } } } else { if(AXE == 1 || py > 0) { if(AXE == 2 || pz > 0) { - p_min = typename K::Point_3(c.xmax(), c.ymin(),c.zmin()); - p_max = typename K::Point_3(c.xmin(), c.ymax(),c.zmax()); + p_min = CGAL::make_array(c.xmax(), c.ymin(),c.zmin()); + p_max = CGAL::make_array(c.xmin(), c.ymax(),c.zmax()); } else { - p_min = typename K::Point_3(c.xmax(), c.ymin(),c.zmax()); - p_max = typename K::Point_3(c.xmin(), c.ymax(),c.zmin()); + p_min = CGAL::make_array(c.xmax(), c.ymin(),c.zmax()); + p_max = CGAL::make_array(c.xmin(), c.ymax(),c.zmin()); } } else { if(AXE == 2 || pz > 0) { - p_min = typename K::Point_3(c.xmax(), c.ymax(),c.zmin()); - p_max = typename K::Point_3(c.xmin(), c.ymin(),c.zmax()); + p_min = CGAL::make_array(c.xmax(), c.ymax(),c.zmin()); + p_max = CGAL::make_array(c.xmin(), c.ymin(),c.zmax()); } else { - p_min = typename K::Point_3(c.xmax(), c.ymax(),c.zmax()); - p_max = typename K::Point_3(c.xmin(), c.ymin(),c.zmin()); + p_min = CGAL::make_array(c.xmax(), c.ymax(),c.zmax()); + p_max = CGAL::make_array(c.xmin(), c.ymin(),c.zmin()); } } } } -template +template inline -typename K::FT -do_axis_intersect_aux(const typename K::FT& alpha, - const typename K::FT& beta, - const typename K::Vector_3* sides) +Uncertain +do_axis_intersect_aux_A0(const FT& alpha, + const FT& beta, + const std::array, 3>& sides, + Fct do_axis_intersect_aux_impl) { - switch ( AXE ) - { - case 0: - return -sides[SIDE].z()*alpha + sides[SIDE].y()*beta; - case 1: - return sides[SIDE].z()*alpha - sides[SIDE].x()*beta; - case 2: - return -sides[SIDE].y()*alpha + sides[SIDE].x()*beta; - default: - CGAL_error(); - return typename K::FT(0); - } + return do_axis_intersect_aux_impl(alpha, beta, sides[SIDE][2], sides[SIDE][1]); +} + +template +inline +Uncertain +do_axis_intersect_aux_A1(const FT& alpha, + const FT& beta, + const std::array, 3>& sides, + Fct do_axis_intersect_aux_impl) +{ + return do_axis_intersect_aux_impl(beta, alpha, sides[SIDE][0], sides[SIDE][2]); +} + +template +inline +Uncertain +do_axis_intersect_aux_A2(const FT& alpha, + const FT& beta, + const std::array, 3>& sides, + Fct do_axis_intersect_aux_impl) +{ + return do_axis_intersect_aux_impl(alpha, beta, sides[SIDE][1], sides[SIDE][0]); } //given a vector checks whether it is collinear to a base vector //of the orthonormal frame. return -1 otherwise -template +template inline int -collinear_axis(typename K::Vector_3 side) +collinear_axis(const std::array& side) { if(certainly(side[0] == 0)) { @@ -177,201 +188,243 @@ collinear_axis(typename K::Vector_3 side) if(certainly(side[1] == 0) && certainly(side[2] == 0)) return 0; } - return -1; } -template +template inline -Uncertain do_axis_intersect(const typename K::Triangle_3& triangle, - const typename K::Vector_3* sides, - const Box3& bbox) +Uncertain do_axis_intersect(const std::array, 3>& triangle, + const std::array, 3>& sides, + const Box3& bbox, + Fct do_axis_intersect_aux_impl) { - const typename K::Point_3* j = & triangle.vertex(SIDE); - const typename K::Point_3* k = & triangle.vertex((SIDE+2)%3); - - typename K::Point_3 p_min, p_max; - get_min_max(AXE==0? 0: AXE==1? sides[SIDE].z(): -sides[SIDE].y(), - AXE==0? -sides[SIDE].z(): AXE==1? 0: sides[SIDE].x(), - AXE==0? sides[SIDE].y(): AXE==1? -sides[SIDE].x(): - typename K::FT(0), bbox, p_min, p_max); - - switch ( AXE ) + const std::array& j = triangle[SIDE]; + const std::array& k = triangle[(SIDE+2)%3]; + const std::array* j_ptr = &j; + const std::array* k_ptr = &k; + + std::array p_min, p_max; + get_min_max(AXE==0? 0: AXE==1? sides[SIDE][2]: -sides[SIDE][1], + AXE==0? -sides[SIDE][2]: AXE==1? 0: sides[SIDE][0], + AXE==0? sides[SIDE][1]: AXE==1? -sides[SIDE][0]: + FT(0), bbox, p_min, p_max); + + switch(AXE) { - case 0: { + case 0: + { // t_max >= t_min - Uncertain b = do_axis_intersect_aux(k->y()-j->y(), k->z()-j->z(), sides) >= 0; - if (is_indeterminate(b)) + Uncertain b = do_axis_intersect_aux_A0(k[1]-j[1], k[2]-j[2], sides, do_axis_intersect_aux_impl) != NEGATIVE; + if(is_indeterminate(b)) return b; - if(b) std::swap(j,k); - return CGAL_AND((do_axis_intersect_aux(p_min.y()-j->y(), p_min.z()-j->z(), sides) <= 0), - (do_axis_intersect_aux(p_max.y()-k->y(), p_max.z()-k->z(), sides) >= 0) ); + if(b) + std::swap(j_ptr, k_ptr); + + return CGAL_AND((do_axis_intersect_aux_A0(p_min[1]-(*j_ptr)[1], p_min[2]-(*j_ptr)[2], + sides, do_axis_intersect_aux_impl) != POSITIVE), + (do_axis_intersect_aux_A0(p_max[1]-(*k_ptr)[1], p_max[2]-(*k_ptr)[2], + sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } - case 1: { + case 1: + { // t_max >= t_min - Uncertain b = do_axis_intersect_aux(k->x()-j->x(), k->z()-j->z(), sides) >= 0; - if (is_indeterminate(b)) + Uncertain b = do_axis_intersect_aux_A1(k[0]-j[0], k[2]-j[2], sides, do_axis_intersect_aux_impl) != NEGATIVE; + if(is_indeterminate(b)) return b; - if(b) std::swap(j,k); - return CGAL_AND((do_axis_intersect_aux(p_min.x()-j->x(), p_min.z()-j->z(), sides) <= 0), - (do_axis_intersect_aux(p_max.x()-k->x(), p_max.z()-k->z(), sides) >= 0) ); + if(b) + std::swap(j_ptr, k_ptr); + return CGAL_AND((do_axis_intersect_aux_A1(p_min[0]-(*j_ptr)[0], p_min[2]-(*j_ptr)[2], + sides, do_axis_intersect_aux_impl) != POSITIVE), + (do_axis_intersect_aux_A1(p_max[0]-(*k_ptr)[0], p_max[2]-(*k_ptr)[2], + sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } - case 2: { + case 2: + { // t_max >= t_min - Uncertain b = do_axis_intersect_aux(k->x()-j->x(), k->y()-j->y(), sides) >= 0; - if ( is_indeterminate(b)) + Uncertain b = do_axis_intersect_aux_A2(k[0]-j[0], k[1]-j[1], sides, do_axis_intersect_aux_impl) != NEGATIVE; + if(is_indeterminate(b)) return b; - if(b) std::swap(j,k); - return CGAL_AND((do_axis_intersect_aux(p_min.x()-j->x(), p_min.y()-j->y(), sides) <= 0), - (do_axis_intersect_aux(p_max.x()-k->x(), p_max.y()-k->y(), sides) >= 0) ); + if(b) + std::swap(j_ptr, k_ptr); + + return CGAL_AND((do_axis_intersect_aux_A2(p_min[0]-(*j_ptr)[0], p_min[1]-(*j_ptr)[1], + sides, do_axis_intersect_aux_impl) != POSITIVE), + (do_axis_intersect_aux_A2(p_max[0]-(*k_ptr)[0], p_max[1]-(*k_ptr)[1], + sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } - default: - // Should not happen + default: // Should not happen CGAL_error(); return make_uncertain(false); } } -template -bool do_intersect_bbox_or_iso_cuboid(const typename K::Triangle_3& a_triangle, - const Box3& a_bbox, - const K& k) +template +Uncertain +do_intersect_bbox_or_iso_cuboid_impl(const std::array< std::array, 3>& triangle, + const Box3& bbox, + Fct do_axis_intersect_aux_impl) { - - if(certainly_not( do_bbox_intersect(a_triangle, a_bbox) )) - return false; - - if(certainly_not( do_intersect(a_triangle.supporting_plane(), a_bbox, k) )) - return false; - - typename K::Vector_3 sides[3]; - sides[0] = a_triangle[1] - a_triangle[0]; - sides[1] = a_triangle[2] - a_triangle[1]; - sides[2] = a_triangle[0] - a_triangle[2]; - int forbidden_axis=-1; - int forbidden_size=-1; + std::array< std::array, 3> sides = + {{ + {triangle[1][0] - triangle[0][0], triangle[1][1] - triangle[0][1], triangle[1][2] - triangle[0][2]}, + {triangle[2][0] - triangle[1][0], triangle[2][1] - triangle[1][1], triangle[2][2] - triangle[1][2]}, + {triangle[0][0] - triangle[2][0], triangle[0][1] - triangle[2][1], triangle[0][2] - triangle[2][2]} + }}; + + int forbidden_axis = -1; + int forbidden_size = -1; //determine whether one vector is collinear with an axis - int tmp=collinear_axis(sides[0]); - if ( tmp!= -1){ - forbidden_axis=tmp; - forbidden_size=0; + int tmp = collinear_axis(sides[0]); + if(tmp != -1) + { + forbidden_axis = tmp; + forbidden_size = 0; } - else{ - tmp=collinear_axis(sides[1]); - if ( tmp!= -1){ - forbidden_axis=tmp; - forbidden_size=1; + else + { + tmp = collinear_axis(sides[1]); + if(tmp != -1) + { + forbidden_axis = tmp; + forbidden_size = 1; } - else{ - tmp=collinear_axis(sides[2]); - if ( tmp!= -1){ - forbidden_axis=tmp; - forbidden_size=2; + else + { + tmp = collinear_axis(sides[2]); + if(tmp != -1) + { + forbidden_axis = tmp; + forbidden_size = 2; } } } -#if 0 - typename K::Point_3 p(a_bbox.xmin(), a_bbox.ymin(), a_bbox.zmin()); - typename K::Point_3 q(a_bbox.xmax(), a_bbox.ymax(), a_bbox.zmax()); - - typename K::Point_3 m = CGAL::midpoint(p,q); - typename K::Vector_3 v = m - CGAL::ORIGIN; - - typename K::Triangle_3 triangle(a_triangle[0]-v, a_triangle[1]-v, a_triangle[2]-v); - - Box3 bbox( (p-v).bbox() + (q-v).bbox()); - -#else - const typename K::Triangle_3& triangle = a_triangle; - const Box3& bbox = a_bbox; -#endif - // Create a "certainly true" Uncertain ind_or_true = make_uncertain(true); - if (forbidden_axis!=0){ - if (forbidden_size!=0){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); - if(is_indeterminate(b)){ + if(forbidden_axis != 0) + { + if(forbidden_size != 0) + { + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); + if(is_indeterminate(b)) ind_or_true = b; - } else if(! b){ + else if(!b) return false; - } } - if (forbidden_size!=1){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); - if(is_indeterminate(b)){ + + if(forbidden_size != 1) + { + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); + if(is_indeterminate(b)) ind_or_true = b; - } else if(! b){ + else if(!b) return false; - } } - if (forbidden_size!=2){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); - if(is_indeterminate(b)){ + + if(forbidden_size != 2) + { + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); + if(is_indeterminate(b)) ind_or_true = b; - } else if(! b){ + else if(!b) return false; - } } } - if (forbidden_axis!=1){ - if (forbidden_size!=0){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); - if(is_indeterminate(b)){ + if(forbidden_axis != 1) + { + if(forbidden_size != 0) + { + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); + if(is_indeterminate(b)) ind_or_true = b; - } else if(! b){ + else if(!b) return false; - } } - if (forbidden_size!=1){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); - if(is_indeterminate(b)){ + + if(forbidden_size != 1) + { + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); + if(is_indeterminate(b)) ind_or_true = b; - } else if(! b){ + else if(!b) return false; - } } - if (forbidden_size!=2){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); - if(is_indeterminate(b)){ + + if(forbidden_size != 2) + { + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); + if(is_indeterminate(b)) ind_or_true = b; - } else if(! b){ + else if(!b) return false; - } } } - if (forbidden_axis!=2){ - if (forbidden_size!=0){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); - if(is_indeterminate(b)){ + if(forbidden_axis != 2) + { + if(forbidden_size != 0) + { + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); + if(is_indeterminate(b)) ind_or_true = b; - } else if(! b){ + else if(!b) return false; - } } - if (forbidden_size!=1){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); - if(is_indeterminate(b)){ + + if(forbidden_size != 1) + { + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); + if(is_indeterminate(b)) ind_or_true = b; - } else if(! b){ + else if(!b) return false; - } } - if (forbidden_size!=2){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); - if(is_indeterminate(b)){ + + if(forbidden_size != 2) + { + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); + if(is_indeterminate(b)) ind_or_true = b; - } else if(! b){ + else if(!b) return false; - } } } - return ind_or_true; // throws exception in case it is indeterminate + + return ind_or_true; +} + +template +bool do_intersect_bbox_or_iso_cuboid(const typename K::Triangle_3& a_triangle, + const Box3& a_bbox, + const K& k) +{ + if(certainly_not(do_bbox_intersect(a_triangle, a_bbox))) + return false; + + if(certainly_not(do_intersect(a_triangle.supporting_plane(), a_bbox, k))) + return false; + + typedef typename K::FT FT; + auto do_axis_intersect_aux_impl = [](const FT& alpha, + const FT& beta, + const FT& c_alpha, + const FT& c_beta) -> Uncertain + { + return CGAL::sign(- c_alpha * alpha + c_beta * beta); + }; + + std::array< std::array, 3> triangle = + {{ + { a_triangle[0][0], a_triangle[0][1], a_triangle[0][2] }, + { a_triangle[1][0], a_triangle[1][1], a_triangle[1][2] }, + { a_triangle[2][0], a_triangle[2][1], a_triangle[2][2] } + }}; + + // exception will be thrown in case the output is indeterminate + return do_intersect_bbox_or_iso_cuboid_impl(triangle, a_bbox, do_axis_intersect_aux_impl); } template