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..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 @@ -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 @@ -16,320 +16,195 @@ #include #include +#include -#include -#include -#include - -#include +#include +#include namespace CGAL { namespace Intersections { namespace internal { -template -void filter_points(std::vector& input, - std::vector& output) -{ - std::sort(input.begin(), input.end()); - auto last = std::unique(input.begin(), input.end()); - input.erase(last, input.end()); - output = std::move(input); -} - -//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; + typename K::Oriented_side_3 oriented_side = k.oriented_side_3_object(); - 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; + const std::array orientations { { + 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]) + } }; - std::vector edges; - edges.reserve(12); + // description of faces of the bbox + 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 } }; - //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)); + 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 } }; - //get all intersections between pl and cub edges - std::vector segments; - std::vector points; + std::array edge_ipt_id; + edge_ipt_id.fill(-1); - for(int i=0; i < 12; ++i) + auto inter_pt_index = + [&plane, &corners, &edge_ipt_id](int i, int j, int edge_id) { - // 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) + if (edge_ipt_id[edge_id]==-1) { - 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); + 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])); } - } - if(segments.empty() && points.empty()) - return result_type(); + return edge_ipt_id[edge_id]; + }; - switch(segments.size()) + bool all_in = true; + bool all_out = true; + + std::vector> neighbor_ids(14, {-1,-1}); + + auto add_neighbor = [&neighbor_ids](int i, int j) { - case 0: - //points dealt with later - break; - case 1: - { - //adj to an edge - if(points.size() == 4) + if (neighbor_ids[i][0] == -1 ) { + neighbor_ids[i][0] = j; + } + else { + if (neighbor_ids[i][0]!=j && neighbor_ids[i][1]==-1) { - return result_type(std::forward(segments.front())); + neighbor_ids[i][1] = j; } - //plane intersecting through an edge (not 2) - 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) + 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) + { + + 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]) + { + case ON_NEGATIVE_SIDE: { - if(!k.equal_3_object()(entry_seg.source(), p) - && ! k.equal_3_object()(entry_seg.target(), p)) + all_out = false; + // check for intersection of the edge + if (orientations[next_id] == ON_POSITIVE_SIDE) { - if(!p1_found) - { - p1 = p; - p1_found = true; - } - else { - p2 = p; - p2_found = true; - break; - } + ids.push_back( + inter_pt_index(current_id, next_id, edge_id)); } - } - 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]); - } - - 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 = 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); - } - - 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) - { - first_found = true; - first_p = p; - } - else - { - tmp_segs.emplace_back(first_p, p); 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); + } } } + + 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]; + add_neighbor(ids[0], ids[1]); + add_neighbor(ids[1], ids[0]); + break; + } + case 1: + solo_id = ids[0]; + default: + break; + } } - if(tmp_segs.size() < 3) - return result_type(); + if (all_in || all_out) return boost::none; + if (start_id == -1) return { result_type(corners[solo_id]) }; - 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 prv_id = -1; + int cur_id = start_id; + std::vector res; + res.reserve(6); + 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 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..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 @@ -32,101 +32,170 @@ 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 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(); - 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: + { + Triangle_3 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 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..ea6dc18cd0f --- /dev/null +++ b/Intersections_3/test/Intersections_3/issue_5428.cpp @@ -0,0 +1,24 @@ +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK; + +typedef EPICK::Point_3 Point; + +int main() +{ + EPICK::Plane_3 pl(Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0)); + + std::vector pts; + 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()); + + auto result = intersection(cub,pl); + + const std::vector* res = boost::get >(&*result); + for (Point p : *res) { + std::cout << p << std::endl; + } + + return 0; +} diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp index 7247a3c472e..f93026b3d11 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 >(&*res); assert(poly != nullptr); assert(poly->size() == 4); @@ -940,8 +942,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 @@ -1100,7 +1102,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 #include #include +#include namespace CGAL{ namespace Polygon_mesh_processing { @@ -48,24 +49,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 @@ -106,19 +89,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) @@ -127,6 +133,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 ]) { @@ -135,13 +142,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; } @@ -152,15 +159,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()); } } } @@ -182,11 +189,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;