From d71057ce96724f7a12a3c8404613bdd190b5023a Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 5 Feb 2021 12:25:27 +0000 Subject: [PATCH 01/16] Test if there is something on the optional --- .../Iso_cuboid_3_Plane_3_intersection.h | 82 ++++++++++--------- 1 file changed, 45 insertions(+), 37 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h index 977c7dcc584..50562e06dfa 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h @@ -97,11 +97,11 @@ intersection( case 1: { //adj to an edge - if(points.size() == 4) + if(points.size() == 4) // plane outside of cub and just touching at the edge, and hence the endpoints of 4 other edges { return result_type(std::forward(segments.front())); } - //plane intersecting through an edge (not 2) + //plane intersecting through an edge (not 2) This then intersects 2 other edges in points else { Poly res(4); @@ -251,43 +251,51 @@ intersection( CGAL::Plane_3 >::result_type Pl_pl_type; std::vector plane_intersections; - Pl_pl_type pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0), - cub.vertex(1), - cub.vertex(5))); - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0), - cub.vertex(3), - cub.vertex(4))); - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0), - cub.vertex(1), - cub.vertex(3))); - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7), - cub.vertex(6), - cub.vertex(1))); - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7), - cub.vertex(4), - cub.vertex(3))); - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7), - cub.vertex(6), - cub.vertex(4))); - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); + Pl_pl_type pl_inter; + + if(pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0), + cub.vertex(1), + cub.vertex(5)))){ + if(const Line_3* line = boost::get(&*pl_inter)){ + plane_intersections.push_back(*line); + } } + if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0), + cub.vertex(3), + cub.vertex(4)))){ + if(const Line_3* line = boost::get(&*pl_inter)){ + plane_intersections.push_back(*line); + } + } + if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0), + cub.vertex(1), + cub.vertex(3)))){ + if(const Line_3* line = boost::get(&*pl_inter)){ + plane_intersections.push_back(*line); + } + } + if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7), + cub.vertex(6), + cub.vertex(1)))){ + if(const Line_3* line = boost::get(&*pl_inter)){ + plane_intersections.push_back(*line); + } + } + if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7), + cub.vertex(4), + cub.vertex(3)))){ + if(const Line_3* line = boost::get(&*pl_inter)){ + plane_intersections.push_back(*line); + } + } + if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7), + cub.vertex(6), + cub.vertex(4)))){ + if(const Line_3* line = boost::get(&*pl_inter)){ + plane_intersections.push_back(*line); + } + } std::list tmp_segs; for(const auto& line : plane_intersections) { From 6d438a0f945b311f30764dd64510716f2cf4522e Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 9 Feb 2021 15:48:03 +0000 Subject: [PATCH 02/16] WIP --- .../Iso_cuboid_3_Plane_3_intersection.h | 432 ++++++------------ .../test/Intersections_3/issue_5428.cpp | 40 ++ 2 files changed, 186 insertions(+), 286 deletions(-) create mode 100644 Intersections_3/test/Intersections_3/issue_5428.cpp diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h index 50562e06dfa..69e0b6198b9 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h @@ -16,6 +16,7 @@ #include #include +#include #include #include @@ -27,319 +28,178 @@ namespace CGAL { namespace Intersections { namespace internal { -template -void filter_points(std::vector& input, - std::vector& output) +template +int +inter_pt_index(int i, int j, + const Plane_3& plane, + std::vector& points, + std::map, int>& id_map) { - std::sort(input.begin(), input.end()); - auto last = std::unique(input.begin(), input.end()); - input.erase(last, input.end()); - output = std::move(input); + std::pair, int>::iterator, bool> res = + id_map.insert(std::make_pair(make_sorted_pair(i, j), + static_cast (points.size()))); + if (res.second) + points.push_back(typename Geom_traits::Construct_plane_line_intersection_point_3() + (plane, points[i], points[j])); + + return res.first->second; } -//Tetrahedron_3 Line_3 +//Iso_cuboid_3 Plane_3 template typename Intersection_traits::result_type intersection( - const typename K::Iso_cuboid_3 &cub, - const typename K::Plane_3 &pl, - const K& k) + const typename K::Iso_cuboid_3& cub, + const typename K::Plane_3& plane, + const K& k) { typedef typename K::Point_3 Point_3; - typedef typename K::Segment_3 Segment_3; - typedef typename K::Line_3 Line_3; - typedef typename K::Plane_3 Plane_3; typedef std::vector Poly; + typedef typename Intersection_traits::result_type result_type; - typedef typename Intersection_traits, - CGAL::Plane_3 >::result_type result_type; + std::vector corners(8); + corners.reserve(14); // 8 corners + up to 6 polygon points + corners[0] = cub[0]; + corners[1] = cub[3]; + corners[2] = cub[2]; + corners[3] = cub[1]; + corners[4] = cub[5]; + corners[5] = cub[4]; + corners[6] = cub[7]; + corners[7] = cub[6]; - typedef typename Intersection_traits, - CGAL::Plane_3 >::result_type Inter_type; + std::array orientations = { { + plane.oriented_side(corners[0]), + plane.oriented_side(corners[1]), + plane.oriented_side(corners[2]), + plane.oriented_side(corners[3]), + plane.oriented_side(corners[4]), + plane.oriented_side(corners[5]), + plane.oriented_side(corners[6]), + plane.oriented_side(corners[7]) + } }; - std::vector edges; - edges.reserve(12); + // description of faces of the bbox + std::array face_indices = + { { 0, 1, 2, 3, + 2, 1, 5, 6, + 3, 2, 6, 7, + 1, 0, 4, 5, + 4, 0, 3, 7, + 6, 5, 4, 7 } }; - //get all edges of cub - for(int i=0; i< 4; ++i) - edges.emplace_back(cub.vertex(i), cub.vertex((i+1)%4)); - for(int i=0; i < 4; ++i) - edges.emplace_back(cub.vertex(i+4), cub.vertex((i+1)%4+4)); - for(int i=0; i < 4; ++i) - edges.emplace_back(cub.vertex(i), cub.vertex((i+1)%4+4)); + std::map, int> id_map; + bool all_in = true; + bool all_out = true; - //get all intersections between pl and cub edges - std::vector segments; - std::vector points; + std::vector next(14, -1); + std::vector prev(14, -1); - for(int i=0; i < 12; ++i) - { - // Intersect_3 checks the orientation of the segment's extremities to avoid actually computing - // the intersection if possible - Inter_type inter = typename K::Intersect_3()(pl, edges[i]); - if(inter) + int start = -1; + int solo = -1; + // for each face of the bbox, we look for intersection of the plane with its edges + std::vector ids; + for (int i = 0; i < 6; ++i) { - if(const Segment_3* seg = boost::get(&*inter)) - segments.push_back(*seg); - else if(const Point_3* p = boost::get(&*inter)) - points.push_back(*p); - } - } - - if(segments.empty() && points.empty()) - return result_type(); - - switch(segments.size()) - { - case 0: - //points dealt with later - break; - case 1: - { - //adj to an edge - if(points.size() == 4) // plane outside of cub and just touching at the edge, and hence the endpoints of 4 other edges - { - return result_type(std::forward(segments.front())); - } - //plane intersecting through an edge (not 2) This then intersects 2 other edges in points - else - { - Poly res(4); - const Segment_3& entry_seg(segments.front()); - Point_3 p1, p2; - bool p1_found(false), - p2_found(false); - - for(const Point_3& p : points) + ids.clear(); + for (int k = 0; k < 4; ++k) { - if(!k.equal_3_object()(entry_seg.source(), p) - && ! k.equal_3_object()(entry_seg.target(), p)) - { - if(!p1_found) + + int current_id = face_indices[4 * i + k]; + int next_id = face_indices[4 * i + (k + 1) % 4]; + + switch (orientations[current_id]) { - p1 = p; - p1_found = true; + case ON_NEGATIVE_SIDE: + { + all_out = false; + // check for intersection of the edge + if (orientations[next_id] == ON_POSITIVE_SIDE) + { + ids.push_back( + inter_pt_index(current_id, next_id, plane, corners, id_map)); + } + break; + } + case ON_POSITIVE_SIDE: + { + all_in = false; + // check for intersection of the edge + if (orientations[next_id] != ON_NEGATIVE_SIDE) + { + ids.push_back( + inter_pt_index(current_id, next_id, plane, corners, id_map)); + } + break; + } + case ON_ORIENTED_BOUNDARY: + { + all_in = all_out = false; + ids.push_back(current_id); + } } - else { - p2 = p; - p2_found = true; - break; - } - } } - CGAL_USE(p2_found); - CGAL_assertion(p1_found && p2_found); - res[0] = entry_seg.target(); - res[1] = p2; - res[2] = p1; - res[3] = entry_seg.source(); - - typename K::Coplanar_orientation_3 coplanar_orientation = - k.coplanar_orientation_3_object(); - - if( coplanar_orientation(res[0], res[1], res[2]) - != coplanar_orientation(res[0], res[1], res[3])) - { - std::swap(res[1], res[2]); - } - + if (ids.size() == 4){ + std::vector res(4); + res[0] = corners[ids[0]]; + res[1] = corners[ids[1]]; + res[2] = corners[ids[2]]; + res[3] = corners[ids[3]]; return result_type(std::forward(res)); - } - } - break; - - case 2: //intersects diagonally - { - Poly res(4); - Segment_3 &front(segments.front()), - &back(segments.back()); - res[0] = front.target(); - res[1] = back.target(); - res[2] = back.source(); - res[3] = front.source(); - typename K::Coplanar_orientation_3 coplanar_orientation = - k.coplanar_orientation_3_object(); - - if( coplanar_orientation(res[0], res[1], res[2]) - != coplanar_orientation(res[0], res[1], res[3])) - { - std::swap(res[1], res[2]); - } - - return result_type(std::forward(res)); - } - break; - case 4: // intersects a face - { - Poly res; - res.reserve(4); - typename K::Has_on_3 has_on; - if(has_on(pl, cub.vertex(0)) - && has_on(pl, cub.vertex(5)) - && has_on(pl, cub.vertex(4))) - { - res.push_back(cub.vertex(0)); - res.push_back(cub.vertex(5)); - res.push_back(cub.vertex(4)); - res.push_back(cub.vertex(3)); - } - else if(has_on(pl, cub.vertex(0)) - && has_on(pl, cub.vertex(1)) - && has_on(pl, cub.vertex(6))) - { - res.push_back(cub.vertex(0)); - res.push_back(cub.vertex(1)); - res.push_back(cub.vertex(6)); - res.push_back(cub.vertex(5)); - - } - else if(has_on(pl, cub.vertex(1)) - && has_on(pl, cub.vertex(2)) - && has_on(pl, cub.vertex(7))) - { - res.push_back(cub.vertex(1)); - res.push_back(cub.vertex(2)); - res.push_back(cub.vertex(7)); - res.push_back(cub.vertex(6)); - } - else if(has_on(pl, cub.vertex(2)) - && has_on(pl, cub.vertex(3)) - && has_on(pl, cub.vertex(4))) - { - res.push_back(cub.vertex(2)); - res.push_back(cub.vertex(3)); - res.push_back(cub.vertex(4)); - res.push_back(cub.vertex(7)); - } - else if(has_on(pl, cub.vertex(6)) - && has_on(pl, cub.vertex(7)) - && has_on(pl, cub.vertex(4))) - { - res.push_back(cub.vertex(6)); - res.push_back(cub.vertex(7)); - res.push_back(cub.vertex(4)); - res.push_back(cub.vertex(5)); - } - else if(has_on(pl, cub.vertex(2)) - && has_on(pl, cub.vertex(1)) - && has_on(pl, cub.vertex(0))) - { - res.push_back(cub.vertex(2)); - res.push_back(cub.vertex(1)); - res.push_back(cub.vertex(0)); - res.push_back(cub.vertex(3)); - } - return result_type(std::forward(res)); - } - break; - default: - break; - } - - Poly filtered_points; - filter_points(points, filtered_points); - - //adjacent to a vertex - if(filtered_points.size() == 1) - { - return result_type(std::forward(filtered_points.front())); - } - - //get intersections between pl and each face -> line. Foreach line, creates segment with points. Then use helper_function to recover polygon. - typedef typename Intersection_traits, - CGAL::Plane_3 >::result_type Pl_pl_type; - - std::vector plane_intersections; - Pl_pl_type pl_inter; - - if(pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0), - cub.vertex(1), - cub.vertex(5)))){ - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - } - - if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0), - cub.vertex(3), - cub.vertex(4)))){ - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - } - if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0), - cub.vertex(1), - cub.vertex(3)))){ - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - } - if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7), - cub.vertex(6), - cub.vertex(1)))){ - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - } - if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7), - cub.vertex(4), - cub.vertex(3)))){ - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - } - if (pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7), - cub.vertex(6), - cub.vertex(4)))){ - if(const Line_3* line = boost::get(&*pl_inter)){ - plane_intersections.push_back(*line); - } - } - std::list tmp_segs; - for(const auto& line : plane_intersections) - { - bool first_found = false; - Point_3 first_p; - typename K::Has_on_3 has_on; - for(const auto &p : filtered_points) - { - if(has_on(line, p)) - { - if(!first_found) + } else { - first_found = true; - first_p = p; + if (ids.size() == 2) + { + if (start == -1) start = ids[0]; + if (next[ids[0]] == -1) { + next[ids[0]] = ids[1]; + } + else { + prev[ids[0]] = ids[1]; + } + if (next[ids[1]] == -1) { + next[ids[1]] = ids[0]; + } + else { + prev[ids[1]] = ids[0]; + } + + } + else + if (ids.size() == 1) + solo = ids[0]; } - else - { - tmp_segs.emplace_back(first_p, p); - break; - } - } } - } - if(tmp_segs.size() < 3) - return result_type(); + if (all_in || all_out) return boost::none; + if (start == -1) return { result_type(corners[solo]) }; - std::list tmp_pts; - fill_points_list(tmp_segs,tmp_pts); - - Poly res; - for(const auto& p : tmp_pts) - res.push_back(p); - - if(res.size() == 3){ - typename K::Triangle_3 tr(res[0], res[1], res[2]); - return result_type(std::forward(tr)); - } - else - { - return result_type(std::forward(res)); - } + int pre = -1; + int cur = start; + std::vector res; + res.reserve(6); + do { + res.push_back(corners[cur]); + int n = next[cur] == pre ? prev[cur] : next[cur]; + if (n == -1 || n == start){ + if(res.size() == 2){ + typename K::Segment_3 seg(res[0], res[1]); + return result_type(std::forward(seg)); + } + if(res.size() == 3){ + typename K::Triangle_3 tr(res[0], res[1], res[2]); + return result_type(std::forward(tr)); + } + return result_type(std::forward(res));; + } + pre = cur; + cur = n; + } while (true); } + + + template typename Intersection_traits::result_type intersection( @@ -350,6 +210,6 @@ intersection( return intersection(cub, pl, k); } -}}} + }}} #endif // CGAL_INTERNAL_INTERSECTIONS_3_ISO_CUBOID_3_PLANE_3_INTERSECTION_H diff --git a/Intersections_3/test/Intersections_3/issue_5428.cpp b/Intersections_3/test/Intersections_3/issue_5428.cpp new file mode 100644 index 00000000000..7750c069516 --- /dev/null +++ b/Intersections_3/test/Intersections_3/issue_5428.cpp @@ -0,0 +1,40 @@ +#include +#include + +typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK; +typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK; + +typedef CGAL::Cartesian_converter IK_to_EK; + +typedef EPICK::Point_3 Point; + +int main() +{ + IK_to_EK to_exact; + + //EPICK::Plane_3 pl(Point(0,0,10),Point(1,0,10),Point(0,1,10)); + + // EPICK::Plane_3 pl(0.265189, 0.902464, 0.33946, -2.47551); + + EPICK::Plane_3 pl(Point(1, 1, 1), Point(1, 2, 1), Point(1, 2, 2)); + + EPECK::Plane_3 epl = to_exact(pl); + + std::vector pts; + pts.push_back(Point(1,1,1)); + pts.push_back(Point(2,2,2)); + CGAL::Bbox_3 cub = CGAL::bbox_3(pts.begin(), pts.end()); + + EPECK::Iso_cuboid_3 ecub = to_exact(cub); + + auto result = intersection(cub,pl); + + //auto result2 = intersection(ecub,epl); + + const std::vector* res = boost::get >(&*result); + for (Point p : *res) { + std::cout << p << std::endl; + } + + return 0; +} From fb731259d72ada512dfaa955f31007b63ac6095c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 9 Feb 2021 17:51:01 +0100 Subject: [PATCH 03/16] fix comparison --- .../internal/Iso_cuboid_3_Plane_3_intersection.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h index 69e0b6198b9..2cbb0a7f0f2 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h @@ -125,7 +125,7 @@ intersection( { all_in = false; // check for intersection of the edge - if (orientations[next_id] != ON_NEGATIVE_SIDE) + if (orientations[next_id] == ON_NEGATIVE_SIDE) { ids.push_back( inter_pt_index(current_id, next_id, plane, corners, id_map)); From b363977c87e0700914bb10e7d3ad029e60fd17cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 9 Feb 2021 17:58:25 +0100 Subject: [PATCH 04/16] clean up --- .../Iso_cuboid_3_Plane_3_intersection.h | 184 +++++++++--------- 1 file changed, 93 insertions(+), 91 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h index 2cbb0a7f0f2..8ed1fff6400 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h @@ -8,7 +8,7 @@ // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial // // -// Author(s) : Maxime Gimeno +// Author(s) : Sebastien Loriot and Andreas Fabri // #ifndef CGAL_INTERNAL_INTERSECTIONS_3_ISO_CUBOID_3_PLANE_3_INTERSECTION_H @@ -28,7 +28,7 @@ namespace CGAL { namespace Intersections { namespace internal { -template +template int inter_pt_index(int i, int j, const Plane_3& plane, @@ -39,7 +39,7 @@ inter_pt_index(int i, int j, id_map.insert(std::make_pair(make_sorted_pair(i, j), static_cast (points.size()))); if (res.second) - points.push_back(typename Geom_traits::Construct_plane_line_intersection_point_3() + points.push_back(typename K::Construct_plane_line_intersection_point_3() (plane, points[i], points[j])); return res.first->second; @@ -54,8 +54,8 @@ intersection( const K& k) { typedef typename K::Point_3 Point_3; - typedef std::vector Poly; typedef typename Intersection_traits::result_type result_type; + typename K::Oriented_side_3 oriented_side = k.oriented_side_3_object(); std::vector corners(8); corners.reserve(14); // 8 corners + up to 6 polygon points @@ -69,14 +69,14 @@ intersection( corners[7] = cub[6]; std::array orientations = { { - plane.oriented_side(corners[0]), - plane.oriented_side(corners[1]), - plane.oriented_side(corners[2]), - plane.oriented_side(corners[3]), - plane.oriented_side(corners[4]), - plane.oriented_side(corners[5]), - plane.oriented_side(corners[6]), - plane.oriented_side(corners[7]) + oriented_side(plane, corners[0]), + oriented_side(plane, corners[1]), + oriented_side(plane, corners[2]), + oriented_side(plane, corners[3]), + oriented_side(plane, corners[4]), + oriented_side(plane, corners[5]), + oriented_side(plane, corners[6]), + oriented_side(plane, corners[7]) } }; // description of faces of the bbox @@ -95,105 +95,107 @@ intersection( std::vector next(14, -1); std::vector prev(14, -1); - int start = -1; - int solo = -1; + int start_id = -1; + int solo_id = -1; // for each face of the bbox, we look for intersection of the plane with its edges std::vector ids; for (int i = 0; i < 6; ++i) + { + ids.clear(); + for (int k = 0; k < 4; ++k) { - ids.clear(); - for (int k = 0; k < 4; ++k) + + int current_id = face_indices[4 * i + k]; + int next_id = face_indices[4 * i + (k + 1) % 4]; + + switch (orientations[current_id]) + { + case ON_NEGATIVE_SIDE: { - - int current_id = face_indices[4 * i + k]; - int next_id = face_indices[4 * i + (k + 1) % 4]; - - switch (orientations[current_id]) - { - case ON_NEGATIVE_SIDE: - { - all_out = false; - // check for intersection of the edge - if (orientations[next_id] == ON_POSITIVE_SIDE) - { - ids.push_back( - inter_pt_index(current_id, next_id, plane, corners, id_map)); - } - break; - } - case ON_POSITIVE_SIDE: - { - all_in = false; - // check for intersection of the edge - if (orientations[next_id] == ON_NEGATIVE_SIDE) - { - ids.push_back( - inter_pt_index(current_id, next_id, plane, corners, id_map)); - } - break; - } - case ON_ORIENTED_BOUNDARY: - { - all_in = all_out = false; - ids.push_back(current_id); - } - } + all_out = false; + // check for intersection of the edge + if (orientations[next_id] == ON_POSITIVE_SIDE) + { + ids.push_back( + inter_pt_index(current_id, next_id, plane, corners, id_map)); + } + break; } - if (ids.size() == 4){ - std::vector res(4); - res[0] = corners[ids[0]]; - res[1] = corners[ids[1]]; - res[2] = corners[ids[2]]; - res[3] = corners[ids[3]]; - return result_type(std::forward(res)); - } else + case ON_POSITIVE_SIDE: { - if (ids.size() == 2) - { - if (start == -1) start = ids[0]; - if (next[ids[0]] == -1) { - next[ids[0]] = ids[1]; - } - else { - prev[ids[0]] = ids[1]; - } - if (next[ids[1]] == -1) { - next[ids[1]] = ids[0]; - } - else { - prev[ids[1]] = ids[0]; - } - - } - else - if (ids.size() == 1) - solo = ids[0]; + all_in = false; + // check for intersection of the edge + if (orientations[next_id] == ON_NEGATIVE_SIDE) + { + ids.push_back(inter_pt_index(current_id, next_id, plane, corners, id_map)); + } + break; } + case ON_ORIENTED_BOUNDARY: + { + all_in = all_out = false; + ids.push_back(current_id); + } + } } - if (all_in || all_out) return boost::none; - if (start == -1) return { result_type(corners[solo]) }; + switch (ids.size()) + { + case 4: + { + std::vector res({ corners[ids[0]], + corners[ids[1]], + corners[ids[2]], + corners[ids[3]] }); + return result_type(res); + } + case 2: + { + if (start_id == -1) start_id = ids[0]; + if (next[ids[0]] == -1) { + next[ids[0]] = ids[1]; + } + else { + prev[ids[0]] = ids[1]; + } + if (next[ids[1]] == -1) { + next[ids[1]] = ids[0]; + } + else { + prev[ids[1]] = ids[0]; + } + break; + } + case 1: + solo_id = ids[0]; + default: + break; + } + } - int pre = -1; - int cur = start; + if (all_in || all_out) return boost::none; + if (start_id == -1) return { result_type(corners[solo_id]) }; + + int prv_id = -1; + int cur_id = start_id; std::vector res; res.reserve(6); do { - res.push_back(corners[cur]); - int n = next[cur] == pre ? prev[cur] : next[cur]; - if (n == -1 || n == start){ + res.push_back(corners[cur_id]); + int nxt_id = next[cur_id] == prv_id ? prev[cur_id] : next[cur_id]; + if (nxt_id == -1 || nxt_id == start_id){ if(res.size() == 2){ typename K::Segment_3 seg(res[0], res[1]); - return result_type(std::forward(seg)); + return result_type(seg); } if(res.size() == 3){ typename K::Triangle_3 tr(res[0], res[1], res[2]); - return result_type(std::forward(tr)); + return result_type(tr); } - return result_type(std::forward(res));; + return result_type(res); } - pre = cur; - cur = n; + prv_id = cur_id; + cur_id = nxt_id; } while (true); } @@ -210,6 +212,6 @@ intersection( return intersection(cub, pl, k); } - }}} +}}} #endif // CGAL_INTERNAL_INTERSECTIONS_3_ISO_CUBOID_3_PLANE_3_INTERSECTION_H From 95bd626bbdc470d6533a6d73b64b8f26c6de1f16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 9 Feb 2021 18:01:55 +0100 Subject: [PATCH 05/16] use one vector --- .../Iso_cuboid_3_Plane_3_intersection.h | 22 +++++++++---------- 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h index 8ed1fff6400..57d29d2efcc 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h @@ -92,8 +92,7 @@ intersection( bool all_in = true; bool all_out = true; - std::vector next(14, -1); - std::vector prev(14, -1); + std::vector> neighbor_ids(14, {-1,-1}); int start_id = -1; int solo_id = -1; @@ -152,17 +151,17 @@ intersection( case 2: { if (start_id == -1) start_id = ids[0]; - if (next[ids[0]] == -1) { - next[ids[0]] = ids[1]; + if (neighbor_ids[ids[0]][0] == -1) { + neighbor_ids[ids[0]][0] = ids[1]; } else { - prev[ids[0]] = ids[1]; + neighbor_ids[ids[0]][1] = ids[1]; } - if (next[ids[1]] == -1) { - next[ids[1]] = ids[0]; + if (neighbor_ids[ids[1]][0] == -1) { + neighbor_ids[ids[1]][0] = ids[0]; } else { - prev[ids[1]] = ids[0]; + neighbor_ids[ids[1]][1] = ids[0]; } break; } @@ -182,7 +181,9 @@ intersection( res.reserve(6); do { res.push_back(corners[cur_id]); - int nxt_id = next[cur_id] == prv_id ? prev[cur_id] : next[cur_id]; + int nxt_id = neighbor_ids[cur_id][0] == prv_id + ? neighbor_ids[cur_id][1] + : neighbor_ids[cur_id][0]; if (nxt_id == -1 || nxt_id == start_id){ if(res.size() == 2){ typename K::Segment_3 seg(res[0], res[1]); @@ -199,9 +200,6 @@ intersection( } while (true); } - - - template typename Intersection_traits::result_type intersection( From 81a6fe16e835aeec499251e4df6df6d65837f36d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 9 Feb 2021 18:27:58 +0100 Subject: [PATCH 06/16] remove idmap --- .../Iso_cuboid_3_Plane_3_intersection.h | 51 +++++++++++-------- 1 file changed, 29 insertions(+), 22 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h index 57d29d2efcc..dc126efdf8a 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h @@ -28,23 +28,6 @@ namespace CGAL { namespace Intersections { namespace internal { -template -int -inter_pt_index(int i, int j, - const Plane_3& plane, - std::vector& points, - std::map, int>& id_map) -{ - std::pair, int>::iterator, bool> res = - id_map.insert(std::make_pair(make_sorted_pair(i, j), - static_cast (points.size()))); - if (res.second) - points.push_back(typename K::Construct_plane_line_intersection_point_3() - (plane, points[i], points[j])); - - return res.first->second; -} - //Iso_cuboid_3 Plane_3 template typename Intersection_traits::result_type @@ -68,7 +51,7 @@ intersection( corners[6] = cub[7]; corners[7] = cub[6]; - std::array orientations = { { + const std::array orientations { { oriented_side(plane, corners[0]), oriented_side(plane, corners[1]), oriented_side(plane, corners[2]), @@ -80,7 +63,7 @@ intersection( } }; // description of faces of the bbox - std::array face_indices = + constexpr std::array face_indices { { 0, 1, 2, 3, 2, 1, 5, 6, 3, 2, 6, 7, @@ -88,7 +71,30 @@ intersection( 4, 0, 3, 7, 6, 5, 4, 7 } }; - std::map, int> id_map; + constexpr std::array edge_indices + { { 0, 1, 2, 3, + 1, 4, 5, 6, + 2, 6, 7, 8, + 0, 9, 10, 4, + 9, 3, 8, 11, + 5, 10, 11, 7 } }; + + std::array edge_ipt_id; + edge_ipt_id.fill(-1); + + auto inter_pt_index = + [&plane, &corners, &edge_ipt_id](int i, int j, int edge_id) + { + if (edge_ipt_id[edge_id]==-1) + { + edge_ipt_id[edge_id] = static_cast (corners.size()); + corners.push_back(typename K::Construct_plane_line_intersection_point_3() + (plane, corners[i], corners[j])); + } + + return edge_ipt_id[edge_id]; + }; + bool all_in = true; bool all_out = true; @@ -106,6 +112,7 @@ intersection( int current_id = face_indices[4 * i + k]; int next_id = face_indices[4 * i + (k + 1) % 4]; + int edge_id = edge_indices[4 * i + k]; switch (orientations[current_id]) { @@ -116,7 +123,7 @@ intersection( if (orientations[next_id] == ON_POSITIVE_SIDE) { ids.push_back( - inter_pt_index(current_id, next_id, plane, corners, id_map)); + inter_pt_index(current_id, next_id, edge_id)); } break; } @@ -126,7 +133,7 @@ intersection( // check for intersection of the edge if (orientations[next_id] == ON_NEGATIVE_SIDE) { - ids.push_back(inter_pt_index(current_id, next_id, plane, corners, id_map)); + ids.push_back(inter_pt_index(current_id, next_id, edge_id)); } break; } From 81ba9da2929371ab65ed43be33d307e0e5487576 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 9 Feb 2021 18:42:39 +0100 Subject: [PATCH 07/16] remove tabs --- Intersections_3/test/Intersections_3/issue_5428.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Intersections_3/test/Intersections_3/issue_5428.cpp b/Intersections_3/test/Intersections_3/issue_5428.cpp index 7750c069516..844b41943c7 100644 --- a/Intersections_3/test/Intersections_3/issue_5428.cpp +++ b/Intersections_3/test/Intersections_3/issue_5428.cpp @@ -33,7 +33,7 @@ int main() const std::vector* res = boost::get >(&*result); for (Point p : *res) { - std::cout << p << std::endl; + std::cout << p << std::endl; } return 0; From a595e529c547b9d9049b784c81e0e4e16c43f216 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 9 Feb 2021 18:43:05 +0100 Subject: [PATCH 08/16] remove edge map and inter pt set --- .../CGAL/Polygon_mesh_processing/clip.h | 84 +++++++++++-------- 1 file changed, 47 insertions(+), 37 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h index 16a42066d03..74d62e9ab84 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h @@ -33,6 +33,7 @@ #include #include +#include namespace CGAL{ namespace Polygon_mesh_processing { @@ -40,24 +41,6 @@ namespace Polygon_mesh_processing { namespace internal { -template -int -inter_pt_index(int i, int j, - const Plane_3& plane, - std::vector& points, - std::map, int>& id_map) -{ - std::pair, int>::iterator, bool> res = - id_map.insert(std::make_pair(make_sorted_pair(i,j), - static_cast (points.size()))); - if(res.second) - points.push_back( - typename Geom_traits::Construct_plane_line_intersection_point_3() - (plane, points[i], points[j])); - - return res.first->second; -} - template @@ -98,19 +81,42 @@ clip_to_bbox(const Plane_3& plane, }}; // description of faces of the bbox - std::array face_indices = - {{ 0, 1, 2, 3, - 2, 1, 5, 6, - 3, 2, 6, 7, - 1, 0, 4, 5, - 4, 0, 3, 7, - 6, 5, 4, 7 }}; + constexpr std::array face_indices + { { 0, 1, 2, 3, + 2, 1, 5, 6, + 3, 2, 6, 7, + 1, 0, 4, 5, + 4, 0, 3, 7, + 6, 5, 4, 7 } }; + + constexpr std::array edge_indices + { { 0, 1, 2, 3, + 1, 4, 5, 6, + 2, 6, 7, 8, + 0, 9, 10, 4, + 9, 3, 8, 11, + 5, 10, 11, 7 } }; + + std::array edge_ipt_id; + edge_ipt_id.fill(-1); + + auto inter_pt_index = + [&plane, &corners, &edge_ipt_id](int i, int j, int edge_id) + { + if (edge_ipt_id[edge_id]==-1) + { + edge_ipt_id[edge_id] = static_cast (corners.size()); + corners.push_back(typename Geom_traits::Construct_plane_line_intersection_point_3() + (plane, corners[i], corners[j])); + } + + return edge_ipt_id[edge_id]; + }; - std::map, int> id_map; std::vector< std::vector > output_faces(6); bool all_in = true; bool all_out = true; - std::set in_point_ids; // to collect the set of points in the clipped bbox + std::bitset<14> in_point_bits; // to collect the set of points in the clipped bbox // for each face of the bbox, we look for intersection of the plane with its edges for(int i=0; i<6; ++i) @@ -119,6 +125,7 @@ clip_to_bbox(const Plane_3& plane, { int current_id = face_indices[4*i + k]; int next_id = face_indices[4*i + (k+1)%4]; + int edge_id = edge_indices[4 * i + k]; switch(orientations[ current_id ]) { @@ -127,13 +134,13 @@ clip_to_bbox(const Plane_3& plane, all_out=false; // point on or on the negative side output_faces[i].push_back(current_id); - in_point_ids.insert(output_faces[i].back()); + in_point_bits.set(output_faces[i].back()); // check for intersection of the edge if(orientations[ next_id ] == ON_POSITIVE_SIDE) { output_faces[i].push_back( - inter_pt_index(current_id, next_id, plane, corners, id_map)); - in_point_ids.insert(output_faces[i].back()); + inter_pt_index(current_id, next_id, edge_id)); + in_point_bits.set(output_faces[i].back()); } break; } @@ -144,15 +151,15 @@ clip_to_bbox(const Plane_3& plane, if(orientations[ next_id ] == ON_NEGATIVE_SIDE) { output_faces[i].push_back( - inter_pt_index(current_id, next_id, plane, corners, id_map)); - in_point_ids.insert(output_faces[i].back()); + inter_pt_index(current_id, next_id, edge_id)); + in_point_bits.set(output_faces[i].back()); } break; } case ON_ORIENTED_BOUNDARY: { output_faces[i].push_back(current_id); - in_point_ids.insert(output_faces[i].back()); + in_point_bits.set(output_faces[i].back()); } } } @@ -174,11 +181,14 @@ clip_to_bbox(const Plane_3& plane, typedef typename graph_traits::face_descriptor face_descriptor; std::map out_vertices; - for(int i : in_point_ids) + for(int i=0; i<14;++i) { - vertex_descriptor v = add_vertex(tm_out); - out_vertices.insert(std::make_pair(i, v)); - put(vpm_out, v, corners[i]); + if (in_point_bits.test(i)) + { + vertex_descriptor v = add_vertex(tm_out); + out_vertices.insert(std::make_pair(i, v)); + put(vpm_out, v, corners[i]); + } } std::map< std::pair, halfedge_descriptor> hedge_map; From 38142dc101ab984c676a95dccd44d1956639af8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 9 Feb 2021 18:51:38 +0100 Subject: [PATCH 09/16] update include directives --- .../internal/Iso_cuboid_3_Plane_3_intersection.h | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h index dc126efdf8a..311c397ddfa 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h @@ -18,11 +18,8 @@ #include #include -#include -#include -#include - -#include +#include +#include namespace CGAL { namespace Intersections { From 231e9f1b9321d0c6bab3595e2b0e79c504ef6a4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 10 Feb 2021 10:29:02 +0100 Subject: [PATCH 10/16] handle cases when plane contains an edge --- .../Iso_cuboid_3_Plane_3_intersection.h | 27 ++++++++++--------- .../Intersections_3/test_intersections_3.cpp | 7 ++--- 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h index 311c397ddfa..270693f2fb9 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h @@ -97,6 +97,19 @@ intersection( std::vector> neighbor_ids(14, {-1,-1}); + auto add_neighbor = [&neighbor_ids](int i, int j) + { + if (neighbor_ids[i][0] == -1 ) { + neighbor_ids[i][0] = j; + } + else { + if (neighbor_ids[i][0]!=j && neighbor_ids[i][1]==-1) + { + neighbor_ids[i][1] = j; + } + } + }; + int start_id = -1; int solo_id = -1; // for each face of the bbox, we look for intersection of the plane with its edges @@ -155,18 +168,8 @@ intersection( case 2: { if (start_id == -1) start_id = ids[0]; - if (neighbor_ids[ids[0]][0] == -1) { - neighbor_ids[ids[0]][0] = ids[1]; - } - else { - neighbor_ids[ids[0]][1] = ids[1]; - } - if (neighbor_ids[ids[1]][0] == -1) { - neighbor_ids[ids[1]][0] = ids[0]; - } - else { - neighbor_ids[ids[1]][1] = ids[0]; - } + add_neighbor(ids[0], ids[1]); + add_neighbor(ids[1], ids[0]); break; } case 1: diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp index 7247a3c472e..ea68b96e871 100644 --- a/Intersections_3/test/Intersections_3/test_intersections_3.cpp +++ b/Intersections_3/test/Intersections_3/test_intersections_3.cpp @@ -907,7 +907,7 @@ struct Test { //edge check_intersection (cub, Pl(P(1,1,1), P(1,2,1), P(1.5,0,0)), - S(P(1,2,1), P(1,1,1))); + S(P(1,1,1), P(1,2,1))); //face @@ -924,6 +924,7 @@ struct Test { assert(p.x() == 1); } res = CGAL::intersection(cub, Pl(P(1,1,1), P(1,2,1), P(2,2,2))); + poly = boost::get >(&*res); assert(poly != nullptr); assert(poly->size() == 4); @@ -940,8 +941,8 @@ struct Test { check_intersection (cub, Pl(P(2, 1.66, 2), P(1.66,2,2), P(2,2,1.66)), - Tr(P(2, 2, 1.66), - P(1.66,2,2), + Tr(P(1.66,2,2), + P(2, 2, 1.66), P(2,1.66,2))); //other edge From 482db1f0ccbb2c28330da25ffa3bd604c555cfcf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 10 Feb 2021 11:11:01 +0100 Subject: [PATCH 11/16] forget to update bbox case --- .../test/Intersections_3/test_intersections_3.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp index ea68b96e871..deb48c9913f 100644 --- a/Intersections_3/test/Intersections_3/test_intersections_3.cpp +++ b/Intersections_3/test/Intersections_3/test_intersections_3.cpp @@ -1101,7 +1101,7 @@ struct Test { //edge check_intersection (cub, Pl(P(1,1,1), P(1,2,1), P(1.5,0,0)), - S(P(1,2,1), P(1,1,1))); + S(P(1,1,1), P(1,2,1))); //face typedef typename CGAL::Intersection_traits Date: Wed, 10 Feb 2021 15:41:45 +0000 Subject: [PATCH 12/16] WIP: try to do plane/tet intersection. I think the edge indexing is wrong --- .../Tetrahedron_3_Plane_3_intersection.h | 228 ++++++++++++------ 1 file changed, 148 insertions(+), 80 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Plane_3_intersection.h index fa55cb268f2..9bc5c0916a0 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Plane_3_intersection.h @@ -32,101 +32,169 @@ template typename Intersection_traits::result_type intersection( const typename K::Tetrahedron_3 &tet, - const typename K::Plane_3 &pl, - const K&) + const typename K::Plane_3 &plane, + const K& k) { - typedef typename Intersection_traits::result_type result_type; + typedef typename K::Point_3 Point_3; + typedef typename Intersection_traits::result_type result_type; + typename K::Oriented_side_3 oriented_side = k.oriented_side_3_object(); - typedef typename Intersection_traits::result_type Inter_type; + std::vector corners(4); + corners.reserve(8); // 4 corners + up to 4 polygon points + corners[0] = tet[0]; + corners[1] = tet[1]; + corners[2] = tet[2]; + corners[3] = tet[3]; - typedef typename K::Segment_3 Segment_3; - Inter_type intersections[4]; - int p_id = -1; - std::vector points; - std::vector segments; - for(int i = 0; i < 4; ++i) + const std::array orientations { { + oriented_side(plane, corners[0]), + oriented_side(plane, corners[1]), + oriented_side(plane, corners[2]), + oriented_side(plane, corners[3]) + } }; + + // description of faces of the bbox + constexpr std::array face_indices + { { 0, 1, 2, + 0, 1, 3, + 1, 2, 3, + 2, 0, 3 } }; + + constexpr std::array edge_indices + { { 0, 1, 2, + 0, 3, 5, + 1, 4, 3, + 2, 5, 4 } }; + + std::array edge_ipt_id; + edge_ipt_id.fill(-1); + + auto inter_pt_index = + [&plane, &corners, &edge_ipt_id](int i, int j, int edge_id) { - const typename K::Triangle_3 triangle(tet.vertex((i+1)%4), - tet.vertex((i+2)%4), - tet.vertex((i+3)%4)); - intersections[i] = typename K::Intersect_3()(pl, triangle); - if(intersections[i]){ - if(const typename K::Triangle_3* tr = boost::get(&*intersections[i])) + if (edge_ipt_id[edge_id]==-1) + { + edge_ipt_id[edge_id] = static_cast (corners.size()); + corners.push_back(typename K::Construct_plane_line_intersection_point_3() + (plane, corners[i], corners[j])); + } + + return edge_ipt_id[edge_id]; + }; + + bool all_in = true; + bool all_out = true; + + std::vector> neighbor_ids(8, {-1,-1}); + + auto add_neighbor = [&neighbor_ids](int i, int j) + { + if (neighbor_ids[i][0] == -1 ) { + neighbor_ids[i][0] = j; + } + else { + if (neighbor_ids[i][0]!=j && neighbor_ids[i][1]==-1) { - typename K::Triangle_3 res = *tr; - return result_type(std::forward(res)); - } - else if( const Segment_3* s - = boost::get(&*intersections[i])) - { - segments.push_back(*s); - } - else if( const typename K::Point_3* p - = boost::get(&*intersections[i])) - { - points.push_back(*p); - p_id = i; + neighbor_ids[i][1] = j; } } - } - CGAL_assertion(segments.size() != 1); + }; - switch(segments.size()) + int start_id = -1; + int solo_id = -1; + // for each face of the bbox, we look for intersection of the plane with its edges + std::vector ids; + for (int i = 0; i < 4; ++i) { - case 0: - { - if(p_id == -1) - return result_type(); - else + ids.clear(); + for (int k = 0; k < 3; ++k) { - typename K::Point_3 p - = *boost::get(&*intersections[p_id]); - return result_type(std::forward(p)); + int current_id = face_indices[3 * i + k]; + int next_id = face_indices[3 * i + (k + 1) % 3]; + int edge_id = edge_indices[3 * i + k]; + + switch (orientations[current_id]) + { + case ON_NEGATIVE_SIDE: + { + all_out = false; + // check for intersection of the edge + if (orientations[next_id] == ON_POSITIVE_SIDE) + { + ids.push_back( + inter_pt_index(current_id, next_id, edge_id)); + } + break; + } + case ON_POSITIVE_SIDE: + { + all_in = false; + // check for intersection of the edge + if (orientations[next_id] == ON_NEGATIVE_SIDE) + { + ids.push_back(inter_pt_index(current_id, next_id, edge_id)); + } + break; + } + case ON_ORIENTED_BOUNDARY: + { + all_in = all_out = false; + ids.push_back(current_id); + } + } } - } - break; - case 2: - { - return result_type(std::forward(segments.back())); - } - break; - default: - { - std::set all_points; - for (const auto& s : segments) + + switch (ids.size()) { - all_points.insert(s.source()); - all_points.insert(s.target()); - } - if(all_points.size() == 3) - { - auto p_it = all_points.begin(); - ++p_it; - typename K::Point_3 mid_point = *p_it; - ++p_it; - typename K::Triangle_3 result(*all_points.begin(), mid_point, *p_it ); - return result_type(std::forward(result)); - } - else //size = 4 - { - std::list segs(segments.begin(), segments.end()); - std::list tmp; - fill_points_list(segs, tmp); - std::vector res; - for( const auto& p : tmp) - res.push_back(p); - return result_type(std::forward >(res)); + case 3: + { + std::vector res({ corners[ids[0]], + corners[ids[1]], + corners[ids[2]] }); + return result_type(res); + } + case 2: + { + if (start_id == -1) start_id = ids[0]; + add_neighbor(ids[0], ids[1]); + add_neighbor(ids[1], ids[0]); + break; + } + case 1: + solo_id = ids[0]; + default: + break; } } - break; - } - CGAL_assertion(false); - return result_type(); + + if (all_in || all_out) return boost::none; + if (start_id == -1) return { result_type(corners[solo_id]) }; + + int prv_id = -1; + int cur_id = start_id; + std::vector res; + res.reserve(4); + do { + res.push_back(corners[cur_id]); + int nxt_id = neighbor_ids[cur_id][0] == prv_id + ? neighbor_ids[cur_id][1] + : neighbor_ids[cur_id][0]; + if (nxt_id == -1 || nxt_id == start_id){ + if(res.size() == 2){ + typename K::Segment_3 seg(res[0], res[1]); + return result_type(seg); + } + if(res.size() == 3){ + typename K::Triangle_3 tr(res[0], res[1], res[2]); + return result_type(tr); + } + return result_type(res); + } + prv_id = cur_id; + cur_id = nxt_id; + } while (true); + } template From 4bc161a332edf22fdf9643c7d7d948bd14b1759f Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 12 Feb 2021 09:42:52 +0000 Subject: [PATCH 13/16] Return a triangle not a vector --- .../internal/Tetrahedron_3_Plane_3_intersection.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Plane_3_intersection.h index 9bc5c0916a0..498ab69a56e 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Plane_3_intersection.h @@ -36,6 +36,7 @@ intersection( const K& k) { typedef typename K::Point_3 Point_3; + typedef typename K::Triangle_3 Triangle_3; typedef typename Intersection_traits::result_type result_type; typename K::Oriented_side_3 oriented_side = k.oriented_side_3_object(); @@ -149,9 +150,9 @@ intersection( { case 3: { - std::vector res({ corners[ids[0]], - corners[ids[1]], - corners[ids[2]] }); + Triangle_3 res(corners[ids[0]], + corners[ids[1]], + corners[ids[2]]); return result_type(res); } case 2: From e833e86598aec6f4a536165573749deb317c1ab4 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 12 Feb 2021 09:43:14 +0000 Subject: [PATCH 14/16] Change the expected result (a permutation) --- Intersections_3/test/Intersections_3/test_intersections_3.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp index deb48c9913f..97b895b9c1c 100644 --- a/Intersections_3/test/Intersections_3/test_intersections_3.cpp +++ b/Intersections_3/test/Intersections_3/test_intersections_3.cpp @@ -672,7 +672,8 @@ struct Test { P(0,1,0)); check_intersection (tet, Pl(P(0,0.5,0), P(1,0.5,-5), P(0.5,0.5,0.5)), - Tr(P(0.5,0.5,0), P(0,0.5,0), P(0,0.5,0.5))); + Tr(P(0, 0.5, 0), P(0.5,0.5,0), P(0,0.5,0.5))); + Pl pl(P(0,0.9,0), P(0.9,0,0), P(0.9,0.01,0.06)); typedef typename CGAL::Intersection_traits Date: Fri, 12 Feb 2021 10:36:53 +0000 Subject: [PATCH 15/16] remove trailing whitespace --- Intersections_3/test/Intersections_3/test_intersections_3.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp index 97b895b9c1c..f93026b3d11 100644 --- a/Intersections_3/test/Intersections_3/test_intersections_3.cpp +++ b/Intersections_3/test/Intersections_3/test_intersections_3.cpp @@ -673,7 +673,7 @@ struct Test { check_intersection (tet, Pl(P(0,0.5,0), P(1,0.5,-5), P(0.5,0.5,0.5)), Tr(P(0, 0.5, 0), P(0.5,0.5,0), P(0,0.5,0.5))); - + Pl pl(P(0,0.9,0), P(0.9,0,0), P(0.9,0.01,0.06)); typedef typename CGAL::Intersection_traits Date: Fri, 12 Feb 2021 14:06:20 +0000 Subject: [PATCH 16/16] Restore issue to the initial bug report --- .../test/Intersections_3/issue_5428.cpp | 22 +++---------------- 1 file changed, 3 insertions(+), 19 deletions(-) diff --git a/Intersections_3/test/Intersections_3/issue_5428.cpp b/Intersections_3/test/Intersections_3/issue_5428.cpp index 844b41943c7..ea6dc18cd0f 100644 --- a/Intersections_3/test/Intersections_3/issue_5428.cpp +++ b/Intersections_3/test/Intersections_3/issue_5428.cpp @@ -1,36 +1,20 @@ #include -#include -typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK; typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK; -typedef CGAL::Cartesian_converter IK_to_EK; - typedef EPICK::Point_3 Point; int main() { - IK_to_EK to_exact; - - //EPICK::Plane_3 pl(Point(0,0,10),Point(1,0,10),Point(0,1,10)); - - // EPICK::Plane_3 pl(0.265189, 0.902464, 0.33946, -2.47551); - - EPICK::Plane_3 pl(Point(1, 1, 1), Point(1, 2, 1), Point(1, 2, 2)); - - EPECK::Plane_3 epl = to_exact(pl); + EPICK::Plane_3 pl(Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0)); std::vector pts; - pts.push_back(Point(1,1,1)); - pts.push_back(Point(2,2,2)); + pts.push_back(Point(-10,-10,-10)); + pts.push_back(Point(10,10,10)); CGAL::Bbox_3 cub = CGAL::bbox_3(pts.begin(), pts.end()); - EPECK::Iso_cuboid_3 ecub = to_exact(cub); - auto result = intersection(cub,pl); - //auto result2 = intersection(ecub,epl); - const std::vector* res = boost::get >(&*result); for (Point p : *res) { std::cout << p << std::endl;