From bf89c548fdc45608e0ae7ef00936a11e6cd4fd50 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 1 Aug 2022 11:43:27 +0200 Subject: [PATCH 01/10] simplify Tetra/triangle intersection computation Testsuite fails for Kernel with non-exact contructions --- .../Tetrahedron_3_Triangle_3_intersection.h | 586 ++---------------- 1 file changed, 45 insertions(+), 541 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h index 272a07da6c0..263472803ef 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h @@ -15,8 +15,7 @@ #ifndef CGAL_INTERNAL_INTERSECTIONS_3_TETRAHEDRON_3_TRIANGLE_3_INTERSECTIONS_H #define CGAL_INTERNAL_INTERSECTIONS_3_TETRAHEDRON_3_TRIANGLE_3_INTERSECTIONS_H -#include -#include +#include #include @@ -30,255 +29,6 @@ namespace CGAL { namespace Intersections { namespace internal { -template -void filter_segments(std::list& segments) -{ - auto are_equal = [](const Segment& l, const Segment& r) -> bool - { - return (l == r || l == r.opposite()); - }; - - auto it = std::unique(segments.begin(), segments.end(), are_equal); - segments.erase(it, segments.end()); -} - -// plane going through the ref segment's source, a point above (given by the normal of the input -// triangle) and two points (ref_other / query) that need to be ordered -template -bool first_comes_first_pt(const typename K::Point_3& ref_source, - const typename K::Point_3& ref_z, - const typename K::Point_3& ref_other, - const typename K::Point_3& query, - const K& k) -{ - typename K::Orientation_3 orientation = k.orientation_3_object(); - - // points have filtered to remove segments' extremities - CGAL_precondition(ref_other != query); - - const Orientation o = orientation(ref_source, ref_z, ref_other, query); - CGAL_assertion(o != COPLANAR); - - // ref_other comes first <==> query is on the positive side of the plane - return (o == POSITIVE); -} - -template -bool first_comes_first(const typename K::Point_3& ref_source, - const typename K::Point_3& ref_z, - const typename K::Point_3& ref_other, - const SegPtVariant& seg_or_pt, - const K& k) -{ - typedef typename K::Point_3 Point_3; - typedef typename K::Segment_3 Segment_3; - - typedef typename std::list::iterator SCI; - typedef typename std::vector::iterator PCI; - - if(seg_or_pt.which() == 0) - { - const Segment_3& s = *(boost::get(seg_or_pt)); - return first_comes_first_pt(ref_source, ref_z, ref_other, s.source(), k); - } - else - { - CGAL_assertion(seg_or_pt.which() == 1); - - const Point_3& p = *(boost::get(seg_or_pt)); - return first_comes_first_pt(ref_source, ref_z, ref_other, p, k); - } -} - -template -typename Intersection_traits::result_type -build_intersection(const typename K::Tetrahedron_3& /*input_tetrahedron*/, - const typename K::Triangle_3& input_triangle, - PointContainer& points, - SegmentContainer& segments, - const K& k) -{ - typedef typename Intersection_traits::result_type result_type; - - typedef typename K::Point_3 Point_3; - typedef typename K::Segment_3 Segment_3; - typedef typename K::Vector_3 Vector_3; - typedef typename K::Triangle_3 Triangle_3; - typedef std::vector Poly; - - typedef typename SegmentContainer::iterator SCI; - typedef typename PointContainer::iterator PCI; - - // @todo? Could do the 1 segment case with this code too... - CGAL_precondition(segments.size() >= 2 && segments.size() <= 4); - CGAL_precondition(points.size() <= 2); - - // Constructions @fixme avoidable? - const Vector_3 input_triangle_normal = input_triangle.supporting_plane().orthogonal_vector(); - - // remove points that are just segments extremities - auto is_extremity = [&segments](const Point_3& p) -> bool - { - for(const Segment_3& s : segments) - if(p == s.source() || p == s.target()) - return true; - return false; - }; - points.erase(std::remove_if(points.begin(), points.end(), is_extremity), - points.end()); - - // Take the first segment as reference, and order the rest to form a convex polygon - // - // All segments and points involved in the intersection are on the input triangle - // and thus everything is coplanar, at least theoretically (the kernel might not provide - // exact constructions...) - // - // Given an arbitrary segment, the code below sorts the other segments and the points - // in a ccw order. Using a vector because the number of segments and points is bounded - // (max 4 segments and max 2 points) so even if the linear insertion is a little ugly, - // it is not expensive anyway. - // - // Example: - /* - x p0 - - / - / - s1 / \ s2 - / \ - -------------- - s0 - */ - // - // s0 is chosen as the reference segment - // output will be 's0 s2 p0 s1' - - Segment_3& ref_s = segments.front(); - Point_3 ref_z = ref_s.source() + input_triangle_normal; - - // The reference segment should be such that all other intersection parts are - // on the positive side of the plane described by the normal of the triangle and ref_s - bool swapped = false; - for(SCI slit = std::next(segments.begin()); slit != segments.end(); ++slit) - { - const Segment_3& other = *slit; - - if(k.orientation_3_object()(ref_s.source(), ref_z, ref_s.target(), other.source()) == CGAL::NEGATIVE || - k.orientation_3_object()(ref_s.source(), ref_z, ref_s.target(), other.target()) == CGAL::NEGATIVE) - { - ref_s = ref_s.opposite(); - ref_z = ref_s.source() + input_triangle_normal; - swapped = true; - break; - } - } - - if(!swapped) - { - for(PCI plit = points.begin(); plit != points.end(); ++plit) - { - const Point_3& other = *plit; - if(k.orientation_3_object()(ref_s.source(), ref_z, ref_s.target(), other) == CGAL::NEGATIVE) - { - swapped = true; - ref_s = ref_s.opposite(); - ref_z = ref_s.source() + input_triangle_normal; - break; - } - } - } - - const Point_3& ref_sp = ref_s.source(); - const Point_3& ref_tp = ref_s.target(); - - // Now, order the other parts of the intersection - std::list > res_elements; // iterators to the points/segments - res_elements.emplace_back(segments.begin()); - - for(SCI slit = std::next(segments.begin()); slit != segments.end(); ++slit) - { - // first, check if the segment is well oriented, meaning its source comes before its target (ccwly) - Segment_3& curr_s = *slit; - - if(curr_s.source() == ref_sp || curr_s.target() == ref_tp) // consecutive segments must have consistent orientation - { - curr_s = curr_s.opposite(); - } - else if(curr_s.source() == ref_tp || curr_s.target() == ref_sp) - { - // nothing to do here as we know that sp&tp are on the positive side of (normal, ref_s) - } - else if(first_comes_first_pt(ref_sp, ref_z, curr_s.target(), curr_s.source(), k)) - { - curr_s = curr_s.opposite(); - } - - // Find where the current segment fit in the final polygon intersection - for(auto rit = std::next(res_elements.begin()); ; ++rit) - { - // always pick the current segment's source to ensure ref_source != ref_other - if(rit == res_elements.end() || first_comes_first(ref_sp, ref_z, curr_s.source(), *rit, k)) - { - res_elements.insert(rit, slit); - break; - } - } - } - - for(PCI plit = points.begin(); plit != points.end(); ++plit) - { - const Point_3& curr_p = *plit; - - // Find where the current point fits in the boundary of the polygon intersection - for(auto rit = std::next(res_elements.begin()); ; ++rit) - { - if(rit == res_elements.end() || first_comes_first(ref_sp, ref_z, curr_p, *rit, k)) - { - res_elements.insert(rit, plit); - break; - } - } - } - - CGAL_postcondition(res_elements.size() == points.size() + segments.size()); - - // Concatenate points to create the polygonal output - Poly res; - for(const boost::variant& e : res_elements) - { - if(const SCI* sci = boost::get(&e)) - { - const Segment_3& s = **sci; - - if(res.empty() || s.source() != res.back()) // common extremity for consecutive segments - res.push_back(s.source()); - if(res.empty() || s.target() != res.front()) - res.push_back(s.target()); - } - else if(const PCI* pci = boost::get(&e)) - { - res.push_back(**pci); - } - else - { - CGAL_assertion(false); - } - } - - CGAL_assertion(std::set(res.begin(), res.end()).size() == res.size()); - CGAL_assertion(res.size() >= 3); - - if(res.size() == 3) - { - Triangle_3 tr { res[0], res[1], res[2] }; - return result_type(std::forward(tr)); - } - else - { - return result_type(std::forward(res)); - } -} - template typename Intersection_traits::result_type intersection(const typename K::Tetrahedron_3& tet, @@ -286,319 +36,73 @@ intersection(const typename K::Tetrahedron_3& tet, const K& k) { typedef typename Intersection_traits::result_type result_type; - typedef typename Intersection_traits::result_type Inter_type; CGAL_precondition(!tet.is_degenerate()); CGAL_precondition(!tr.is_degenerate()); typedef typename K::Point_3 Point_3; typedef typename K::Segment_3 Segment_3; + typedef typename K::Plane_3 Plane_3; typedef typename K::Triangle_3 Triangle_3; - typedef std::vector Poly; + typedef typename K::Line_3 Line_3; + typedef typename K::Point_3 Point_3; - typename K::Bounded_side_3 bounded_side = k.bounded_side_3_object(); + typename K::Construct_plane_3 plane = k.construct_plane_3_object(); typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object(); typename K::Construct_triangle_3 triangle = k.construct_triangle_3_object(); + typename K::Construct_segment_3 segment = k.construct_segment_3_object(); + typename K::Construct_line_3 line = k.construct_line_3_object(); + typename K::Oriented_side_3 oriented_side = k.oriented_side_3_object(); - std::vector vertex_sides(3); - std::vector points; - int inside_points = 0; - int strictly_inside_points = 0; + std::vector res = { vertex(tr,0), vertex(tr,1), vertex(tr,2) }; - for(int i=0; i<3; ++i) + for (int i=0; i<4; ++i) { - vertex_sides[i] = bounded_side(tet, vertex(tr, i)); + Plane_3 pl = plane(vertex(tet, (1+i)%4), vertex(tet, (2+i)%4),vertex(tet, (3+i)%4)); + if (oriented_side(pl, vertex(tet,i))!=ON_POSITIVE_SIDE) // TODO: this should be precomputed based on the tetra orientation + pl = pl.opposite(); - if(vertex_sides[i] != ON_UNBOUNDED_SIDE) - ++inside_points; + std::vector current; + std::vector orientations(res.size()); + for (std::size_t i=0; i segments; - std::vector seg_ids; - for(std::size_t i = 0; i < 4; ++i) + const bool test_segment = i!=1 || res.size()!=2; + std::size_t j = (i+1)%res.size(); + switch(orientations[j]) { - const Triangle_3 face = triangle(vertex(tet, (i+1)%4), - vertex(tet, (i+2)%4), - vertex(tet, (i+3)%4)); - intersections[i] = intersection(tr, face, k); - if(intersections[i]) - { - // a face is inside the input tr - if(const Triangle_3* t = boost::get(&*intersections[i])) - { - Triangle_3 res = *t; - return result_type(std::forward(res)); - } - else if(const Segment_3* s = boost::get(&*intersections[i])) - { - // get segs and pts to construct poly - segments.push_back(*s); - seg_ids.push_back(i); - } - else if(const Point_3* p = boost::get(&*intersections[i])) - { - points.push_back(*p); - } - else if(const Poly* p = boost::get(&*intersections[i])) - { - // the input triangle is in the supporting plane of a tet face, return the poly. - Poly res = *p; - return result_type(std::forward(res)); - } - } - } - - if(segments.size() > 1) - filter_segments(segments); - - // no segments and no inside points, there can still be an intersecting (tet vertex on - // an edge|face of the triangle) - if(segments.empty()) - { - if(points.empty()) - return result_type(); - - return result_type(std::forward(points.front())); - } - else if(segments.size() == 1) - { - // adjacency to an edge, return resulting segment. - return result_type(segments.front()); - } - else if(segments.size() > 1) - { - return build_intersection(tet, tr, points, segments, k); - } - } - break; - case 1: - case 2: // 1 or 2 inside points - { - Inter_type intersections[4]; - std::list segments; - for(std::size_t i = 0; i < 4; ++i) - { - const Triangle_3 face = triangle(vertex(tet, (i+1)%4), - vertex(tet, (i+2)%4), - vertex(tet, (i+3)%4)); - intersections[i] = intersection(tr, face, k); - if(intersections[i]) - { - if(const Triangle_3* t = boost::get(&*intersections[i])) - { - Triangle_3 res = *t; - return result_type(std::forward(res)); - } - else if(const Segment_3* s = boost::get(&*intersections[i])) - { - segments.push_back(*s); - } - else if(const Point_3* p = boost::get(&*intersections[i])) - { - points.push_back(*p); - } - else if(const Poly* p = boost::get(&*intersections[i])) - { - // the input is in a supporting plane of a face - Poly res = *p; - return result_type(std::forward(res)); - } - } - } - - if(segments.size() > 1) - filter_segments(segments); - - switch(segments.size()) - { - case 0: - { - // there can only be one point of contact, otherwise by convexity - // there would be a full segment on a face (interior segment isn't possible either - // because there are at most 2 inside points and an interior segment would also - // yield at least a segment on the boundary) - return result_type(std::forward(points.front())); - } - case 1: // 1 segment - { - const Segment_3& s = segments.front(); - - if(strictly_inside_points == 1) - { - // Doesn't matter whether there is another (non-strictly) inside point: if there is, - // it is an extremity of the segment - - const int str_inside_pt_pos = - int(std::find(vertex_sides.begin(), vertex_sides.end(), ON_BOUNDED_SIDE) - vertex_sides.begin()); - CGAL_assertion(str_inside_pt_pos >= 0 && str_inside_pt_pos < 3); - - Triangle_3 res_tr = triangle(vertex(tr, str_inside_pt_pos), s.source(), s.target()); - return result_type(std::forward(res_tr)); - } - else if(strictly_inside_points == 2) - { - CGAL_assertion(inside_points == 2); // can't be 3 since we're in the 1&2 switch - - Poly res(4); - - // Grab the 2 strictly inside points - int id = 0; - for(int i=0; i<3; ++i) - if(vertex_sides[i] == ON_BOUNDED_SIDE) - res[id++] = vertex(tr, i); - - CGAL_assertion(id == 2); - - if((res[0] - res[1]) * (s.source() - s.target()) > 0) - { - res[2] = s.target(); - res[3] = s.source(); - } - else - { - res[3] = s.target(); - res[2] = s.source(); - } - - return result_type(std::forward(res)); - } - else if(inside_points == 1) // 1 point on the boundary - { - CGAL_assertion(strictly_inside_points == 0); - - // Grab the inside point - const int boundary_pt_pos = - int(std::find(vertex_sides.begin(), vertex_sides.end(), ON_BOUNDARY) - vertex_sides.begin()); - CGAL_assertion(boundary_pt_pos >= 0 && boundary_pt_pos < 3); - - const Point_3& boundary_pt = vertex(tr, boundary_pt_pos); - if(boundary_pt == s.source() || boundary_pt == s.target()) - { - return result_type(s); - } - else - { - Triangle_3 res_tr = triangle(boundary_pt, s.source(), s.target()); - return result_type(std::forward(res_tr)); - } - } - else // 2 points on the boundary - { - CGAL_assertion(inside_points == 2 && strictly_inside_points == 0); - - // 2 boundary points and 1 segment, have to distinguish between cases - // depending on if the extremities of the segment are triangle extremities - - std::array boundary_pts; - std::array is_boundary_point_an_extremity; - - // Grab the inside points - std::size_t id = 0; - for(int i=0; i<3; ++i) - { - if(vertex_sides[i] == ON_BOUNDARY) - { - boundary_pts[id] = i; - - if(vertex(tr, i) == s.source()) - is_boundary_point_an_extremity[id] = true; - else if(vertex(tr, i) == s.target()) - is_boundary_point_an_extremity[id] = true; - else - is_boundary_point_an_extremity[id] = false; - - ++id; - } - } - - CGAL_assertion(id == 2); - - if(is_boundary_point_an_extremity[0]) - { - if(is_boundary_point_an_extremity[1]) - { - // the segment is composed of the two boundary points - return result_type(s); - } - else // only boundary_pts[0] is an extremity - { - Triangle_3 res_tr = triangle(s.source(), s.target(), vertex(tr, boundary_pts[1])); - return result_type(std::forward(res_tr)); - } - } - else // boundary_pts[0] is not an extremity - { - if(is_boundary_point_an_extremity[1]) // only boundary_pts[1] is an extremity - { - Triangle_3 res_tr = triangle(s.source(), s.target(), vertex(tr, boundary_pts[0])); - return result_type(std::forward(res_tr)); - } - else // neither boundary points are extremities - { - Poly res(4); - res[0] = vertex(tr, boundary_pts[0]); - res[1] = vertex(tr, boundary_pts[1]); - - if((res[0] - res[1]) * (s.source() - s.target()) > 0) - { - res[2] = s.target(); - res[3] = s.source(); - } - else - { - res[3] = s.target(); - res[2] = s.source(); - } - - return result_type(std::forward(res)); - } - } - } - - CGAL_assertion(false); - } - break; - // 2 or 3 segments (and 1 or 2 inside points) - case 2: - case 3: - case 4: - { - // @todo do that for a single segment too? - return build_intersection(tet, tr, points, segments, k); - } - break; + case ON_POSITIVE_SIDE: + if (test_segment && orientations[i]==ON_NEGATIVE_SIDE) + current.push_back(*CGAL::Intersections::internal::intersection_point(pl, line(res[i], res[j]), k)); + current.push_back(res[j]); + break; + case ON_NEGATIVE_SIDE: + if (test_segment && orientations[i]==ON_POSITIVE_SIDE) + current.push_back(*CGAL::Intersections::internal::intersection_point(pl, line(res[i], res[j]), k)); + break; default: - // can't have more than 4 segments (1 per tet face) - CGAL_assertion(false); - break; + current.push_back(res[j]); } } - break; - - case 3: - { - // the input triangle is entirely contained within the tetrahedron - return result_type(tr); - } - break; - default: - CGAL_assertion(false); // never happens (only 3 pts in a tr) - break; + res.swap(current); + if (res.empty()) + return boost::none; } - return result_type(); + switch(res.size()) + { + case 1: + return result_type(res[0]); + case 2: + return result_type(segment(res[0], res[1])); + case 3: + return result_type(triangle(res[0], res[1], res[2])); + default: + return result_type(res); + } } template From 445be90c04669b4623530ab8e286ac87a0b38a39 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 1 Aug 2022 13:47:44 +0200 Subject: [PATCH 02/10] update triangle (permutation) --- .../test/Intersections_3/test_intersections_Tetrahedron_3.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp index 6b83de7d349..5aa3a0db2db 100644 --- a/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp +++ b/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp @@ -166,7 +166,7 @@ public: // face inside tr check_intersection(tet, Tr(p(0,2,0), p(0,-2,0), p(0,0,2)), - Tr(p(0,0,1), p(0,0,0), p(0,1,0))); + Tr(p(0,1,0), p(0,0,0), p(0,0,1))); // on the same plane, polygonal intersection Base::template check_intersection(tet, Tr(p(0,-2,0), p(1,1,0), p(1,2,0))); From 5c10048d3e50358361ddbe7537b62d5fc8c73f6e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 1 Aug 2022 13:48:13 +0200 Subject: [PATCH 03/10] add a bitset to track planes containing polygon points used to avoid computing points and use input tetrahedron vertex +workaound for kernels with inexact constructions --- .../Tetrahedron_3_Triangle_3_intersection.h | 67 +++++++++++++++++-- 1 file changed, 62 insertions(+), 5 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h index 263472803ef..5ff59b61852 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h @@ -24,6 +24,7 @@ #include #include #include +#include namespace CGAL { namespace Intersections { @@ -56,17 +57,38 @@ intersection(const typename K::Tetrahedron_3& tet, std::vector res = { vertex(tr,0), vertex(tr,1), vertex(tr,2) }; + std::vector> supporting_planes(3); // bitset used to indicate when a point is on a plane - for (int i=0; i<4; ++i) + // iteratively clip the triangle using the halfspaces which intersections is `tet` + for (int pid=0; pid<4; ++pid) { - Plane_3 pl = plane(vertex(tet, (1+i)%4), vertex(tet, (2+i)%4),vertex(tet, (3+i)%4)); - if (oriented_side(pl, vertex(tet,i))!=ON_POSITIVE_SIDE) // TODO: this should be precomputed based on the tetra orientation + Plane_3 pl = plane(vertex(tet, (1+pid)%4), vertex(tet, (2+pid)%4),vertex(tet, (3+pid)%4)); + if (oriented_side(pl, vertex(tet,pid))!=ON_POSITIVE_SIDE) // TODO: this should be precomputed based on the tetra orientation pl = pl.opposite(); std::vector current; + std::vector> current_sp; std::vector orientations(res.size()); for (std::size_t i=0; i Date: Mon, 1 Aug 2022 13:53:44 +0200 Subject: [PATCH 04/10] add new test from issue --- .../test_intersections_Tetrahedron_3.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp index 5aa3a0db2db..d1815a20b17 100644 --- a/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp +++ b/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp @@ -307,12 +307,23 @@ public: } } + void issue_6777() + { + Tr tri(P(0.191630, -0.331630, -0.370000), P(-0.124185, -0.385815, -0.185000), P(-0.0700000, -0.0700000, 0.00000)); + Tet tet(P(0, -1, 0), P(-1, 0, 0), P(0, 0, 0), P(0, 0, -1)); + auto res = intersection(tri, tet); + assert(res != boost::none); + const std::vector

