From 17c132640c8fe251a31e8bdf73813ecd49911a90 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Mon, 7 Jun 2021 13:06:47 +0200 Subject: [PATCH 01/30] Filtered_kernel: Do_intersect(Tetrahedron_3,Bbox_3) --- .../internal/Static_filters/Do_intersect_3.h | 96 +++++++++++++++++++ 1 file changed, 96 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 8a5c95b67854..5e45bf5302fb 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 @@ -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; @@ -128,6 +129,101 @@ class Do_intersect_3 + result_type + operator()(const Tetrahedron_3 &t, const Bbox_3& b) const + { + Get_approx get_approx; + bool in = false, out = false; + double px, py, pz; + + { + const Point_3& p = t[0]; + 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) ) + { + if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ + in = true; + }else{ + out = true; + } + }else{ + return Base::operator()(t,b); + } + } + + { + const Point_3& p = t[1]; + 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) ) + { + if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ + if(out){ + return true; + } + in = true; + }else{ + if(in){ + return true; + } + out = true; + } + }else{ + return Base::operator()(t,b); + } + } + + { + const Point_3& p = t[2]; + 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) ) + { + if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ + if(out){ + return true; + } + in = true; + }else{ + if(in){ + return true; + } + out = true; + } + }else{ + return Base::operator()(t,b); + } + } + + + + { + const Point_3& p = t[3]; + 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) ) + { + if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ + if(out){ + return true; + } + in = true; + }else{ + if(in){ + return true; + } + out = true; + } + }else{ + return Base::operator()(t,b); + } + } + + + if(in){ + return true; + } + return Base::operator()(t,b); + } + + result_type operator()(const Bbox_3& b, const Ray_3 &r) const { From b3cfabfe3639205f9aa7c8ff7826726e718e0793 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Mon, 7 Jun 2021 14:04:56 +0200 Subject: [PATCH 02/30] In fact we can return true as soon as one tet vertex is in the bbox --- .../internal/Static_filters/Do_intersect_3.h | 37 +------------------ 1 file changed, 2 insertions(+), 35 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 5e45bf5302fb..7b16659b5c67 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 @@ -133,7 +133,6 @@ class Do_intersect_3 operator()(const Tetrahedron_3 &t, const Bbox_3& b) const { Get_approx get_approx; - bool in = false, out = false; double px, py, pz; { @@ -142,9 +141,7 @@ class Do_intersect_3 fit_in_double(get_approx(p).z(), pz) ) { if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - in = true; - }else{ - out = true; + return true; } }else{ return Base::operator()(t,b); @@ -157,15 +154,7 @@ class Do_intersect_3 fit_in_double(get_approx(p).z(), pz) ) { if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - if(out){ - return true; - } - in = true; - }else{ - if(in){ return true; - } - out = true; } }else{ return Base::operator()(t,b); @@ -178,15 +167,7 @@ class Do_intersect_3 fit_in_double(get_approx(p).z(), pz) ) { if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - if(out){ return true; - } - in = true; - }else{ - if(in){ - return true; - } - out = true; } }else{ return Base::operator()(t,b); @@ -201,25 +182,11 @@ class Do_intersect_3 fit_in_double(get_approx(p).z(), pz) ) { if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - if(out){ - return true; - } - in = true; - }else{ - if(in){ return true; - } - out = true; } - }else{ - return Base::operator()(t,b); - } + } } - - if(in){ - return true; - } return Base::operator()(t,b); } From 2489e16220ad7a45021ea27074b7179f6c058082 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Mon, 7 Jun 2021 15:00:05 +0200 Subject: [PATCH 03/30] Add check that bboxes overlap --- .../include/CGAL/internal/Static_filters/Do_intersect_3.h | 3 +++ 1 file changed, 3 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 7b16659b5c67..69ec9a615cd4 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 @@ -135,6 +135,9 @@ class Do_intersect_3 Get_approx get_approx; double px, py, pz; + if(! do_overlap(t.bbox(),b)){ + return false; + } { const Point_3& p = t[0]; if (fit_in_double(get_approx(p).x(), px) && fit_in_double(get_approx(p).y(), py) && From ce0037c37918469b1bfb82bbeddf96812fca2784 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 11 Jun 2021 13:35:17 +0200 Subject: [PATCH 04/30] Add the same for Triangle_3 --- AABB_tree/include/CGAL/AABB_traits.h | 95 ++++++++++++++++++- .../internal/Static_filters/Do_intersect_3.h | 58 +++++++++++ 2 files changed, 150 insertions(+), 3 deletions(-) diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index e3ffcacf50ce..2ca991bcbab1 100644 --- a/AABB_tree/include/CGAL/AABB_traits.h +++ b/AABB_tree/include/CGAL/AABB_traits.h @@ -33,6 +33,23 @@ namespace CGAL { +template +struct Tet { + typedef typename Kernel_traits::Kernel::Point_3 Point_3; + typedef typename Kernel_traits::Kernel::Triangle_3 Triangle_3; + + Tet() + {} + + Tet(const T& tet, bool b0, bool b1, bool b2, bool b3) + : tet(tet), bbox(tet.bbox()), b0(b0), b1(b1), b2(b2), b3(b3) + {} + + T tet; + Bbox_3 bbox; + bool b0=false, b1=false, b2=false, b3=false; +}; + namespace internal{ namespace AABB_tree { template @@ -324,19 +341,91 @@ class AABB_traits /// In the case the query is a `CGAL::AABB_tree`, the `do_intersect()` /// function of this tree is used. class Do_intersect { + typedef typename GeomTraits::Ray_3 Ray_3; const AABB_traits& m_traits; public: Do_intersect(const AABB_traits& traits) :m_traits(traits) {} template - bool operator()(const Query& q, const Bounding_box& bbox) const + bool operator()(const Tet& q, const Bounding_box& bbox) const { - return CGAL::do_intersect(q, bbox); + // bool expected = CGAL::do_intersect(q.tet, bbox); + // return expected; + bool retur = false; + if(! do_overlap(q.bbox, bbox)){ + retur = false; + // if(retur != expected){ std::cout << "oops 1" << std::endl; } + return retur; + } + typedef typename Tet::Triangle_3 Triangle_3; + if(! q.b0){ + Triangle_3 t(q.tet[1], q.tet[2],q.tet[3]); + if(do_intersect(t,bbox)){ + retur = true; + // if(retur != expected){ std::cout << "oops 2" << std::endl; } + return retur; + } + } + if(! q.b1){ + Triangle_3 t(q.tet[2], q.tet[3],q.tet[0]); + if(do_intersect(t,bbox)){ + retur = true; + // if(retur != expected){ std::cout << "oops 3" << std::endl; } + return retur; + } + } + if(! q.b2){ + Triangle_3 t(q.tet[3], q.tet[0],q.tet[1]); + if(do_intersect(t,bbox)){ + retur = true; + // if(retur != expected){ std::cout << "oops 4" << std::endl; } + return retur; + } + } + if(! q.b3){ + Triangle_3 t(q.tet[0], q.tet[1],q.tet[2]); + if(do_intersect(t,bbox)){ + retur = true; + // if(retur != expected){ std::cout << "oops 5" << std::endl; } + return retur; + } + } + Point_3 bbp(bbox.xmin(),bbox.ymin(),bbox.zmin()); + if(q.tet.has_on_bounded_side(bbp)){ + retur = true; + // if(retur != expected){ std::cout << "oops 6" << std::endl; } + return retur; + } + + return false; } template - bool operator()(const Query& q, const Primitive& pr) const + bool operator()(const Tet& q, const Primitive& pr) const + { + return GeomTraits().do_intersect_3_object()(q.tet, internal::Primitive_helper::get_datum(pr,m_traits)); + } + + + bool operator()(const Ray_3& q, const Bounding_box& bbox) const + { + return CGAL::do_intersect(q, bbox); + } + + + bool operator()(const Ray_3& q, const Primitive& pr) const + { + return GeomTraits().do_intersect_3_object()(q, internal::Primitive_helper::get_datum(pr,m_traits)); + } + + bool operator()(const Sphere_3& q, const Bounding_box& bbox) const + { + return CGAL::do_intersect(q, bbox); + } + + + bool operator()(const Sphere_3& q, const Primitive& pr) const { return GeomTraits().do_intersect_3_object()(q, internal::Primitive_helper::get_datum(pr,m_traits)); } 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 69ec9a615cd4..1febbafac2b6 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 @@ -235,6 +235,64 @@ 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); + } + + + result_type + operator()(const Triangle_3 &t, const Bbox_3& b) const + { + Get_approx get_approx; + double px, py, pz; + + if(! do_overlap(t.bbox(),b)){ + return false; + } + { + const Point_3& p = t[0]; + 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) ) + { + if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ + return true; + } + }else{ + return Base::operator()(t,b); + } + } + + { + const Point_3& p = t[1]; + 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) ) + { + if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ + return true; + } + }else{ + return Base::operator()(t,b); + } + } + + { + const Point_3& p = t[2]; + 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) ) + { + if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ + return true; + } + }else{ + return Base::operator()(t,b); + } + } + return Base::operator()(t,b); + } + result_type operator()(const Bbox_3& b, const Sphere_3 &s) const { From e5c72e5bedfae9e8b6dfc6575c3a55cba23ce148 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 11 Jun 2021 14:20:38 +0200 Subject: [PATCH 05/30] Undo accidental change --- AABB_tree/include/CGAL/AABB_traits.h | 95 +--------------------------- 1 file changed, 3 insertions(+), 92 deletions(-) diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index 2ca991bcbab1..e3ffcacf50ce 100644 --- a/AABB_tree/include/CGAL/AABB_traits.h +++ b/AABB_tree/include/CGAL/AABB_traits.h @@ -33,23 +33,6 @@ namespace CGAL { -template -struct Tet { - typedef typename Kernel_traits::Kernel::Point_3 Point_3; - typedef typename Kernel_traits::Kernel::Triangle_3 Triangle_3; - - Tet() - {} - - Tet(const T& tet, bool b0, bool b1, bool b2, bool b3) - : tet(tet), bbox(tet.bbox()), b0(b0), b1(b1), b2(b2), b3(b3) - {} - - T tet; - Bbox_3 bbox; - bool b0=false, b1=false, b2=false, b3=false; -}; - namespace internal{ namespace AABB_tree { template @@ -341,91 +324,19 @@ class AABB_traits /// In the case the query is a `CGAL::AABB_tree`, the `do_intersect()` /// function of this tree is used. class Do_intersect { - typedef typename GeomTraits::Ray_3 Ray_3; const AABB_traits& m_traits; public: Do_intersect(const AABB_traits& traits) :m_traits(traits) {} template - bool operator()(const Tet& q, const Bounding_box& bbox) const - { - // bool expected = CGAL::do_intersect(q.tet, bbox); - // return expected; - bool retur = false; - if(! do_overlap(q.bbox, bbox)){ - retur = false; - // if(retur != expected){ std::cout << "oops 1" << std::endl; } - return retur; - } - typedef typename Tet::Triangle_3 Triangle_3; - if(! q.b0){ - Triangle_3 t(q.tet[1], q.tet[2],q.tet[3]); - if(do_intersect(t,bbox)){ - retur = true; - // if(retur != expected){ std::cout << "oops 2" << std::endl; } - return retur; - } - } - if(! q.b1){ - Triangle_3 t(q.tet[2], q.tet[3],q.tet[0]); - if(do_intersect(t,bbox)){ - retur = true; - // if(retur != expected){ std::cout << "oops 3" << std::endl; } - return retur; - } - } - if(! q.b2){ - Triangle_3 t(q.tet[3], q.tet[0],q.tet[1]); - if(do_intersect(t,bbox)){ - retur = true; - // if(retur != expected){ std::cout << "oops 4" << std::endl; } - return retur; - } - } - if(! q.b3){ - Triangle_3 t(q.tet[0], q.tet[1],q.tet[2]); - if(do_intersect(t,bbox)){ - retur = true; - // if(retur != expected){ std::cout << "oops 5" << std::endl; } - return retur; - } - } - Point_3 bbp(bbox.xmin(),bbox.ymin(),bbox.zmin()); - if(q.tet.has_on_bounded_side(bbp)){ - retur = true; - // if(retur != expected){ std::cout << "oops 6" << std::endl; } - return retur; - } - - return false; - } - - template - bool operator()(const Tet& q, const Primitive& pr) const - { - return GeomTraits().do_intersect_3_object()(q.tet, internal::Primitive_helper::get_datum(pr,m_traits)); - } - - - bool operator()(const Ray_3& q, const Bounding_box& bbox) const - { - return CGAL::do_intersect(q, bbox); - } - - - bool operator()(const Ray_3& q, const Primitive& pr) const - { - return GeomTraits().do_intersect_3_object()(q, internal::Primitive_helper::get_datum(pr,m_traits)); - } - - bool operator()(const Sphere_3& q, const Bounding_box& bbox) const + bool operator()(const Query& q, const Bounding_box& bbox) const { return CGAL::do_intersect(q, bbox); } - - bool operator()(const Sphere_3& q, const Primitive& pr) const + template + bool operator()(const Query& q, const Primitive& pr) const { return GeomTraits().do_intersect_3_object()(q, internal::Primitive_helper::get_datum(pr,m_traits)); } From 46b112dd544df390b138c7041cc97b5699e8e160 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 15 Jun 2021 13:22:34 +0200 Subject: [PATCH 06/30] split the function to avoid the switch --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 59 +++++++++++-------- 1 file changed, 34 insertions(+), 25 deletions(-) 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 a46ee15c7ebd..89e3a497ca71 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 @@ -141,25 +141,34 @@ namespace internal { } - 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) + do_axis_intersect_aux_A0(const typename K::FT& alpha, + const typename K::FT& beta, + const typename K::Vector_3* sides) { - 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 -sides[SIDE].z()*alpha + sides[SIDE].y()*beta; + } + + template + inline + typename K::FT + do_axis_intersect_aux_A1(const typename K::FT& alpha, + const typename K::FT& beta, + const typename K::Vector_3* sides) + { + return sides[SIDE].z()*alpha - sides[SIDE].x()*beta; + } + + template + inline + typename K::FT + do_axis_intersect_aux_A2(const typename K::FT& alpha, + const typename K::FT& beta, + const typename K::Vector_3* sides) + { + return -sides[SIDE].y()*alpha + sides[SIDE].x()*beta; } //given a vector checks whether it is collinear to a base vector @@ -198,31 +207,31 @@ namespace internal { { case 0: { // t_max >= t_min - Uncertain b = do_axis_intersect_aux(k->y()-j->y(), k->z()-j->z(), sides) >= 0; + Uncertain b = do_axis_intersect_aux_A0(k->y()-j->y(), k->z()-j->z(), sides) >= 0; 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) ); + return CGAL_AND((do_axis_intersect_aux_A0(p_min.y()-j->y(), p_min.z()-j->z(), sides) <= 0), + (do_axis_intersect_aux_A0(p_max.y()-k->y(), p_max.z()-k->z(), sides) >= 0) ); } case 1: { // t_max >= t_min - Uncertain b = do_axis_intersect_aux(k->x()-j->x(), k->z()-j->z(), sides) >= 0; + Uncertain b = do_axis_intersect_aux_A1(k->x()-j->x(), k->z()-j->z(), sides) >= 0; 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) ); + return CGAL_AND((do_axis_intersect_aux_A1(p_min.x()-j->x(), p_min.z()-j->z(), sides) <= 0), + (do_axis_intersect_aux_A1(p_max.x()-k->x(), p_max.z()-k->z(), sides) >= 0) ); } case 2: { // t_max >= t_min - Uncertain b = do_axis_intersect_aux(k->x()-j->x(), k->y()-j->y(), sides) >= 0; + Uncertain b = do_axis_intersect_aux_A2(k->x()-j->x(), k->y()-j->y(), sides) >= 0; 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) ); + return CGAL_AND((do_axis_intersect_aux_A2(p_min.x()-j->x(), p_min.y()-j->y(), sides) <= 0), + (do_axis_intersect_aux_A2(p_max.x()-k->x(), p_max.y()-k->y(), sides) >= 0) ); } default: // Should not happen From 65d630536b521d41af913a9f0ebb26b9e4fc0197 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 15 Jun 2021 14:36:26 +0200 Subject: [PATCH 07/30] factorize --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) 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 89e3a497ca71..5245c40ff445 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 @@ -140,6 +140,16 @@ namespace internal { } } + template + inline + Sign + do_axis_intersect_aux_impl(const typename K::FT& alpha, + const typename K::FT& beta, + const typename K::FT& c_alpha, + const typename K::FT& c_beta) + { + return - c_alpha * alpha + c_beta * beta; + } template inline @@ -148,7 +158,7 @@ namespace internal { const typename K::FT& beta, const typename K::Vector_3* sides) { - return -sides[SIDE].z()*alpha + sides[SIDE].y()*beta; + return do_axis_intersect_aux_impl(alpha, beta, sides[SIDE].z(), sides[SIDE].y()); } template @@ -158,7 +168,7 @@ namespace internal { const typename K::FT& beta, const typename K::Vector_3* sides) { - return sides[SIDE].z()*alpha - sides[SIDE].x()*beta; + return do_axis_intersect_aux_impl(beta, alpha, sides[SIDE].x(), sides[SIDE].z()); } template @@ -168,7 +178,7 @@ namespace internal { const typename K::FT& beta, const typename K::Vector_3* sides) { - return -sides[SIDE].y()*alpha + sides[SIDE].x()*beta; + return do_axis_intersect_aux_impl(alpha, beta, sides[SIDE].y(), sides[SIDE].x()); } //given a vector checks whether it is collinear to a base vector From 93476f0b43a2a3ef311d9574cb77bcf65e8a1ffb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 15 Jun 2021 15:03:02 +0200 Subject: [PATCH 08/30] use sign --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) 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 5245c40ff445..9de6856dcaad 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 @@ -142,18 +142,18 @@ namespace internal { template inline - Sign + Uncertain do_axis_intersect_aux_impl(const typename K::FT& alpha, const typename K::FT& beta, const typename K::FT& c_alpha, const typename K::FT& c_beta) { - return - c_alpha * alpha + c_beta * beta; + return sign(- c_alpha * alpha + c_beta * beta); } template inline - typename K::FT + Sign do_axis_intersect_aux_A0(const typename K::FT& alpha, const typename K::FT& beta, const typename K::Vector_3* sides) @@ -163,7 +163,7 @@ namespace internal { template inline - typename K::FT + Sign do_axis_intersect_aux_A1(const typename K::FT& alpha, const typename K::FT& beta, const typename K::Vector_3* sides) @@ -173,7 +173,7 @@ namespace internal { template inline - typename K::FT + Sign do_axis_intersect_aux_A2(const typename K::FT& alpha, const typename K::FT& beta, const typename K::Vector_3* sides) @@ -217,31 +217,31 @@ namespace internal { { case 0: { // t_max >= t_min - Uncertain b = do_axis_intersect_aux_A0(k->y()-j->y(), k->z()-j->z(), sides) >= 0; + Uncertain b = do_axis_intersect_aux_A0(k->y()-j->y(), k->z()-j->z(), sides) != NEGATIVE; if (is_indeterminate(b)) return b; if(b) std::swap(j,k); - return CGAL_AND((do_axis_intersect_aux_A0(p_min.y()-j->y(), p_min.z()-j->z(), sides) <= 0), - (do_axis_intersect_aux_A0(p_max.y()-k->y(), p_max.z()-k->z(), sides) >= 0) ); + return CGAL_AND((do_axis_intersect_aux_A0(p_min.y()-j->y(), p_min.z()-j->z(), sides) != POSITIVE), + (do_axis_intersect_aux_A0(p_max.y()-k->y(), p_max.z()-k->z(), sides) != NEGATIVE) ); } case 1: { // t_max >= t_min - Uncertain b = do_axis_intersect_aux_A1(k->x()-j->x(), k->z()-j->z(), sides) >= 0; + Uncertain b = do_axis_intersect_aux_A1(k->x()-j->x(), k->z()-j->z(), sides) != NEGATIVE; if (is_indeterminate(b)) return b; if(b) std::swap(j,k); - return CGAL_AND((do_axis_intersect_aux_A1(p_min.x()-j->x(), p_min.z()-j->z(), sides) <= 0), - (do_axis_intersect_aux_A1(p_max.x()-k->x(), p_max.z()-k->z(), sides) >= 0) ); + return CGAL_AND((do_axis_intersect_aux_A1(p_min.x()-j->x(), p_min.z()-j->z(), sides) != POSITIVE), + (do_axis_intersect_aux_A1(p_max.x()-k->x(), p_max.z()-k->z(), sides) != NEGATIVE) ); } case 2: { // t_max >= t_min - Uncertain b = do_axis_intersect_aux_A2(k->x()-j->x(), k->y()-j->y(), sides) >= 0; + Uncertain b = do_axis_intersect_aux_A2(k->x()-j->x(), k->y()-j->y(), sides) != NEGATIVE; if ( is_indeterminate(b)) return b; if(b) std::swap(j,k); - return CGAL_AND((do_axis_intersect_aux_A2(p_min.x()-j->x(), p_min.y()-j->y(), sides) <= 0), - (do_axis_intersect_aux_A2(p_max.x()-k->x(), p_max.y()-k->y(), sides) >= 0) ); + return CGAL_AND((do_axis_intersect_aux_A2(p_min.x()-j->x(), p_min.y()-j->y(), sides) != POSITIVE), + (do_axis_intersect_aux_A2(p_max.x()-k->x(), p_max.y()-k->y(), sides) != NEGATIVE) ); } default: // Should not happen From 6fe3eef000f56638f9a8056cdaf0f4aba9e74387 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 15 Jun 2021 15:30:22 +0200 Subject: [PATCH 09/30] use a for-loop --- .../internal/Static_filters/Do_intersect_3.h | 39 +++++-------------- 1 file changed, 9 insertions(+), 30 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 1febbafac2b6..1cf0993cfd96 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 @@ -249,47 +249,26 @@ class Do_intersect_3 Get_approx get_approx; double px, py, pz; + // check overlaps between bboxes if(! do_overlap(t.bbox(),b)){ return false; } + for (int i=0; i<3; ++i) { - const Point_3& p = t[0]; + 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) ) - { - if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - return true; - } - }else{ - return Base::operator()(t,b); + { + if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ + return true; + } } - } - - { - const Point_3& p = t[1]; - 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) ) - { - if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - return true; - } - }else{ + else + { return Base::operator()(t,b); } } - { - const Point_3& p = t[2]; - 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) ) - { - if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - return true; - } - }else{ - return Base::operator()(t,b); - } - } return Base::operator()(t,b); } From 5f2c6fdbdffa0b92976b40b811deee293ab8f8ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 16 Jun 2021 11:27:21 +0200 Subject: [PATCH 10/30] refactor the code prior to plug static filters --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 231 +++++++++--------- 1 file changed, 117 insertions(+), 114 deletions(-) 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 9de6856dcaad..80dd430c21fd 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 @@ -82,111 +82,101 @@ namespace internal { // 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 - inline - Uncertain - do_axis_intersect_aux_impl(const typename K::FT& alpha, - const typename K::FT& beta, - const typename K::FT& c_alpha, - const typename K::FT& c_beta) - { - return sign(- c_alpha * alpha + c_beta * beta); - } - template + template inline Sign - do_axis_intersect_aux_A0(const typename K::FT& alpha, - const typename K::FT& beta, - const typename K::Vector_3* sides) + do_axis_intersect_aux_A0(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].z(), sides[SIDE].y()); + return do_axis_intersect_aux_impl(alpha, beta, sides[SIDE][2], sides[SIDE][1]); } - template + template inline Sign - do_axis_intersect_aux_A1(const typename K::FT& alpha, - const typename K::FT& beta, - const typename K::Vector_3* sides) + 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].x(), sides[SIDE].z()); + return do_axis_intersect_aux_impl(beta, alpha, sides[SIDE][0], sides[SIDE][2]); } - template + template inline Sign - do_axis_intersect_aux_A2(const typename K::FT& alpha, - const typename K::FT& beta, - const typename K::Vector_3* sides) + 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].y(), sides[SIDE].x()); + 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) ){ if ( certainly(side[1]==0) ) return 2; if ( certainly(side[2]==0) ) return 1; @@ -198,50 +188,52 @@ namespace internal { 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< std::array , 3>& triangle, + const std::array< 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); + // TODO: no longer a ptr, is it an issue for swap? + std::array j = triangle[SIDE]; + std::array k = triangle[(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); + 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: { // t_max >= t_min - Uncertain b = do_axis_intersect_aux_A0(k->y()-j->y(), k->z()-j->z(), sides) != NEGATIVE; + 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_A0(p_min.y()-j->y(), p_min.z()-j->z(), sides) != POSITIVE), - (do_axis_intersect_aux_A0(p_max.y()-k->y(), p_max.z()-k->z(), sides) != NEGATIVE) ); + return CGAL_AND((do_axis_intersect_aux_A0(p_min[1]-j[1], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), + (do_axis_intersect_aux_A0(p_max[1]-k[1], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } case 1: { // t_max >= t_min - Uncertain b = do_axis_intersect_aux_A1(k->x()-j->x(), k->z()-j->z(), sides) != NEGATIVE; + 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_A1(p_min.x()-j->x(), p_min.z()-j->z(), sides) != POSITIVE), - (do_axis_intersect_aux_A1(p_max.x()-k->x(), p_max.z()-k->z(), sides) != NEGATIVE) ); + return CGAL_AND((do_axis_intersect_aux_A1(p_min[0]-j[0], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), + (do_axis_intersect_aux_A1(p_max[0]-k[0], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } case 2: { // t_max >= t_min - Uncertain b = do_axis_intersect_aux_A2(k->x()-j->x(), k->y()-j->y(), sides) != NEGATIVE; + 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_A2(p_min.x()-j->x(), p_min.y()-j->y(), sides) != POSITIVE), - (do_axis_intersect_aux_A2(p_max.x()-k->x(), p_max.y()-k->y(), sides) != NEGATIVE) ); + return CGAL_AND((do_axis_intersect_aux_A2(p_min[0]-j[0], p_min[1]-j[1], sides, do_axis_intersect_aux_impl) != POSITIVE), + (do_axis_intersect_aux_A2(p_max[0]-k[0], p_max[1]-k[1], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } default: // Should not happen @@ -250,38 +242,35 @@ namespace internal { } } - template - bool do_intersect_bbox_or_iso_cuboid(const typename K::Triangle_3& a_triangle, - const Box3& a_bbox, - const K& k) + template + bool do_intersect_bbox_or_iso_cuboid_impl(const std::array< std::array, 3>& triangle, + const Box3& bbox, + const K& k, + Fct do_axis_intersect_aux_impl) { + 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]} + }}; - 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; //determine whether one vector is collinear with an axis - int tmp=collinear_axis(sides[0]); + int tmp=collinear_axis(sides[0]); if ( tmp!= -1){ forbidden_axis=tmp; forbidden_size=0; } else{ - tmp=collinear_axis(sides[1]); + tmp=collinear_axis(sides[1]); if ( tmp!= -1){ forbidden_axis=tmp; forbidden_size=1; } else{ - tmp=collinear_axis(sides[2]); + tmp=collinear_axis(sides[2]); if ( tmp!= -1){ forbidden_axis=tmp; forbidden_size=2; @@ -289,28 +278,12 @@ namespace internal { } } -#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); + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); if(is_indeterminate(b)){ ind_or_true = b; } else if(! b){ @@ -318,7 +291,7 @@ namespace internal { } } if (forbidden_size!=1){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); if(is_indeterminate(b)){ ind_or_true = b; } else if(! b){ @@ -326,7 +299,7 @@ namespace internal { } } if (forbidden_size!=2){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); if(is_indeterminate(b)){ ind_or_true = b; } else if(! b){ @@ -337,7 +310,7 @@ namespace internal { if (forbidden_axis!=1){ if (forbidden_size!=0){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); if(is_indeterminate(b)){ ind_or_true = b; } else if(! b){ @@ -345,7 +318,7 @@ namespace internal { } } if (forbidden_size!=1){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); if(is_indeterminate(b)){ ind_or_true = b; } else if(! b){ @@ -353,7 +326,7 @@ namespace internal { } } if (forbidden_size!=2){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); if(is_indeterminate(b)){ ind_or_true = b; } else if(! b){ @@ -364,7 +337,7 @@ namespace internal { if (forbidden_axis!=2){ if (forbidden_size!=0){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); if(is_indeterminate(b)){ ind_or_true = b; } else if(! b){ @@ -372,7 +345,7 @@ namespace internal { } } if (forbidden_size!=1){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); if(is_indeterminate(b)){ ind_or_true = b; } else if(! b){ @@ -380,7 +353,7 @@ namespace internal { } } if (forbidden_size!=2){ - Uncertain b = do_axis_intersect(triangle, sides, bbox); + Uncertain b = do_axis_intersect(triangle, sides, bbox, do_axis_intersect_aux_impl); if(is_indeterminate(b)){ ind_or_true = b; } else if(! b){ @@ -391,6 +364,36 @@ namespace internal { return ind_or_true; // throws exception in case it is indeterminate } + 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 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] } + }}; + + return do_intersect_bbox_or_iso_cuboid_impl(triangle, a_bbox, k, do_axis_intersect_aux_impl); + } + template bool do_intersect(const typename K::Triangle_3& triangle, const CGAL::Bbox_3& bbox, From cdaabd3d72b0f678a2967fb0452ba26969cb9f36 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 16 Jun 2021 13:56:08 +0200 Subject: [PATCH 11/30] static filter for do_intersect bbox/triangle --- .../internal/Static_filters/Do_intersect_3.h | 105 +++++++++++++++++- .../internal/Bbox_3_Triangle_3_do_intersect.h | 16 +-- 2 files changed, 107 insertions(+), 14 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 1cf0993cfd96..3250b398cea4 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 @@ -246,20 +246,22 @@ class Do_intersect_3 result_type operator()(const Triangle_3 &t, const Bbox_3& b) const { - Get_approx get_approx; - double px, py, pz; - // check overlaps between bboxes if(! do_overlap(t.bbox(),b)){ return false; } + + // 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(), px) && fit_in_double(get_approx(p).y(), py) && - fit_in_double(get_approx(p).z(), pz) ) + 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]) ) { - if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ + 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; } } @@ -269,6 +271,97 @@ class Do_intersect_3 } } + // 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 + { + Sign int_tmp_result; + double double_tmp_result; + double eps; + double_tmp_result = ((-c_alpha * alpha) + (c_beta * beta)); + double max1 = fabs(c_alpha); + if( (max1 < fabs(c_beta)) ) + { + max1 = fabs(c_beta); + } + double max2 = fabs(alpha); + if( (max2 < fabs(beta)) ) + { + max2 = fabs(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) ) + { + return Uncertain(); + } + else + { + if( (upper_bound_1 > 1.67597599124282389316e+153) ) + { + return Uncertain(); + } + 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 + { + return Uncertain(); + } + } + } + return int_tmp_result; + }; + + auto do_intersect_plane_bbox = [](const Triangle_3& t, const Bbox_3& bbox) + { + typename SFK::Orientation_3 orient = SFK().orientation_3_object(); + CGAL::Orientation side = orient(t[0], t[1],t[2], Point_3(bbox.xmin(), bbox.ymin(),bbox.zmin())); + if(side == COPLANAR) return true; + CGAL::Oriented_side s = orient(t[0], t[1],t[2], Point_3(bbox.xmax(), bbox.ymax(),bbox.zmax())); + if(s != side) return true; + s = orient(t[0], t[1],t[2], Point_3(bbox.xmin(), bbox.ymin(),bbox.zmax())); + if(s != side) return true; + s = orient(t[0], t[1],t[2], Point_3(bbox.xmax(), bbox.ymax(),bbox.zmin())); + if(s != side) return true; + s = orient(t[0], t[1],t[2], Point_3(bbox.xmin(), bbox.ymax(),bbox.zmin())); + if(s != side) return true; + s = orient(t[0], t[1],t[2], Point_3(bbox.xmax(), bbox.ymin(),bbox.zmax())); + if(s != side) return true; + s = orient(t[0], t[1],t[2], Point_3(bbox.xmin(), bbox.ymax(),bbox.zmax())); + if(s != side) return true; + s = orient(t[0], t[1],t[2], Point_3(bbox.xmax(), bbox.ymin(),bbox.zmin())); + if(s != side) return true; + return false; + }; + + if ( !do_intersect_plane_bbox(t,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); } 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 80dd430c21fd..6af9fe95a17c 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 @@ -242,11 +242,11 @@ namespace internal { } } - template - bool do_intersect_bbox_or_iso_cuboid_impl(const std::array< std::array, 3>& triangle, - const Box3& bbox, - const K& k, - Fct do_axis_intersect_aux_impl) + 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) { std::array< std::array, 3> sides = {{ @@ -361,7 +361,7 @@ namespace internal { } } } - return ind_or_true; // throws exception in case it is indeterminate + return ind_or_true; } template @@ -390,8 +390,8 @@ namespace internal { { a_triangle[1][0], a_triangle[1][1], a_triangle[1][2] }, { a_triangle[2][0], a_triangle[2][1], a_triangle[2][2] } }}; - - return do_intersect_bbox_or_iso_cuboid_impl(triangle, a_bbox, k, do_axis_intersect_aux_impl); + // 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 From 602b93e30f98185931dde99162f8da103bf410df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 16 Jun 2021 16:39:38 +0200 Subject: [PATCH 12/30] improve intersection predicate between supporting plane and bbox quick and dirty with copy/paste from other files --- .../internal/Static_filters/Do_intersect_3.h | 271 ++++++++++++++++-- 1 file changed, 248 insertions(+), 23 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 3250b398cea4..5850056ef433 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 @@ -242,6 +242,252 @@ class Do_intersect_3 return this->operator()(t, b); } + /// TODO: move in utils + Uncertain sign_of_minor(double px, double py, double qx, double qy, double rx, double ry) const + { + 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 = fabs(qx_px); + if( (max1 < fabs(rx_px)) ) + { + max1 = fabs(rx_px); + } + double max2 = fabs(ry_py); + if( (max2 < fabs(qy_py)) ) + { + max2 = fabs(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) ) + { + return Uncertain(); + } + else + { + if( (upper_bound_1 > 1.67597599124282389316e+153) ) + { + return Uncertain(); + } + 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 + { + return Uncertain(); + } + } + } + return int_tmp_result; + }; + + bool get_cross_product_sign(const std::array< std::array, 3 >& pts, std::array& signs) const + { + 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)) return false; + signs[0]=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)) return false; + signs[1]=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)) return false; + signs[2]=s; + return true; + } + + bool do_intersect_supporting_plane_bbox(const Triangle_3& t, const std::array< std::array, 3>& pts, const Bbox_3& bbox) const + { + auto my_orient = [](const std::array< std::array, 3>& pts, double x, double y, double z) -> Uncertain + { + 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; + } + return Uncertain(); + }; + + auto orient = [my_orient](const Triangle_3& t, const std::array< std::array, 3>& pts, + double x, double y, double z) -> Orientation + { + Uncertain res = my_orient(pts, x, y, z); + if (!is_indeterminate(res)) return make_certain(res); + typename SFK::Orientation_3 orient = SFK().orientation_3_object(); + //TODO: avoid calling the static filter again... + 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 + 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 ! (orient(t, pts, p_max[0], p_max[1], p_max[2]) == ON_NEGATIVE_SIDE || + orient(t, pts, p_min[0], p_min[1], p_min[2]) == ON_POSITIVE_SIDE); + } + + CGAL::Orientation side = orient(t, pts, bbox.xmin(), bbox.ymin(),bbox.zmin()); + if(side == COPLANAR) return true; + CGAL::Oriented_side s = orient(t, pts, bbox.xmax(), bbox.ymax(),bbox.zmax()); + if(s != side) return true; + s = orient(t, pts, bbox.xmin(), bbox.ymin(),bbox.zmax()); + if(s != side) return true; + s = orient(t, pts, bbox.xmax(), bbox.ymax(),bbox.zmin()); + if(s != side) return true; + s = orient(t, pts, bbox.xmin(), bbox.ymax(),bbox.zmin()); + if(s != side) return true; + s = orient(t, pts, bbox.xmax(), bbox.ymin(),bbox.zmax()); + if(s != side) return true; + s = orient(t, pts, bbox.xmin(), bbox.ymax(),bbox.zmax()); + if(s != side) return true; + s = orient(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 @@ -274,6 +520,7 @@ class Do_intersect_3 // 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 { + // TODO: shall we add the case 0 if one of the alpha's and beta's are 0? Sign int_tmp_result; double double_tmp_result; double eps; @@ -333,29 +580,7 @@ class Do_intersect_3 return int_tmp_result; }; - auto do_intersect_plane_bbox = [](const Triangle_3& t, const Bbox_3& bbox) - { - typename SFK::Orientation_3 orient = SFK().orientation_3_object(); - CGAL::Orientation side = orient(t[0], t[1],t[2], Point_3(bbox.xmin(), bbox.ymin(),bbox.zmin())); - if(side == COPLANAR) return true; - CGAL::Oriented_side s = orient(t[0], t[1],t[2], Point_3(bbox.xmax(), bbox.ymax(),bbox.zmax())); - if(s != side) return true; - s = orient(t[0], t[1],t[2], Point_3(bbox.xmin(), bbox.ymin(),bbox.zmax())); - if(s != side) return true; - s = orient(t[0], t[1],t[2], Point_3(bbox.xmax(), bbox.ymax(),bbox.zmin())); - if(s != side) return true; - s = orient(t[0], t[1],t[2], Point_3(bbox.xmin(), bbox.ymax(),bbox.zmin())); - if(s != side) return true; - s = orient(t[0], t[1],t[2], Point_3(bbox.xmax(), bbox.ymin(),bbox.zmax())); - if(s != side) return true; - s = orient(t[0], t[1],t[2], Point_3(bbox.xmin(), bbox.ymax(),bbox.zmax())); - if(s != side) return true; - s = orient(t[0], t[1],t[2], Point_3(bbox.xmax(), bbox.ymin(),bbox.zmin())); - if(s != side) return true; - return false; - }; - - if ( !do_intersect_plane_bbox(t,b) ) + 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); From 77b77adbc44a6e8dbdea60347983f272f769f04e Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 17 Jun 2021 09:03:33 +0200 Subject: [PATCH 13/30] fabs -> CGAL::abs Does the compiler optimize two identical calls of CGAL::abs away? --- .../internal/Static_filters/Do_intersect_3.h | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 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 5850056ef433..1c8a571a0b91 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 @@ -253,15 +253,15 @@ class Do_intersect_3 double double_tmp_result; double eps; double_tmp_result = ((qx_px * ry_py) - (rx_px * qy_py)); - double max1 = fabs(qx_px); - if( (max1 < fabs(rx_px)) ) + double max1 = CGAL::abs(qx_px); + if( (max1 < CGAL::abs(rx_px)) ) { - max1 = fabs(rx_px); + max1 = CGAL::abs(rx_px); } - double max2 = fabs(ry_py); - if( (max2 < fabs(qy_py)) ) + double max2 = CGAL::abs(ry_py); + if( (max2 < CGAL::abs(qy_py)) ) { - max2 = fabs(qy_py); + max2 = CGAL::abs(qy_py); } double lower_bound_1; double upper_bound_1; @@ -525,15 +525,15 @@ class Do_intersect_3 double double_tmp_result; double eps; double_tmp_result = ((-c_alpha * alpha) + (c_beta * beta)); - double max1 = fabs(c_alpha); - if( (max1 < fabs(c_beta)) ) + double max1 = CGAL::abs(c_alpha); + if( (max1 < CGAL::abs(c_beta)) ) { - max1 = fabs(c_beta); + max1 = CGAL::abs(c_beta); } - double max2 = fabs(alpha); - if( (max2 < fabs(beta)) ) + double max2 = CGAL::abs(alpha); + if( (max2 < CGAL::abs(beta)) ) { - max2 = fabs(beta); + max2 = CGAL::abs(beta); } double lower_bound_1; double upper_bound_1; From 46e2042df46fca91373b152980d243a99eb71d91 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 17 Jun 2021 12:07:47 +0200 Subject: [PATCH 14/30] Add profiler --- .../internal/Static_filters/Do_intersect_3.h | 23 ++++++++++++++++--- 1 file changed, 20 insertions(+), 3 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 1c8a571a0b91..f3e9b1a2fa86 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 @@ -245,6 +245,8 @@ class Do_intersect_3 /// TODO: move in utils 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); @@ -280,12 +282,14 @@ class Do_intersect_3 } if( (lower_bound_1 < 5.00368081960964746551e-147) ) { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); return Uncertain(); } else { if( (upper_bound_1 > 1.67597599124282389316e+153) ) { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); return Uncertain(); } eps = (8.88720573725927976811e-16 * (max1 * max2)); @@ -301,24 +305,37 @@ class Do_intersect_3 } else { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); return Uncertain(); } } } + 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)) return false; + if (is_indeterminate(s)){ + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + return false; + } signs[0]=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)) return false; + if (is_indeterminate(s)){ + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + return false; + } signs[1]=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)) return false; + if (is_indeterminate(s)){ + CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + return false; + } signs[2]=s; + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return true; } From 8844a94d42d5a46901dc8c93a3f074b8f2c1c760 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 17 Jun 2021 12:44:20 +0200 Subject: [PATCH 15/30] Use Uncertain<>::indeterminante() --- .../CGAL/internal/Static_filters/Do_intersect_3.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 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 f3e9b1a2fa86..0c179d2aba33 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 @@ -283,14 +283,14 @@ class Do_intersect_3 if( (lower_bound_1 < 5.00368081960964746551e-147) ) { CGAL_BRANCH_PROFILER_BRANCH_1(tmp); - return Uncertain(); + return Uncertain::indeterminate(); } else { if( (upper_bound_1 > 1.67597599124282389316e+153) ) { CGAL_BRANCH_PROFILER_BRANCH_1(tmp); - return Uncertain(); + return Uncertain::indeterminate(); } eps = (8.88720573725927976811e-16 * (max1 * max2)); if( (double_tmp_result > eps) ) @@ -306,7 +306,7 @@ class Do_intersect_3 else { CGAL_BRANCH_PROFILER_BRANCH_1(tmp); - return Uncertain(); + return Uncertain::indeterminate(); } } } @@ -418,7 +418,7 @@ class Do_intersect_3 if (det > eps) return POSITIVE; if (det < -eps) return NEGATIVE; } - return Uncertain(); + return Uncertain::indeterminate(); }; auto orient = [my_orient](const Triangle_3& t, const std::array< std::array, 3>& pts, @@ -569,13 +569,13 @@ class Do_intersect_3 } if( (lower_bound_1 < 5.00368081960964746551e-147) ) { - return Uncertain(); + return Uncertain::indeterminate(); } else { if( (upper_bound_1 > 1.67597599124282389316e+153) ) { - return Uncertain(); + return Uncertain::indeterminate(); } eps = (8.88720573725927976811e-16 * (max1 * max2)); if( (double_tmp_result > eps) ) @@ -590,7 +590,7 @@ class Do_intersect_3 } else { - return Uncertain(); + return Uncertain::indeterminate(); } } } From 4c16735f85ac088658cca7a3aac932201562e168 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 17 Jun 2021 12:55:53 +0200 Subject: [PATCH 16/30] make_certain() --- .../include/CGAL/internal/Static_filters/Do_intersect_3.h | 6 +++--- 1 file changed, 3 insertions(+), 3 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 0c179d2aba33..74f8dee1828b 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 @@ -322,19 +322,19 @@ class Do_intersect_3 CGAL_BRANCH_PROFILER_BRANCH_1(tmp); return false; } - signs[0]=s; + 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]=s; + 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]=s; + signs[2]=make_certain(s); CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return true; } From 9bfe21144e8fbee1da42b1c66bb222965c57bb92 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 17 Jun 2021 12:56:32 +0200 Subject: [PATCH 17/30] fix incorrect return type --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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 6af9fe95a17c..a8bb60ff3a2f 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 @@ -140,7 +140,7 @@ namespace internal { template inline - Sign + Uncertain do_axis_intersect_aux_A0(const FT& alpha, const FT& beta, const std::array, 3>& sides, @@ -151,7 +151,7 @@ namespace internal { template inline - Sign + Uncertain do_axis_intersect_aux_A1(const FT& alpha, const FT& beta, const std::array, 3>& sides, @@ -162,7 +162,7 @@ namespace internal { template inline - Sign + Uncertain do_axis_intersect_aux_A2(const FT& alpha, const FT& beta, const std::array, 3>& sides, From 60559d799ae4f3291e18b15f0d1feee7c493a11b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 17 Jun 2021 13:08:34 +0200 Subject: [PATCH 18/30] do no call twice the failing static filter --- .../include/CGAL/internal/Static_filters/Do_intersect_3.h | 3 +-- 1 file changed, 1 insertion(+), 2 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 74f8dee1828b..e851bf9cfe03 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 @@ -426,8 +426,7 @@ class Do_intersect_3 { Uncertain res = my_orient(pts, x, y, z); if (!is_indeterminate(res)) return make_certain(res); - typename SFK::Orientation_3 orient = SFK().orientation_3_object(); - //TODO: avoid calling the static filter again... + typename K_base::Orientation_3 orient = K_base().orientation_3_object(); return orient(t[0], t[1], t[2], Point_3(x,y,z)); }; //// From d13c88710df84e543f933adddb8056d56b471989 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 17 Jun 2021 16:02:10 +0200 Subject: [PATCH 19/30] clean up --- .../internal/Static_filters/Do_intersect_3.h | 43 ++++++++++--------- 1 file changed, 23 insertions(+), 20 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 e851bf9cfe03..cdfeb12fb9e2 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 @@ -339,9 +339,13 @@ class Do_intersect_3 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 { - auto my_orient = [](const std::array< std::array, 3>& pts, double x, double y, double z) -> Uncertain + // 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 { double pqx = pts[1][0] - pts[0][0]; double pqy = pts[1][1] - pts[0][1]; @@ -421,21 +425,23 @@ class Do_intersect_3 return Uncertain::indeterminate(); }; - auto orient = [my_orient](const Triangle_3& t, const std::array< std::array, 3>& pts, - double x, double y, double z) -> Orientation + 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 = my_orient(pts, x, y, z); + 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(); + 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 + // 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) { @@ -480,31 +486,29 @@ class Do_intersect_3 } } - return ! (orient(t, pts, p_max[0], p_max[1], p_max[2]) == ON_NEGATIVE_SIDE || - orient(t, pts, p_min[0], p_min[1], p_min[2]) == ON_POSITIVE_SIDE); + 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 = orient(t, pts, bbox.xmin(), bbox.ymin(),bbox.zmin()); + CGAL::Orientation side = orientation(t, pts, bbox.xmin(), bbox.ymin(),bbox.zmin()); if(side == COPLANAR) return true; - CGAL::Oriented_side s = orient(t, pts, bbox.xmax(), bbox.ymax(),bbox.zmax()); + CGAL::Oriented_side s = orientation(t, pts, bbox.xmax(), bbox.ymax(),bbox.zmax()); if(s != side) return true; - s = orient(t, pts, bbox.xmin(), bbox.ymin(),bbox.zmax()); + s = orientation(t, pts, bbox.xmin(), bbox.ymin(),bbox.zmax()); if(s != side) return true; - s = orient(t, pts, bbox.xmax(), bbox.ymax(),bbox.zmin()); + s = orientation(t, pts, bbox.xmax(), bbox.ymax(),bbox.zmin()); if(s != side) return true; - s = orient(t, pts, bbox.xmin(), bbox.ymax(),bbox.zmin()); + s = orientation(t, pts, bbox.xmin(), bbox.ymax(),bbox.zmin()); if(s != side) return true; - s = orient(t, pts, bbox.xmax(), bbox.ymin(),bbox.zmax()); + s = orientation(t, pts, bbox.xmax(), bbox.ymin(),bbox.zmax()); if(s != side) return true; - s = orient(t, pts, bbox.xmin(), bbox.ymax(),bbox.zmax()); + s = orientation(t, pts, bbox.xmin(), bbox.ymax(),bbox.zmax()); if(s != side) return true; - s = orient(t, pts, bbox.xmax(), bbox.ymin(),bbox.zmin()); + 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 { @@ -536,7 +540,6 @@ class Do_intersect_3 // 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 { - // TODO: shall we add the case 0 if one of the alpha's and beta's are 0? Sign int_tmp_result; double double_tmp_result; double eps; From 3d04ab1e9ebfa6c84d13724a8efc6bbe1d8ed5db Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 13 Aug 2021 10:27:36 +0100 Subject: [PATCH 20/30] Explicitely use CGAL::sign This only fixes some compilation errors, but not the runtime errors of AABB_tree for VC++ --- .../Intersections_3/internal/Bbox_3_Triangle_3_do_intersect.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 a8bb60ff3a2f..9cff1b4b690d 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 @@ -381,7 +381,7 @@ namespace internal { const FT& c_alpha, const FT& c_beta) -> Uncertain { - return sign(- c_alpha * alpha + c_beta * beta); + return CGAL::sign(- c_alpha * alpha + c_beta * beta); }; std::array< std::array, 3> triangle = From 7d2321cecd791812de3142e70da2165959bf9b35 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 25 Aug 2021 16:05:03 +0100 Subject: [PATCH 21/30] Work around CGAL_AND as explained in the manual --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) 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 9cff1b4b690d..e6b2b9d097f8 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 @@ -213,8 +213,13 @@ namespace internal { if (is_indeterminate(b)) return b; if(b) std::swap(j,k); - return CGAL_AND((do_axis_intersect_aux_A0(p_min[1]-j[1], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), - (do_axis_intersect_aux_A0(p_max[1]-k[1], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); + Uncertain tmp = do_axis_intersect_aux_A0(p_min[1]-j[1], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE; + + return certainly_not(tmp) ? make_uncertain(false) : tmp & (do_axis_intersect_aux_A0(p_max[1]-k[1], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE); + + + //return CGAL_AND((do_axis_intersect_aux_A0(p_min[1]-j[1], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), + // (do_axis_intersect_aux_A0(p_max[1]-k[1], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } case 1: { // t_max >= t_min @@ -222,8 +227,9 @@ namespace internal { if (is_indeterminate(b)) return b; if(b) std::swap(j,k); - return CGAL_AND((do_axis_intersect_aux_A1(p_min[0]-j[0], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), - (do_axis_intersect_aux_A1(p_max[0]-k[0], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); + + Uncertain tmp = do_axis_intersect_aux_A1(p_min[0]-j[0], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE; + return certainly_not(tmp) ? make_uncertain(false) : tmp & (do_axis_intersect_aux_A1(p_max[0]-k[0], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE); } case 2: { @@ -232,8 +238,9 @@ namespace internal { if ( is_indeterminate(b)) return b; if(b) std::swap(j,k); - return CGAL_AND((do_axis_intersect_aux_A2(p_min[0]-j[0], p_min[1]-j[1], sides, do_axis_intersect_aux_impl) != POSITIVE), - (do_axis_intersect_aux_A2(p_max[0]-k[0], p_max[1]-k[1], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); + + Uncertain tmp = do_axis_intersect_aux_A2(p_min[0]-j[0], p_min[1]-j[1], sides, do_axis_intersect_aux_impl) != POSITIVE; + return certainly_not(tmp) ? make_uncertain(false) : tmp &(do_axis_intersect_aux_A2(p_max[0]-k[0], p_max[1]-k[1], sides, do_axis_intersect_aux_impl) != NEGATIVE); } default: // Should not happen From 10a74957efb3f0ec0627e2e2b10538fdb865178f Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 25 Aug 2021 16:15:55 +0100 Subject: [PATCH 22/30] Remove the test against the bbox of the triangle/tetrahedron. This should be done before, for example in the AABB tree --- .../include/CGAL/internal/Static_filters/Do_intersect_3.h | 8 -------- 1 file changed, 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 cdfeb12fb9e2..3a5e2c1bd8ca 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 @@ -135,9 +135,6 @@ class Do_intersect_3 Get_approx get_approx; double px, py, pz; - if(! do_overlap(t.bbox(),b)){ - return false; - } { const Point_3& p = t[0]; if (fit_in_double(get_approx(p).x(), px) && fit_in_double(get_approx(p).y(), py) && @@ -512,11 +509,6 @@ class Do_intersect_3 result_type operator()(const Triangle_3 &t, const Bbox_3& b) const { - // check overlaps between bboxes - if(! do_overlap(t.bbox(),b)){ - return false; - } - // check if at least one triangle point is inside the bbox Get_approx get_approx; std::array< std::array, 3> pts; From a95858fc555f8ee22fd7f3536da9f5af62cd3390 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 31 Aug 2021 08:56:56 +0100 Subject: [PATCH 23/30] Revert "Work around CGAL_AND as explained in the manual" This reverts commit 7d2321cecd791812de3142e70da2165959bf9b35. --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) 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 e6b2b9d097f8..9cff1b4b690d 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 @@ -213,13 +213,8 @@ namespace internal { if (is_indeterminate(b)) return b; if(b) std::swap(j,k); - Uncertain tmp = do_axis_intersect_aux_A0(p_min[1]-j[1], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE; - - return certainly_not(tmp) ? make_uncertain(false) : tmp & (do_axis_intersect_aux_A0(p_max[1]-k[1], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE); - - - //return CGAL_AND((do_axis_intersect_aux_A0(p_min[1]-j[1], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), - // (do_axis_intersect_aux_A0(p_max[1]-k[1], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); + return CGAL_AND((do_axis_intersect_aux_A0(p_min[1]-j[1], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), + (do_axis_intersect_aux_A0(p_max[1]-k[1], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } case 1: { // t_max >= t_min @@ -227,9 +222,8 @@ namespace internal { if (is_indeterminate(b)) return b; if(b) std::swap(j,k); - - Uncertain tmp = do_axis_intersect_aux_A1(p_min[0]-j[0], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE; - return certainly_not(tmp) ? make_uncertain(false) : tmp & (do_axis_intersect_aux_A1(p_max[0]-k[0], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE); + return CGAL_AND((do_axis_intersect_aux_A1(p_min[0]-j[0], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), + (do_axis_intersect_aux_A1(p_max[0]-k[0], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } case 2: { @@ -238,9 +232,8 @@ namespace internal { if ( is_indeterminate(b)) return b; if(b) std::swap(j,k); - - Uncertain tmp = do_axis_intersect_aux_A2(p_min[0]-j[0], p_min[1]-j[1], sides, do_axis_intersect_aux_impl) != POSITIVE; - return certainly_not(tmp) ? make_uncertain(false) : tmp &(do_axis_intersect_aux_A2(p_max[0]-k[0], p_max[1]-k[1], sides, do_axis_intersect_aux_impl) != NEGATIVE); + return CGAL_AND((do_axis_intersect_aux_A2(p_min[0]-j[0], p_min[1]-j[1], sides, do_axis_intersect_aux_impl) != POSITIVE), + (do_axis_intersect_aux_A2(p_max[0]-k[0], p_max[1]-k[1], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); } default: // Should not happen From 7049d667e5b1574c0f91d0c827574fff3f3988f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 10 Sep 2021 18:12:30 +0200 Subject: [PATCH 24/30] Avoid copying arrays --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 92 ++++++++++--------- 1 file changed, 48 insertions(+), 44 deletions(-) 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 5b87ed6678af..ccc948e18735 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 @@ -191,14 +191,15 @@ collinear_axis(const std::array& side) template inline -Uncertain do_axis_intersect(const std::array< std::array , 3>& triangle, - const std::array< std::array , 3>& sides, +Uncertain do_axis_intersect(const std::array, 3>& triangle, + const std::array, 3>& sides, const Box3& bbox, Fct do_axis_intersect_aux_impl) { - // TODO: no longer a ptr, is it an issue for swap? - std::array j = triangle[SIDE]; - std::array k = triangle[(SIDE+2)%3]; + 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], @@ -208,48 +209,51 @@ Uncertain do_axis_intersect(const std::array< std::array , 3>& tria switch(AXE) { - case 0: + case 0: { - // t_max >= t_min - 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_A0(p_min[1]-j[1], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), - (do_axis_intersect_aux_A0(p_max[1]-k[1], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); - } - case 1: + // t_max >= t_min + 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_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: { - // t_max >= t_min - 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_A1(p_min[0]-j[0], p_min[2]-j[2], sides, do_axis_intersect_aux_impl) != POSITIVE), - (do_axis_intersect_aux_A1(p_max[0]-k[0], p_max[2]-k[2], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); - - } - case 2: + // t_max >= t_min + 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_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: { - // t_max >= t_min - 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_A2(p_min[0]-j[0], p_min[1]-j[1], sides, do_axis_intersect_aux_impl) != POSITIVE), - (do_axis_intersect_aux_A2(p_max[0]-k[0], p_max[1]-k[1], sides, do_axis_intersect_aux_impl) != NEGATIVE) ); - } - default: - // Should not happen - CGAL_error(); - return make_uncertain(false); + // t_max >= t_min + 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_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 + CGAL_error(); + return make_uncertain(false); } } From b67589e16b9db85d0a66bbeb6e2c4024e67a03bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 14 Sep 2021 08:41:34 +0200 Subject: [PATCH 25/30] workaround sign issue --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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 ccc948e18735..31a861df29ce 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 @@ -24,6 +24,8 @@ #include #include +#include + namespace CGAL { namespace Intersections { namespace internal { @@ -411,7 +413,8 @@ bool do_intersect_bbox_or_iso_cuboid(const typename K::Triangle_3& a_triangle, const FT& c_alpha, const FT& c_beta) -> Uncertain { - return CGAL::sign(- c_alpha * alpha + c_beta * beta); + Sgn sign; + return sign(- c_alpha * alpha + c_beta * beta); }; std::array< std::array, 3> triangle = From 1e9cad98a864bf33a4919413d4709bb3e8c1b2c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 14 Sep 2021 09:32:40 +0200 Subject: [PATCH 26/30] better fix --- .../internal/Bbox_3_Triangle_3_do_intersect.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) 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 31a861df29ce..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 @@ -24,7 +24,7 @@ #include #include -#include +#include namespace CGAL { namespace Intersections { @@ -413,8 +413,7 @@ bool do_intersect_bbox_or_iso_cuboid(const typename K::Triangle_3& a_triangle, const FT& c_alpha, const FT& c_beta) -> Uncertain { - Sgn sign; - return sign(- c_alpha * alpha + c_beta * beta); + return CGAL::sign(- c_alpha * alpha + c_beta * beta); }; std::array< std::array, 3> triangle = From c8d2761f45f9fd1b24cec4276e3f64fa94a31348 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 29 Sep 2021 15:14:54 +0100 Subject: [PATCH 27/30] partial work after Mael's review --- .../Filtered_kernel/internal/Do_intersect_3.h | 50 ++++--------------- 1 file changed, 9 insertions(+), 41 deletions(-) diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h index 36a87dd83965..cdaec52c3318 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h @@ -128,15 +128,14 @@ class Do_intersect_3 } - result_type operator()(const Tetrahedron_3 &t, const Bbox_3& b) const { Get_approx get_approx; double px, py, pz; - { - const Point_3& p = t[0]; + 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) ) { @@ -148,49 +147,18 @@ class Do_intersect_3 } } - { - const Point_3& p = t[1]; - 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) ) - { - if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - return true; - } - }else{ - return Base::operator()(t,b); - } - } - - { - const Point_3& p = t[2]; - 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) ) - { - if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - return true; - } - }else{ - return Base::operator()(t,b); - } - } - - + return Base::operator()(t,b); + } - { - const Point_3& p = t[3]; - 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) ) - { - if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - return true; - } - } - } - return Base::operator()(t,b); + result_type + operator()(const Bbox_3& b, const Tetrahedron_3 &t) const + { + return this->operator()(t, b); } + result_type operator()(const Bbox_3& b, const Ray_3 &r) const { From 37f25419426aae2aae4b49b9c9daa5449ac94d41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 8 Oct 2021 20:07:50 +0200 Subject: [PATCH 28/30] Add profiler macros --- .../Filtered_kernel/internal/Do_intersect_3.h | 37 ++++++++++++------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h index cdaec52c3318..499ea481b738 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h @@ -127,22 +127,40 @@ 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) { + 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()) ) { - if( (px >= b.xmin() && px <= b.xmax()) && (py >= b.ymin() && py <= b.ymax()) && (pz >= b.zmin() && pz <= b.zmax()) ){ - return true; - } - }else{ + return true; + } + + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + } + else + { return Base::operator()(t,b); } } @@ -150,15 +168,6 @@ class Do_intersect_3 return Base::operator()(t,b); } - - result_type - operator()(const Bbox_3& b, const Tetrahedron_3 &t) const - { - return this->operator()(t, b); - } - - - result_type operator()(const Bbox_3& b, const Ray_3 &r) const { From 180f9c8e8d5d549f174f24703efd1bc9a1d7a645 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 8 Oct 2021 20:08:01 +0200 Subject: [PATCH 29/30] Remove obsolete todo --- .../include/CGAL/Filtered_kernel/internal/Do_intersect_3.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h index 499ea481b738..b3baa9630296 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Do_intersect_3.h @@ -216,7 +216,6 @@ class Do_intersect_3 return this->operator()(t, b); } - /// TODO: move in utils 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); From 381ff61796eff0b2a974ef69af131a848d1405ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 8 Oct 2021 21:52:32 +0200 Subject: [PATCH 30/30] More profiler macros --- .../internal/Static_filters/Do_intersect_3.h | 24 ++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) 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 15c66af9ff01..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 @@ -320,6 +320,9 @@ class Do_intersect_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]; @@ -395,6 +398,8 @@ class Do_intersect_3 if (det > eps) return POSITIVE; if (det < -eps) return NEGATIVE; } + + CGAL_BRANCH_PROFILER_BRANCH(tmp); return Uncertain::indeterminate(); }; @@ -485,6 +490,9 @@ class Do_intersect_3 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; @@ -495,9 +503,16 @@ class Do_intersect_3 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]) ) { - 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()) ){ + 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 { @@ -508,6 +523,8 @@ class Do_intersect_3 // 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; @@ -539,12 +556,14 @@ class Do_intersect_3 } 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)); @@ -560,10 +579,13 @@ class Do_intersect_3 } else { + CGAL_BRANCH_PROFILER_BRANCH_1(tmp2); return Uncertain::indeterminate(); } } } + + CGAL_BRANCH_PROFILER_BRANCH_2(tmp2); return int_tmp_result; };