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..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 @@ -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 @@ -25,260 +24,12 @@ #include #include #include +#include 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 +37,131 @@ 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::Triangle_3 Triangle_3; - typedef std::vector Poly; + typedef typename K::Plane_3 Plane_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(); + typename K::Orientation_3 orientation = k.orientation_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) }; + std::vector> supporting_planes(3); // bitset used to indicate when a point is on a plane - for(int i=0; i<3; ++i) + // 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) { - vertex_sides[i] = bounded_side(tet, vertex(tr, i)); + 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); - if(vertex_sides[i] != ON_UNBOUNDED_SIDE) - ++inside_points; - - if(vertex_sides[i] == ON_BOUNDED_SIDE) + std::vector current; + std::vector> current_sp; + 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) + orientations[i]=oriented_side(pl, res[i]); + if (orientations[i]==ON_ORIENTED_BOUNDARY) { - 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]) + supporting_planes[i].set(pid); + //workaround for kernels with inexact constructions + //-- + if (supporting_planes[i].count()==3) { - // a face is inside the input tr - if(const Triangle_3* t = boost::get(&*intersections[i])) + for (int b=0; i<4; ++b) { - 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 (!supporting_planes[i].test(b)) + { + res[i] = vertex(tet, b); + break; + } } } - } - - 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 + + for (std::size_t i=0; i segments; - 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]) - { - if(const Triangle_3* t = boost::get(&*intersections[i])) + case ON_POSITIVE_SIDE: + if (test_segment && orientations[i]==ON_NEGATIVE_SIDE) { - 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) + current_sp.push_back(supporting_planes[i] & supporting_planes[j]); + current_sp.back().set(pid); + if (current_sp.back().count()==3) { - res[2] = s.target(); - res[3] = s.source(); + for (int b=0; i<4; ++b) + if (!current_sp.back().test(b)) + { + current.push_back(vertex(tet, b)); + break; + } } else - { - res[3] = s.target(); - res[2] = s.source(); - } - - return result_type(std::forward(res)); + current.push_back(*CGAL::Intersections::internal::intersection_point(pl, line(res[i], res[j]), k)); } - else if(inside_points == 1) // 1 point on the boundary + current.push_back(res[j]); + current_sp.push_back(supporting_planes[j]); + break; + case ON_NEGATIVE_SIDE: + if (test_segment && orientations[i]==ON_POSITIVE_SIDE) { - 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()) + current_sp.push_back(supporting_planes[i] & supporting_planes[j]); + current_sp.back().set(pid); + if (current_sp.back().count()==3) { - return result_type(s); + for (int b=0; i<4; ++b) + if (!current_sp.back().test(b)) + { + current.push_back(vertex(tet, b)); + break; + } } else - { - Triangle_3 res_tr = triangle(boundary_pt, s.source(), s.target()); - return result_type(std::forward(res_tr)); - } + current.push_back(*CGAL::Intersections::internal::intersection_point(pl, line(res[i], res[j]), k)); } - 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; + break; default: - // can't have more than 4 segments (1 per tet face) - CGAL_assertion(false); - break; + { + CGAL_assertion(supporting_planes[j].test(pid)); + current.push_back(res[j]); + current_sp.push_back(supporting_planes[j]); + } } } - break; + res.swap(current); + supporting_planes.swap(current_sp); - 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; + 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 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..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)), @@ -166,7 +168,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))); @@ -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))); @@ -307,12 +317,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(); } };