*vps = boost::get>(&*res); + assert(vps!=nullptr); + } + void run() { std::cout << "3D Tetrahedron Intersection tests\n"; Tet_Tet(); Tet_Tr(); + issue_6777(); } }; From 2193551f42c62b2e768ba904f5a052ea440ff658 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 1 Aug 2022 14:08:22 +0200 Subject: [PATCH 05/10] precompute plane vertex indices --- .../internal/Tetrahedron_3_Triangle_3_intersection.h | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h index 5ff59b61852..832ed0142c8 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h @@ -54,17 +54,20 @@ intersection(const typename K::Tetrahedron_3& tet, typename K::Construct_segment_3 segment = k.construct_segment_3_object(); typename K::Construct_line_3 line = k.construct_line_3_object(); typename K::Oriented_side_3 oriented_side = k.oriented_side_3_object(); + typename K::Orientation_3 orientation = k.orientation_3_object(); std::vector res = { vertex(tr,0), vertex(tr,1), vertex(tr,2) }; std::vector> supporting_planes(3); // bitset used to indicate when a point is on a plane // iteratively clip the triangle using the halfspaces which intersections is `tet` + static constexpr std::array vids = { 1, 2, 3, 0, 3, 2, 0, 1, 3, 1, 0, 2 }; for (int pid=0; pid<4; ++pid) { - Plane_3 pl = plane(vertex(tet, (1+pid)%4), vertex(tet, (2+pid)%4),vertex(tet, (3+pid)%4)); - if (oriented_side(pl, vertex(tet,pid))!=ON_POSITIVE_SIDE) // TODO: this should be precomputed based on the tetra orientation - pl = pl.opposite(); + Plane_3 pl = orientation(tet)==POSITIVE + ? plane(vertex(tet, vids[pid*3]), vertex(tet, vids[pid*3+2]),vertex(tet, vids[pid*3+1])) + : plane(vertex(tet, vids[pid*3]), vertex(tet, vids[pid*3+1]),vertex(tet, vids[pid*3+2])); + CGAL_assertion(oriented_side(pl, vertex(tet,pid))==ON_POSITIVE_SIDE); std::vector current; std::vector> current_sp; From 5d32639dc0bfe55161a1122583f99a9035f6b623 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 1 Aug 2022 14:22:43 +0200 Subject: [PATCH 06/10] fix comment --- .../internal/Tetrahedron_3_Triangle_3_intersection.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h index 832ed0142c8..56aeca15132 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h @@ -60,7 +60,7 @@ intersection(const typename K::Tetrahedron_3& tet, std::vector res = { vertex(tr,0), vertex(tr,1), vertex(tr,2) }; std::vector> supporting_planes(3); // bitset used to indicate when a point is on a plane - // iteratively clip the triangle using the halfspaces which intersections is `tet` + // iteratively clip `tr` with the halfspaces which intersection is `tet` static constexpr std::array vids = { 1, 2, 3, 0, 3, 2, 0, 1, 3, 1, 0, 2 }; for (int pid=0; pid<4; ++pid) { From 77a50cfd7c6a68f2f1b09aedf17df8d490f4a178 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 1 Aug 2022 14:47:25 +0200 Subject: [PATCH 07/10] compute it once --- .../internal/Tetrahedron_3_Triangle_3_intersection.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h index 56aeca15132..82bd3831dd4 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h @@ -62,9 +62,10 @@ intersection(const typename K::Tetrahedron_3& tet, // iteratively clip `tr` with the halfspaces which intersection is `tet` static constexpr std::array vids = { 1, 2, 3, 0, 3, 2, 0, 1, 3, 1, 0, 2 }; + const bool tet_ori_positive = (orientation(tet)==POSITIVE); for (int pid=0; pid<4; ++pid) { - Plane_3 pl = orientation(tet)==POSITIVE + Plane_3 pl = tet_ori_positive ? plane(vertex(tet, vids[pid*3]), vertex(tet, vids[pid*3+2]),vertex(tet, vids[pid*3+1])) : plane(vertex(tet, vids[pid*3]), vertex(tet, vids[pid*3+1]),vertex(tet, vids[pid*3+2])); CGAL_assertion(oriented_side(pl, vertex(tet,pid))==ON_POSITIVE_SIDE); From bf6e2a29299dae4f277f2bea648b8ffff5d365ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 2 Aug 2022 07:32:24 +0200 Subject: [PATCH 08/10] remove unused typedefs --- .../internal/Tetrahedron_3_Triangle_3_intersection.h | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h index 82bd3831dd4..2c39be3faf1 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h @@ -42,11 +42,7 @@ intersection(const typename K::Tetrahedron_3& tet, CGAL_precondition(!tr.is_degenerate()); typedef typename K::Point_3 Point_3; - typedef typename K::Segment_3 Segment_3; typedef typename K::Plane_3 Plane_3; - typedef typename K::Triangle_3 Triangle_3; - typedef typename K::Line_3 Line_3; - typedef typename K::Point_3 Point_3; typename K::Construct_plane_3 plane = k.construct_plane_3_object(); typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object(); From bc9aa692c33b411f46782354ce527f36969dad4e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Wed, 7 Sep 2022 16:00:59 +0200 Subject: [PATCH 09/10] Add a few more tests to Tr-Tet --- .../test_intersections_Tetrahedron_3.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp index d1815a20b17..6f45f6200e6 100644 --- a/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp +++ b/Intersections_3/test/Intersections_3/test_intersections_Tetrahedron_3.cpp @@ -143,6 +143,8 @@ public: // edge shared, 3rd point outside check_intersection(tet, Tr(p(0,1,0), p(1,0,0), P(0.5,0,-100)), S(p(0,1,0), p(1,0,0))); + check_intersection(tet, Tr(P(0.75,0.25,0), p(10,10,10), P(0.25,0.75,0)), + S(P(0.75,0.25,0), P(0.25,0.75,0))); // shared edge, 3rd point inside check_intersection(tet, Tr(p(0,1,0), p(1,0,0), P(0.25,0.25,0.25)), @@ -189,6 +191,14 @@ public: // vertex on edge & triangle inside, double segment non-incident Base::template check_intersection(tet, Tr(P(0.25,0,0.25), P(-1,0.5,0.25), P(1.5,0.5,0.25))); + // vertex on face, triangle outside & point intersection + Base::check_intersection(tet, Tr(P(-1,1,0.25), P(-1,0,0.25), P(0,0.25,0.25)), + P(0, 0.25, 0.25)); + Base::check_intersection(tet, Tr(P(-1,0,0.25), P(-1,1,0.25), P(0,0.25,0.25)), + P(0, 0.25, 0.25)); + Base::check_intersection(tet, Tr(P(0,0.25,0.25), P(-1,1,0.25), P(-1,0,0.25)), + P(0, 0.25, 0.25)); + // vertex on face, triangle outside & segment intersection Base::check_intersection(tet, Tr(P(0.5,0,-0.25), P(0.5,0,0.25), P(0.5,-0.5,0)), S(P(0.5, 0, 0.25), P(0.5, 0, 0))); From d8969404709a33b2a4b3b20381e8410690629b1e Mon Sep 17 00:00:00 2001 From: Sebastien Loriot Date: Wed, 7 Sep 2022 16:08:58 +0200 Subject: [PATCH 10/10] Cosmetic changes Co-authored-by: Mael --- .../internal/Tetrahedron_3_Triangle_3_intersection.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h index 2c39be3faf1..49bc573489e 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Triangle_3_intersection.h @@ -56,8 +56,8 @@ intersection(const typename K::Tetrahedron_3& tet, std::vector res = { vertex(tr,0), vertex(tr,1), vertex(tr,2) }; std::vector> supporting_planes(3); // bitset used to indicate when a point is on a plane - // iteratively clip `tr` with the halfspaces which intersection is `tet` - static constexpr std::array vids = { 1, 2, 3, 0, 3, 2, 0, 1, 3, 1, 0, 2 }; + // iteratively clip `tr` with the halfspaces whose intersection form `tet` + static constexpr std::array vids = { 1,2,3, 0,3,2, 0,1,3, 1,0,2 }; const bool tet_ori_positive = (orientation(tet)==POSITIVE); for (int pid=0; pid<4; ++pid) { @@ -80,11 +80,13 @@ intersection(const typename K::Tetrahedron_3& tet, if (supporting_planes[i].count()==3) { for (int b=0; i<4; ++b) + { if (!supporting_planes[i].test(b)) { res[i] = vertex(tet, b); break; } + } } //-- }