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