From 3b003535c717960c6038ef172d6b4a494dbbf0d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 17 Jan 2024 10:31:19 +0100 Subject: [PATCH 01/15] WIP --- .../Polygon_mesh_processing/autorefinement.h | 421 ++++++++++++------ 1 file changed, 274 insertions(+), 147 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index a3606f816a7b..1cb0e3aaa4d8 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -114,7 +114,6 @@ Segment_inter_type do_coplanar_segments_intersect(std::size_t pi, std::size_t qi, std::size_t ri, std::size_t si, const std::vector& points, - const typename K::Vector_3& /* plane_normal */, const K& k = K()) { typename K::Collinear_are_ordered_along_line_3 cln_order = k.collinear_are_ordered_along_line_3_object(); @@ -249,11 +248,13 @@ do_coplanar_segments_intersect(std::size_t pi, std::size_t qi, return NO_INTERSECTION; } + // imported from Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h +#warning TODO fill the indices correctly template void coplanar_intersections(const std::array& t1, const std::array& t2, - std::vector& inter_pts) + std::vector>& inter_pts) { const typename K::Point_3& p1 = t1[0], q1 = t1[1], r1 = t1[2]; const typename K::Point_3& p2 = t2[0], q2 = t2[1], r2 = t2[2]; @@ -270,15 +271,14 @@ void coplanar_intersections(const std::array& t1, Intersections::internal::intersection_coplanar_triangles_cutoff(r1,p1,q1,2,p2,q2,r2,k,l_inter_pts); //line r1p1 for (const Intersections::internal::Point_on_triangle& pot : l_inter_pts) - inter_pts.push_back( pot.point(p1,q1,r1,p2,q2,r2,k) ); + inter_pts.emplace_back(pot.point(p1,q1,r1,p2,q2,r2,k), 0, 0); } -// imported from Polygon_mesh_processing/internal/Corefinement/intersect_triangle_segment_3.h +// imported from Polygon_mesh_processing/internal/Corefinement/intersect_triangle_and_segment_3.h template -void +std::optional> find_intersection(const typename K::Point_3& p, const typename K::Point_3& q, //segment const typename K::Point_3& a, const typename K::Point_3& b, const typename K::Point_3& c, //triangle - std::vector& inter_pts, bool is_p_coplanar=false, bool is_q_coplanar=false) // note that in coref this was wrt a halfedge not p/q { Orientation ab=orientation(p,q,a,b); @@ -286,54 +286,56 @@ find_intersection(const typename K::Point_3& p, const typename K::Point_3& q, / Orientation ca=orientation(p,q,c,a); if ( ab==POSITIVE || bc==POSITIVE || ca==POSITIVE ) - return; + return std::nullopt; int nb_coplanar=(ab==COPLANAR?1:0) + (bc==COPLANAR?1:0) + (ca==COPLANAR?1:0); - if (is_p_coplanar) - { - inter_pts.push_back(p); - return; - } - if (is_q_coplanar) - { - inter_pts.push_back(q); - return; - } - - if (nb_coplanar!=2) - { - inter_pts.push_back( - typename K::Construct_plane_line_intersection_point_3()(a, b, c, p, q) - ); - } - else + if (nb_coplanar==2) { + // even if a common point is not new it is still needed to be reported so + // that the intersection segment is known. if (ab!=COPLANAR) - { // intersection is c - inter_pts.push_back(c); - return; - } - - if (bc!=COPLANAR) - { + return std::make_pair(c, -3); + else if (bc!=COPLANAR) // intersection is a - inter_pts.push_back(a); - return; + return std::make_pair(a, -1); + else + { + CGAL_assertion(ca!=COPLANAR); + // intersection is b + return std::make_pair(b, -2); } - CGAL_assertion(ca!=COPLANAR); - // intersection is b - inter_pts.push_back(b); } + + typename K::Point_3 ipt = is_p_coplanar ? p : + is_q_coplanar ? q : + typename K::Construct_plane_line_intersection_point_3() + (a, b, c, p, q); + + if (nb_coplanar == 0) + return std::make_pair(ipt, 0); + + + CGAL_assertion(nb_coplanar==1); + + if (ab==COPLANAR) + // intersection is ab + return std::make_pair(ipt, 1); + if (bc==COPLANAR) + // intersection is bc + return std::make_pair(ipt, 2); + CGAL_assertion(ca==COPLANAR); + // intersection is ca + return std::make_pair(ipt, 3); } template -void test_edge(const typename K::Point_3& p, const typename K::Point_3& q, - const typename K::Point_3& a, const typename K::Point_3& b, const typename K::Point_3& c, - const Orientation abcp, - const Orientation abcq, - std::vector& inter_pts) +std::optional> +test_edge(const typename K::Point_3& p, const typename K::Point_3& q, + const typename K::Point_3& a, const typename K::Point_3& b, const typename K::Point_3& c, + const Orientation abcp, + const Orientation abcq) { switch ( abcp ) { case POSITIVE: @@ -341,46 +343,40 @@ void test_edge(const typename K::Point_3& p, const typename K::Point_3& q, case POSITIVE: // the segment lies in the positive open halfspaces defined by the // triangle's supporting plane - break; + return std::nullopt; case NEGATIVE: // p sees the triangle in counterclockwise order - find_intersection(p,q,a,b,c,inter_pts); - break; + return find_intersection(p,q,a,b,c); //case COPLANAR: default: // q belongs to the triangle's supporting plane // p sees the triangle in counterclockwise order - find_intersection(p,q,a,b,c,inter_pts,false,true); + return find_intersection(p,q,a,b,c,false,true); } - break; case NEGATIVE: switch ( abcq ) { case POSITIVE: // q sees the triangle in counterclockwise order - find_intersection(q,p,a,b,c,inter_pts); - break; + return find_intersection(q,p,a,b,c); case NEGATIVE: // the segment lies in the negative open halfspaces defined by the // triangle's supporting plane - break; + return std::nullopt; // case COPLANAR: default: // q belongs to the triangle's supporting plane // p sees the triangle in clockwise order - find_intersection(q,p,a,b,c,inter_pts,true,false); + return find_intersection(q,p,a,b,c,true,false); } - break; default: //case COPLANAR: // p belongs to the triangle's supporting plane switch ( abcq ) { case POSITIVE: // q sees the triangle in counterclockwise order - find_intersection(q,p,a,b,c,inter_pts,false, true); - break; + return find_intersection(q,p,a,b,c,false, true); case NEGATIVE: // q sees the triangle in clockwise order - find_intersection(p,q,a,b,c,inter_pts,true); - break; + return find_intersection(p,q,a,b,c,true); //case COPLANAR: default: // the segment is coplanar with the triangle's supporting plane @@ -392,7 +388,7 @@ void test_edge(const typename K::Point_3& p, const typename K::Point_3& q, // nothing done as coplanar case handle in collect_intersections // and other intersection points will be collected with non-coplanar edges //} - break; + return std::nullopt; } } } @@ -400,7 +396,7 @@ void test_edge(const typename K::Point_3& p, const typename K::Point_3& q, template bool collect_intersections(const std::array& t1, const std::array& t2, - std::vector& inter_pts) + std::vector>& inter_pts) { // test edges of t1 vs t2 std::array ori; @@ -409,6 +405,7 @@ bool collect_intersections(const std::array& t1, if (ori[0]== COPLANAR && ori[1]==COPLANAR && ori[2]==COPLANAR) { + coplanar_intersections(t1, t2, inter_pts); #ifdef CGAL_AUTOREF_DEBUG_DEPTH for (auto p : inter_pts) @@ -421,7 +418,17 @@ bool collect_intersections(const std::array& t1, for (int i=0; i<3; ++i) { int j=(i+1)%3; - test_edge(t1[i], t1[j], t2[0], t2[1], t2[2], ori[i], ori[j], inter_pts); + std::optional > opt = + test_edge(t1[i], t1[j], t2[0], t2[1], t2[2], ori[i], ori[j]); + if (opt) + { + if (ori[i]==COPLANAR) + inter_pts.emplace_back(opt->first, -(i+1), opt->second); + else if (ori[j]==COPLANAR) + inter_pts.emplace_back(opt->first, -(j+1), opt->second); + else + inter_pts.emplace_back(opt->first, i+1 , opt->second); + } } // test edges of t2 vs t1 @@ -430,12 +437,23 @@ bool collect_intersections(const std::array& t1, for (int i=0; i<3; ++i) { int j=(i+1)%3; - test_edge(t2[i], t2[j], t1[0], t1[1], t1[2], ori[i], ori[j], inter_pts); + std::optional > opt = + test_edge(t2[i], t2[j], t1[0], t1[1], t1[2], ori[i], ori[j]); + if (opt) + { + if (ori[i]==COPLANAR) + inter_pts.emplace_back(opt->first, opt->second, -(i+1)); + else if (ori[j]==COPLANAR) + inter_pts.emplace_back(opt->first, opt->second, -(j+1)); + else + inter_pts.emplace_back(opt->first, opt->second, i+1 ); + } } + #warning TODO get rid of sort and unique calls // because we don't handle intersection type and can have edge-edge edge-vertex duplicates - std::sort(inter_pts.begin(), inter_pts.end()); - auto last = std::unique(inter_pts.begin(), inter_pts.end()); + std::sort(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return get<0>(p)(q);}); + auto last = std::unique(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return get<0>(p)==get<0>(q);}); inter_pts.erase(last, inter_pts.end()); #ifdef CGAL_AUTOREF_DEBUG_DEPTH @@ -446,45 +464,24 @@ bool collect_intersections(const std::array& t1, return false; } +//TODO: rename struct +struct Triangle_data +{ + std::vector points; + std::vector> segments; + std::vector segment_input_triangle_ids; +}; + template void generate_subtriangles(std::size_t ti, - std::vector>& segments, - std::vector& points, - const std::vector& in_triangle_ids, + Triangle_data& triangle_data, const std::set >& intersecting_triangles, const std::set >& coplanar_triangles, const std::vector>& triangles, PointVector& new_triangles ) { - typedef CGAL::Projection_traits_3 P_traits; - typedef CGAL::No_constraint_intersection_tag Itag; - - typedef CGAL::Constrained_Delaunay_triangulation_2 CDT_2; - //typedef CGAL::Constrained_triangulation_plus_2 CDT; - typedef CDT_2 CDT; - - const std::array& t = triangles[ti]; - - // positive triangle normal - typename EK::Vector_3 n = normal(t[0], t[1], t[2]); - typename EK::Point_3 o(CGAL::ORIGIN); - - bool orientation_flipped = false; - if ( typename EK::Less_xyz_3()(o+n,o) ) - { - n=-n; - orientation_flipped = true; - } - - P_traits cdt_traits(n); - CDT cdt(cdt_traits); - cdt.insert_outside_affine_hull(t[0]); - cdt.insert_outside_affine_hull(t[1]); - typename CDT::Vertex_handle v = cdt.tds().insert_dim_up(cdt.infinite_vertex(), orientation_flipped); - v->set_point(t[2]); - #ifdef CGAL_AUTOREFINE_DEBUG_COUNTERS struct Counter { @@ -518,7 +515,13 @@ void generate_subtriangles(std::size_t ti, #define CGAL_AUTOREF_COUNTER_INSTRUCTION(X) #endif - // pre-compute segment intersections + #warning TODO + + std::vector& points=triangle_data.points; + std::vector>& segments=triangle_data.segments; + std::vector& in_triangle_ids=triangle_data.segment_input_triangle_ids; + +// pre-compute segment intersections if (!segments.empty()) { std::size_t nbs = segments.size(); @@ -528,7 +531,7 @@ void generate_subtriangles(std::size_t ti, CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer1.start();) std::map point_id_map; - +//TODO: we already have sorted the points while deduplicating segments! for (std::size_t pid=0; pid(segments[i].first, segments[i].second, segments[j].first, segments[j].second, - points, n); + points); CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer5.stop();) switch(seg_inter_type) @@ -739,13 +742,63 @@ void generate_subtriangles(std::size_t ti, auto last = std::unique(segments.begin(), segments.end()); segments.erase(last, segments.end()); +// init CDT + insert points and constraints CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer3.start();) - if (segments.empty()) - cdt.insert(points.begin(), points.end()); - else - cdt.insert_constraints(points.begin(), points.end(), segments.begin(), segments.end()); + typedef CGAL::Projection_traits_3 P_traits; + typedef CGAL::No_constraint_intersection_tag Itag; + + typedef CGAL::Constrained_Delaunay_triangulation_2 CDT_2; + //typedef CGAL::Constrained_triangulation_plus_2 CDT; + typedef CDT_2 CDT; + + const std::array& t = triangles[ti]; + + // positive triangle normal + typename EK::Vector_3 n = normal(t[0], t[1], t[2]); + typename EK::Point_3 o(CGAL::ORIGIN); + + bool orientation_flipped = false; + if ( typename EK::Less_xyz_3()(o+n,o) ) + { + n=-n; + orientation_flipped = true; + } + + std::vector vhandles(triangle_data.points.size()); + P_traits cdt_traits(n); + CDT cdt(cdt_traits); + vhandles[0]=cdt.insert_outside_affine_hull(t[0]); + vhandles[1]=cdt.insert_outside_affine_hull(t[1]); + vhandles[2] = cdt.tds().insert_dim_up(cdt.infinite_vertex(), orientation_flipped); + vhandles[2]->set_point(t[2]); + + // insert points and fill vhandles + std::vector indices(triangle_data.points.size()-3); + std::iota(indices.begin(), indices.end(), 3); + typedef typename Pointer_property_map::type Pmap; + typedef Spatial_sort_traits_adapter_2 Search_traits; + spatial_sort(indices.begin(), indices.end(), + Search_traits(make_property_map(points), cdt_traits)); + + typename CDT::Face_handle hint; + for (std::size_t i : indices) + { + vhandles[i] = cdt.insert(points[i], hint); + hint=vhandles[i]->face(); + } + + for (const std::pair& ids : triangle_data.segments) + { + //TODO remove me + CGAL_assertion(ids.first < vhandles.size()); + CGAL_assertion(ids.second < vhandles.size()); + CGAL_assertion( vhandles[ids.first]!= typename CDT::Vertex_handle() ); + CGAL_assertion( vhandles[ids.second]!= typename CDT::Vertex_handle() ); + cdt.insert_constraint(vhandles[ids.first], vhandles[ids.second]); + } CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer3.stop();) +// extract new triangles for (typename CDT::Face_handle fh : cdt.finite_face_handles()) { if (orientation_flipped) @@ -838,7 +891,8 @@ void autorefine_triangle_soup(PointRange& soup_points, Visitor visitor(choose_parameter(get_parameter(np, internal_np::visitor))); - constexpr bool parallel_execution = std::is_same_v; + //~ constexpr bool parallel_execution = std::is_same_v; + constexpr bool parallel_execution = false; #ifndef CGAL_LINKED_WITH_TBB static_assert (!parallel_execution, @@ -895,21 +949,24 @@ void autorefine_triangle_soup(PointRange& soup_points, // init the vector of triangles used for the autorefinement of triangles typedef CGAL::Exact_predicates_exact_constructions_kernel EK; - std::vector< std::array > triangles(tiid+1); + std::vector< std::array > triangles(tiid+1); // TODO get rid of triangles and use all_triangle_data + // vector of data for refining triangles + std::vector all_triangle_data(triangles.size()); Cartesian_converter to_exact; for(Input_TID f : intersected_faces) { - triangles[tri_inter_ids[f]]= CGAL::make_array( + std::size_t tid=tri_inter_ids[f]; + triangles[tid]= CGAL::make_array( to_exact( get(pm, soup_points[soup_triangles[f][0]]) ), to_exact( get(pm, soup_points[soup_triangles[f][1]]) ), to_exact( get(pm, soup_points[soup_triangles[f][2]]) ) ); + all_triangle_data[tid].points.resize(3); + all_triangle_data[tid].points[0]=triangles[tri_inter_ids[f]][0]; + all_triangle_data[tid].points[1]=triangles[tri_inter_ids[f]][1]; + all_triangle_data[tid].points[2]=triangles[tri_inter_ids[f]][2]; } - std::vector< std::vector > > all_segments(triangles.size()); - std::vector< std::vector > all_points(triangles.size()); - std::vector< std::vector > all_in_triangle_ids(triangles.size()); - CGAL_PMP_AUTOREFINE_VERBOSE("compute intersections"); #ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS Real_timer t; @@ -928,7 +985,7 @@ void autorefine_triangle_soup(PointRange& soup_points, const std::array& t1 = triangles[i1]; const std::array& t2 = triangles[i2]; - std::vector inter_pts; + std::vector> inter_pts; bool triangles_are_coplanar = autorefine_impl::collect_intersections(t1, t2, inter_pts); CGAL_assertion( @@ -937,27 +994,73 @@ void autorefine_triangle_soup(PointRange& soup_points, if (!inter_pts.empty()) { + auto get_ipt_id1 = [](const std::tuple& tpl) + { + if (get<1>(tpl)<0) return -get<1>(tpl)-1; + }; + auto get_ipt_id2 = [](const std::tuple& tpl) + { + if (get<2>(tpl)<0) return -get<2>(tpl)-1; + }; + + std::size_t nbi = inter_pts.size(); switch(nbi) { + #warning TODO use inter pt info case 1: - all_points[i1].push_back(inter_pts[0]); - all_points[i2].push_back(inter_pts[0]); + if (get<1>(inter_pts[0])>=0) + all_triangle_data[i1].points.push_back(get<0>(inter_pts[0])); + if (get<2>(inter_pts[0])>=0) + all_triangle_data[i2].points.push_back(get<0>(inter_pts[0])); break; case 2: - all_segments[i1].push_back(CGAL::make_array(inter_pts[0], inter_pts[1])); - all_segments[i2].push_back(CGAL::make_array(inter_pts[0], inter_pts[1])); - all_in_triangle_ids[i1].push_back(i2); - all_in_triangle_ids[i2].push_back(i1); + { + std::size_t src_id=get<1>(inter_pts[0])<0?(-get<1>(inter_pts[0])-1):all_triangle_data[i1].points.size(); + if (get<1>(inter_pts[0])>=0) + all_triangle_data[i1].points.push_back(get<0>(inter_pts[0])); + std::size_t tgt_id=get<1>(inter_pts[1])<0?(-get<1>(inter_pts[1])-1):all_triangle_data[i1].points.size(); + if (get<1>(inter_pts[1])>=0) + all_triangle_data[i1].points.push_back(get<0>(inter_pts[1])); + all_triangle_data[i1].segments.emplace_back(src_id, tgt_id); + all_triangle_data[i1].segment_input_triangle_ids.push_back(i2); + // + src_id=get<2>(inter_pts[0])<0?(-get<2>(inter_pts[0])-1):all_triangle_data[i2].points.size(); + if (get<2>(inter_pts[0])>=0) + all_triangle_data[i2].points.push_back(get<0>(inter_pts[0])); + tgt_id=get<2>(inter_pts[1])<0?(-get<2>(inter_pts[1])-1):all_triangle_data[i2].points.size(); + if (get<2>(inter_pts[1])>=0) + all_triangle_data[i2].points.push_back(get<0>(inter_pts[1])); + all_triangle_data[i2].segments.emplace_back(src_id, tgt_id); + all_triangle_data[i2].segment_input_triangle_ids.push_back(i1); + } break; default: + { + std::vector ipt_ids1(nbi+1), ipt_ids2(nbi+1); + all_triangle_data[i1].segment_input_triangle_ids.insert( + all_triangle_data[i1].segment_input_triangle_ids.end(), nbi, i2); + all_triangle_data[i2].segment_input_triangle_ids.insert( + all_triangle_data[i2].segment_input_triangle_ids.end(), nbi, i1); + for (std::size_t i=0;i(inter_pts[i])<0?(-get<1>(inter_pts[i])-1):all_triangle_data[i1].points.size(); + std::size_t id2=get<2>(inter_pts[i])<0?(-get<2>(inter_pts[i])-1):all_triangle_data[i2].points.size(); + if (get<1>(inter_pts[i])>=0) all_triangle_data[i1].points.push_back(get<0>(inter_pts[i])); + if (get<2>(inter_pts[i])>=0) all_triangle_data[i2].points.push_back(get<0>(inter_pts[i])); + ipt_ids1[i]=id1; + ipt_ids2[i]=id2; } + ipt_ids1.back()=ipt_ids1.front(); + ipt_ids2.back()=ipt_ids2.front(); + + for (std::size_t i=0;i>> all_segments_ids(all_segments.size()); - auto deduplicate_inserted_segments = [&](std::size_t ti) { - if (!all_segments[ti].empty()) + if (!all_triangle_data[ti].segments.empty()) { - std::map point_id_map; + std::vector& points=all_triangle_data[ti].points; + std::vector>& segments=all_triangle_data[ti].segments; + std::vector indices(points.size()-3); + std::iota(indices.begin(), indices.end(),3); - auto get_point_id = [&](const EK::Point_3& pt) - { - auto insert_res = point_id_map.insert(std::make_pair(pt, all_points[ti].size())); - if (insert_res.second) - all_points[ti].push_back(pt); - return insert_res.first->second; - }; + std::sort(indices.begin(), indices.end(), [&points](std::size_t i, std::size_t j) + { return points[i] id_map(points.size()); + id_map[0]=0; + id_map[1]=1; + id_map[2]=2; + std::vector unique_ids; + unique_ids.reserve(indices.size()); - if (!all_points[ti].empty()) + //make points unique + create mapping between indices + for (std::size_t i=0; i tmp; - tmp.swap(all_points[ti]); - for (const EPoint_3& pt : tmp) - get_point_id(pt); + std::size_t new_id=unique_ids.size()+3; + unique_ids.push_back(indices[i]); + id_map[indices[i]]=new_id; + while(i+1!=indices.size() && points[indices[i]]==points[indices[i+1]]) + { + id_map[indices[++i]]=new_id; + } } - std::size_t nbs = all_segments[ti].size(); - std::vector> filtered_segments; + //~ if (unique_ids.size()==indices.size()) + //~ // TODO: do we want to keep points sorted? if yes always swap twice + //~ return; // no duplicates + + // now make points unique + using EPoint_3 = EK::Point_3; // workaround for MSVC 2022 bug + std::vector tmp; + tmp.reserve(unique_ids.size()+3); + tmp.push_back(points[0]); + tmp.push_back(points[1]); + tmp.push_back(points[2]); + for(std::size_t i : unique_ids) + tmp.push_back(points[i]); + tmp.swap(points); + + // now make segments unique + std::size_t nbs = segments.size(); + std::vector> filtered_segments; std::vector filtered_in_triangle_ids; filtered_segments.reserve(nbs); + filtered_in_triangle_ids.reserve(nbs); std::set> segset; for (std::size_t si=0; sisegments[si].first); + CGAL_assertion(id_map.size()>segments[si].second); + std::pair seg_ids = + CGAL::make_sorted_pair(id_map[segments[si].first], id_map[segments[si].second]); + if (segset.insert(seg_ids).second) { - all_segments_ids[ti].emplace_back(src_id, tgt_id); - filtered_in_triangle_ids.push_back(all_in_triangle_ids[ti][si]); + filtered_segments.push_back(seg_ids); + filtered_in_triangle_ids.push_back(all_triangle_data[ti].segment_input_triangle_ids[si]); } } - if (all_segments_ids[ti].size()!=nbs) - filtered_in_triangle_ids.swap(all_in_triangle_ids[ti]); + filtered_in_triangle_ids.swap(all_triangle_data[ti].segment_input_triangle_ids); + filtered_segments.swap(segments); + + CGAL_assertion(points.size()==std::set(points.begin(), points.end()).size()); } }; @@ -1062,11 +1189,11 @@ void autorefine_triangle_soup(PointRange& soup_points, auto refine_triangles = [&](std::size_t ti) { - if (all_segments[ti].empty() && all_points[ti].empty()) + if (all_triangle_data[ti].points.empty()) new_triangles.push_back({triangles[ti], ti}); else { - autorefine_impl::generate_subtriangles(ti, all_segments_ids[ti], all_points[ti], all_in_triangle_ids[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); } #ifdef CGAL_AUTOREF_USE_PROGRESS_DISPLAY From dbe61d2feb3621491ae09ba8abd8a95a9187e424 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 17 Jan 2024 19:33:31 +0100 Subject: [PATCH 02/15] update intersection computation code to always report the complete combinatoric of the intersection points --- .../Triangle_3_Triangle_3_intersection.h | 237 ++++++------------ 1 file changed, 76 insertions(+), 161 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h index 269c86c3e606..37fc2aa19148 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h @@ -39,11 +39,15 @@ struct Point_on_triangle { // triangle points are not stored in this class but are expected // to always be passed in the same order. For a triangle pqr, - // edge 0 is pq, edge 1 qr and edge 2 rp. Point 0 is p, 1 is q and 2 is r. + // edge 1 is pq, edge 2 qr and edge 3 rp. Point -1 is p, -2 is q and -3 is r. + // + // (0,-) vertex of t2 inside t1 + // (-,-) common point of t1 and t2 + // (-,0) vertex of t1 inside t2 + // (-,+) vertex of t1 on an edge of t2 + // (+,+) edge of t1 intersecting edge of t2 + // (+,-) vertex of t2 on an edge of t1 // - // (id, -1) point on t1 - // (-1, id) point on t2 - // (id1, id2) intersection of edges std::pair t1_t2_ids; boost::container::flat_set extra_t1; // store other ids of edges containing the point typename Kernel::FT alpha; // @@ -60,16 +64,16 @@ struct Point_on_triangle { switch(id) { - case 0: + case -1: return p; - case 1: + case -2: return q; default: return r; } } - Point_on_triangle(int i1, int i2=-1, typename Kernel::FT alpha = 0.) // TODO add global zero()? + Point_on_triangle(int i1, int i2, typename Kernel::FT alpha = 0.) // TODO add global zero()? : t1_t2_ids(i1,i2) , alpha(alpha) {} @@ -85,20 +89,20 @@ struct Point_on_triangle const typename Kernel::Point_3& r2, const Kernel& k) const { - if (t1_t2_ids.first!=-1) + if (t1_t2_ids.first!=0 && t1_t2_ids.second >=0) { - if (t1_t2_ids.second==-1) + if (t1_t2_ids.first<0) return (edge_id1==t1_t2_ids.first || (edge_id1+1)%3==t1_t2_ids.first) ? ZERO:POSITIVE; // it is a point on t1 // this is an intersection point if (t1_t2_ids.first==edge_id1) return ZERO; - if (t1_t2_ids.first==(edge_id1+1)%3) + if (t1_t2_ids.first-1==(edge_id1)%3) { if (alpha==0) return ZERO; return alpha>=0 ? POSITIVE:NEGATIVE; } - CGAL_assertion((t1_t2_ids.first+1)%3==edge_id1); + CGAL_assertion((t1_t2_ids.first)%3==edge_id1-1); if (alpha==1) return ZERO; return alpha<=1?POSITIVE:NEGATIVE; } @@ -111,9 +115,13 @@ struct Point_on_triangle } } + int id1() const { return t1_t2_ids.first; } int id2() const { return t1_t2_ids.second; } + void set_id1(int i) { t1_t2_ids.first=i; } + void set_id2(int i) { t1_t2_ids.second=i; } + // construct the intersection point from the info stored typename Kernel::Point_3 point(const typename Kernel::Point_3& p1, @@ -124,12 +132,13 @@ struct Point_on_triangle const typename Kernel::Point_3& r2, const Kernel& k) const { - if (t1_t2_ids.first==-1) + if (t1_t2_ids.second<0) return point_from_id(p2,q2,r2,t1_t2_ids.second); - if (t1_t2_ids.second==-1) + if (t1_t2_ids.first<0) return point_from_id(p1,q1,r1,t1_t2_ids.first); - return k.construct_barycenter_3_object()(point_from_id(p1,q1,r1,(t1_t2_ids.first+1)%3), alpha, point_from_id(p1,q1,r1,t1_t2_ids.first)) ; + return k.construct_barycenter_3_object()(point_from_id(p1,q1,r1,-(t1_t2_ids.first)%3), alpha, + point_from_id(p1,q1,r1,-t1_t2_ids.first)) ; } }; @@ -160,153 +169,50 @@ intersection(const Point_on_triangle& p, typename Kernel::Compute_alpha_for_coplanar_triangle_intersection_3 compute_alpha = k.compute_alpha_for_coplanar_triangle_intersection_3_object(); typedef Point_on_triangle Pot; - switch(p.id1()) + + auto common_edge_id = [](int pid, int qid) { - case -1: + if (pid==0 || qid==0) return 0; + if (pid<0) { - switch(q.id1()) + if (qid<0) { - case -1: // A: (-1, ip2) - (-1, iq2) - { - CGAL_assertion((p.id2()+1)%3 == q.id2() || (q.id2()+1)%3 == p.id2()); -// CGAL_assertion(p.extra_t1.empty() && q.extra_t1.empty()); // TMP to see if it's worth implementing special case -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 1\n"; -#endif - typename Kernel::FT alpha = compute_alpha(p1, q1, - Pot::point_from_id(p2, q2, r2, p.id2()), - Pot::point_from_id(p2, q2, r2, q.id2())); - int id2 = (p.id2()+1)%3 == q.id2() ? p.id2() : q.id2(); - return Point_on_triangle(edge_id_t1, id2, alpha); // intersection with an original edge of t2 - } - default: - if (q.id2()!=-1) // B: (-1, ip2) - (iq1, iq2) - { - if (p.id2() == q.id2() || p.id2() == (q.id2()+1)%3) - { -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 2\n"; -#endif - // points are on the same edge of t2 --> we shorten an already cut edge - typename Kernel::FT alpha = compute_alpha(p1, q1, - Pot::point_from_id(p2, q2, r2, q.id2()), - Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3)); - - return Point_on_triangle(edge_id_t1, q.id2(), alpha); - } -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 3\n"; -#endif - // point of t1: look for an edge of t1 containing both points - CGAL_assertion( p.extra_t1.count(q.id1())!=0 || p.extra_t1.count(3-q.id1()-edge_id_t1)!=0 ); - int eid1 = p.extra_t1.count(q.id1())!=0 ? q.id1() : 3-q.id1()-edge_id_t1; - return Point_on_triangle((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1 - } - // C: (-1, ip2) - (iq1, -1) - //vertex of t1, special case t1 edge passed thru a vertex of t2 - CGAL_assertion(edge_id_t1 == 2); -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 4\n"; -#endif - CGAL_assertion(q.id1()==1); - CGAL_assertion(!p.extra_t1.empty()); - return Point_on_triangle(p.extra_t1.count(0)==1?0:2,-1); + return (-pid)%3==-qid-1 ? -pid:-qid; + } + else + { + return (-pid==qid || (qid)%3==-pid-1) ? qid : 0; } } - default: + else { - switch(p.id2()) + if (qid<0) { - case -1: - { - switch(q.id1()) - { - case -1: // G: (ip1, -1) - (-1, iq2) - //vertex of t1, special case t1 edge passed thru a vertex of t2 -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 5\n"; -#endif - CGAL_assertion(edge_id_t1 == 2); - CGAL_assertion(p.id1()==1); - CGAL_assertion(!q.extra_t1.empty()); - return Point_on_triangle(q.extra_t1.count(0)==1?0:2,-1); - default: - { -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 6\n"; -#endif - CGAL_assertion(q.id2()!=-1); // I: (ip1, -1) - (iq2, -1) - //H: (ip1,-1), (iq1, iq2) - CGAL_assertion(edge_id_t1==2); - // p and q are on the same edge of t1 - CGAL_assertion(p.id1()==q.id1() || p.id1()==(q.id1()+1)%3); - return Point_on_triangle((q.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1); - } - } - } - default: - { - switch(q.id1()) - { - case -1: // D: (ip1, ip2) - (-1, iq2) - { - if (q.id2() == p.id2() || q.id2() == (p.id2()+1)%3) - { -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 7\n"; -#endif - // points are on the same edge of t2 --> we shorten an already cut edge - typename Kernel::FT alpha = compute_alpha(p1, q1, - Pot::point_from_id(p2, q2, r2, p.id2()), - Pot::point_from_id(p2, q2, r2, (p.id2()+1)%3)); - - return Point_on_triangle(edge_id_t1, p.id2(), alpha); - } -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 8\n"; -#endif - // point of t1 - //std::cout << "q.extra_t1: "; for(int qet1 : q.extra_t1) std::cout << " " << qet1; std::cout << "\n"; - CGAL_assertion( q.extra_t1.count(p.id1())!=0 || q.extra_t1.count(3-p.id1()-edge_id_t1)!=0 ); - int eid1 = q.extra_t1.count(p.id1())!=0 ? p.id1() : 3-p.id1()-edge_id_t1; - return Point_on_triangle((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1 - } - default: - { - switch(q.id2()) - { - case -1: // F: (ip1, ip2) - (iq1, -1) - { - // p and q are on the same edge of t1 - CGAL_assertion(q.id1()==p.id1() || q.id1()==(p.id1()+1)%3); -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 9\n"; -#endif - return Point_on_triangle((p.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1); - } - default: // E: (ip1, ip2) - (iq1, iq2) - { - if (p.id2()==q.id2()) - { -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 10\n"; -#endif - typename Kernel::FT alpha = compute_alpha(p1, q1, - Pot::point_from_id(p2, q2, r2, q.id2()), - Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3)); - return Point_on_triangle(edge_id_t1, q.id2(), alpha); - } - // we are intersecting an edge of t1 - CGAL_assertion(p.id1()==q.id1() || edge_id_t1==2); - int eid1 = p.id1()==q.id1() ? p.id1() : 1; - return Point_on_triangle((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1 - } - } - } - } - } + return (-qid==pid || (pid)%3==-qid-1) ? pid : 0; + } + else + { + return pid==qid ? pid : 0; } } + }; + + int common_eid1 = common_edge_id(p.id1(), q.id1()); + int common_eid2 = common_edge_id(p.id2(), q.id2()); + + if (common_eid1!=0) + { + CGAL_assertion( (common_eid1)%3==edge_id_t1-1 || (edge_id_t1)%3==common_eid1-1 ); + int pid1 = (common_eid1)%3==edge_id_t1-1? -edge_id_t1:-common_eid1; + return Point_on_triangle(pid1, common_eid2); + } + else + { + CGAL_assertion(common_eid2!=0); + typename Kernel::FT alpha = compute_alpha(p1, q1, + Pot::point_from_id(p2, q2, r2, -common_eid2), + Pot::point_from_id(p2, q2, r2, -(common_eid2)%3-1)); + return Point_on_triangle(edge_id_t1, common_eid2, alpha); } } @@ -336,8 +242,17 @@ void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p1, for (Point_on_triangle& pot : inter_pts) { orientations[ &pot ]=pot.orientation(p1,q1,r1,edge_id,p2,q2,r2,k); - if (pot.id1()==-1 && orientations[ &pot ]==COLLINEAR) - pot.extra_t1.insert(edge_id); + if (orientations[ &pot ]==COLLINEAR) + { + if (pot.id1()==0) + pot.set_id1(edge_id); + else + { + CGAL_assertion(pot.id1()>0); + CGAL_assertion((edge_id)%3==pot.id1()-1 || (pot.id1())%3==edge_id-1); + pot.set_id1( (edge_id)%3==pot.id1()-1 ? -pot.id1(): -edge_id ); + } + } } #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION @@ -411,30 +326,30 @@ intersection_coplanar_triangles(const typename K::Triangle_3& t1, r2 = t2.vertex(2); std::list> inter_pts; - inter_pts.push_back(Point_on_triangle(-1,0)); - inter_pts.push_back(Point_on_triangle(-1,1)); - inter_pts.push_back(Point_on_triangle(-1,2)); + inter_pts.push_back(Point_on_triangle(0,-1)); + inter_pts.push_back(Point_on_triangle(0,-2)); + inter_pts.push_back(Point_on_triangle(0,-3)); #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION auto print_points = [&]() { - for(auto p : inter_pts) std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) "; std::cout <<"\n"; + for(auto p : inter_pts) {std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) ";} std::cout <<"\n"; }; std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif //intersect t2 with the three half planes which intersection defines t1 - intersection_coplanar_triangles_cutoff(p1,q1,r1,0,p2,q2,r2,k,inter_pts); //line pq + intersection_coplanar_triangles_cutoff(p1,q1,r1,1,p2,q2,r2,k,inter_pts); //line pq #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif - intersection_coplanar_triangles_cutoff(q1,r1,p1,1,p2,q2,r2,k,inter_pts); //line qr + intersection_coplanar_triangles_cutoff(q1,r1,p1,2,p2,q2,r2,k,inter_pts); //line qr #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif - intersection_coplanar_triangles_cutoff(r1,p1,q1,2,p2,q2,r2,k,inter_pts); //line rp + intersection_coplanar_triangles_cutoff(r1,p1,q1,3,p2,q2,r2,k,inter_pts); //line rp #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); From b9232b36772dd7ebb1a3cad8773474b02ccb6d20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 18 Jan 2024 10:56:36 +0100 Subject: [PATCH 03/15] fix bad indices --- .../internal/Triangle_3_Triangle_3_intersection.h | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h index 37fc2aa19148..471904cd97ab 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h @@ -137,7 +137,7 @@ struct Point_on_triangle if (t1_t2_ids.first<0) return point_from_id(p1,q1,r1,t1_t2_ids.first); - return k.construct_barycenter_3_object()(point_from_id(p1,q1,r1,-(t1_t2_ids.first)%3), alpha, + return k.construct_barycenter_3_object()(point_from_id(p1,q1,r1,-(t1_t2_ids.first)%3-1), alpha, point_from_id(p1,q1,r1,-t1_t2_ids.first)) ; } }; @@ -164,7 +164,7 @@ intersection(const Point_on_triangle& p, #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION std::cout << " calling intersection: "; std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) -"; - std::cout << " (" << q.id1() << "," << q.id2() << ",[" << q.alpha << "]) || e" << edge_id_t1; + std::cout << " (" << q.id1() << "," << q.id2() << ",[" << q.alpha << "]) || e" << edge_id_t1 << "\n"; #endif typename Kernel::Compute_alpha_for_coplanar_triangle_intersection_3 compute_alpha = k.compute_alpha_for_coplanar_triangle_intersection_3_object(); @@ -333,25 +333,22 @@ intersection_coplanar_triangles(const typename K::Triangle_3& t1, #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION auto print_points = [&]() { + std::cout << " ipts size: " << inter_pts.size() << "\n"; for(auto p : inter_pts) {std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) ";} std::cout <<"\n"; }; - std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif //intersect t2 with the three half planes which intersection defines t1 intersection_coplanar_triangles_cutoff(p1,q1,r1,1,p2,q2,r2,k,inter_pts); //line pq #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif intersection_coplanar_triangles_cutoff(q1,r1,p1,2,p2,q2,r2,k,inter_pts); //line qr #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif intersection_coplanar_triangles_cutoff(r1,p1,q1,3,p2,q2,r2,k,inter_pts); //line rp #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif From 61461b14b3981d3f2582808ddc2ef9f3c06447a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 18 Jan 2024 10:57:17 +0100 Subject: [PATCH 04/15] more debug and fix indices --- .../Polygon_mesh_processing/autorefinement.h | 95 ++++++++++++++----- 1 file changed, 73 insertions(+), 22 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index 1cb0e3aaa4d8..50e7c5b361e5 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -250,7 +250,6 @@ do_coplanar_segments_intersect(std::size_t pi, std::size_t qi, // imported from Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h -#warning TODO fill the indices correctly template void coplanar_intersections(const std::array& t1, const std::array& t2, @@ -260,18 +259,37 @@ void coplanar_intersections(const std::array& t1, const typename K::Point_3& p2 = t2[0], q2 = t2[1], r2 = t2[2]; std::list> l_inter_pts; - l_inter_pts.push_back(Intersections::internal::Point_on_triangle(-1,0)); - l_inter_pts.push_back(Intersections::internal::Point_on_triangle(-1,1)); - l_inter_pts.push_back(Intersections::internal::Point_on_triangle(-1,2)); + l_inter_pts.push_back(Intersections::internal::Point_on_triangle(0,-1)); + l_inter_pts.push_back(Intersections::internal::Point_on_triangle(0,-2)); + l_inter_pts.push_back(Intersections::internal::Point_on_triangle(0,-3)); +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + auto print_points = [&]() + { + std::cout << " ipts size: " << l_inter_pts.size() << "\n"; + for(auto p : l_inter_pts) {std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) ";} std::cout <<"\n"; + }; +#endif //intersect t2 with the three half planes which intersection defines t1 K k; - Intersections::internal::intersection_coplanar_triangles_cutoff(p1,q1,r1,0,p2,q2,r2,k,l_inter_pts); //line p1q1 - Intersections::internal::intersection_coplanar_triangles_cutoff(q1,r1,p1,1,p2,q2,r2,k,l_inter_pts); //line q1r1 - Intersections::internal::intersection_coplanar_triangles_cutoff(r1,p1,q1,2,p2,q2,r2,k,l_inter_pts); //line r1p1 +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + print_points(); +#endif + Intersections::internal::intersection_coplanar_triangles_cutoff(p1,q1,r1,1,p2,q2,r2,k,l_inter_pts); //line p1q1 +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + print_points(); +#endif + Intersections::internal::intersection_coplanar_triangles_cutoff(q1,r1,p1,2,p2,q2,r2,k,l_inter_pts); //line q1r1 +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + print_points(); +#endif + Intersections::internal::intersection_coplanar_triangles_cutoff(r1,p1,q1,3,p2,q2,r2,k,l_inter_pts); //line r1p1 +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + print_points(); +#endif for (const Intersections::internal::Point_on_triangle& pot : l_inter_pts) - inter_pts.emplace_back(pot.point(p1,q1,r1,p2,q2,r2,k), 0, 0); + inter_pts.emplace_back(pot.point(p1,q1,r1,p2,q2,r2,k), pot.id1(), pot.id2()); } // imported from Polygon_mesh_processing/internal/Corefinement/intersect_triangle_and_segment_3.h @@ -891,8 +909,7 @@ void autorefine_triangle_soup(PointRange& soup_points, Visitor visitor(choose_parameter(get_parameter(np, internal_np::visitor))); - //~ constexpr bool parallel_execution = std::is_same_v; - constexpr bool parallel_execution = false; + constexpr bool parallel_execution = std::is_same_v; #ifndef CGAL_LINKED_WITH_TBB static_assert (!parallel_execution, @@ -994,43 +1011,63 @@ void autorefine_triangle_soup(PointRange& soup_points, if (!inter_pts.empty()) { - auto get_ipt_id1 = [](const std::tuple& tpl) - { - if (get<1>(tpl)<0) return -get<1>(tpl)-1; - }; - auto get_ipt_id2 = [](const std::tuple& tpl) - { - if (get<2>(tpl)<0) return -get<2>(tpl)-1; - }; - - std::size_t nbi = inter_pts.size(); switch(nbi) { #warning TODO use inter pt info case 1: if (get<1>(inter_pts[0])>=0) + { all_triangle_data[i1].points.push_back(get<0>(inter_pts[0])); + CGAL_assertion(get<0>(inter_pts[0])!=t1[0]); + CGAL_assertion(get<0>(inter_pts[0])!=t1[1]); + CGAL_assertion(get<0>(inter_pts[0])!=t1[2]); + } if (get<2>(inter_pts[0])>=0) + { all_triangle_data[i2].points.push_back(get<0>(inter_pts[0])); + CGAL_assertion(get<0>(inter_pts[0])!=t2[0]); + CGAL_assertion(get<0>(inter_pts[0])!=t2[1]); + CGAL_assertion(get<0>(inter_pts[0])!=t2[2]); + } break; case 2: { std::size_t src_id=get<1>(inter_pts[0])<0?(-get<1>(inter_pts[0])-1):all_triangle_data[i1].points.size(); if (get<1>(inter_pts[0])>=0) + { all_triangle_data[i1].points.push_back(get<0>(inter_pts[0])); + CGAL_assertion(get<0>(inter_pts[0])!=t1[0]); + CGAL_assertion(get<0>(inter_pts[0])!=t1[1]); + CGAL_assertion(get<0>(inter_pts[0])!=t1[2]); + } std::size_t tgt_id=get<1>(inter_pts[1])<0?(-get<1>(inter_pts[1])-1):all_triangle_data[i1].points.size(); if (get<1>(inter_pts[1])>=0) + { all_triangle_data[i1].points.push_back(get<0>(inter_pts[1])); + CGAL_assertion(get<0>(inter_pts[1])!=t1[0]); + CGAL_assertion(get<0>(inter_pts[1])!=t1[1]); + CGAL_assertion(get<0>(inter_pts[1])!=t1[2]); + } all_triangle_data[i1].segments.emplace_back(src_id, tgt_id); all_triangle_data[i1].segment_input_triangle_ids.push_back(i2); // src_id=get<2>(inter_pts[0])<0?(-get<2>(inter_pts[0])-1):all_triangle_data[i2].points.size(); if (get<2>(inter_pts[0])>=0) + { all_triangle_data[i2].points.push_back(get<0>(inter_pts[0])); + CGAL_assertion(get<0>(inter_pts[0])!=t2[0]); + CGAL_assertion(get<0>(inter_pts[0])!=t2[1]); + CGAL_assertion(get<0>(inter_pts[0])!=t2[2]); + } tgt_id=get<2>(inter_pts[1])<0?(-get<2>(inter_pts[1])-1):all_triangle_data[i2].points.size(); if (get<2>(inter_pts[1])>=0) + { all_triangle_data[i2].points.push_back(get<0>(inter_pts[1])); + CGAL_assertion(get<0>(inter_pts[1])!=t2[0]); + CGAL_assertion(get<0>(inter_pts[1])!=t2[1]); + CGAL_assertion(get<0>(inter_pts[1])!=t2[2]); + } all_triangle_data[i2].segments.emplace_back(src_id, tgt_id); all_triangle_data[i2].segment_input_triangle_ids.push_back(i1); } @@ -1047,8 +1084,21 @@ void autorefine_triangle_soup(PointRange& soup_points, { std::size_t id1=get<1>(inter_pts[i])<0?(-get<1>(inter_pts[i])-1):all_triangle_data[i1].points.size(); std::size_t id2=get<2>(inter_pts[i])<0?(-get<2>(inter_pts[i])-1):all_triangle_data[i2].points.size(); - if (get<1>(inter_pts[i])>=0) all_triangle_data[i1].points.push_back(get<0>(inter_pts[i])); - if (get<2>(inter_pts[i])>=0) all_triangle_data[i2].points.push_back(get<0>(inter_pts[i])); + if (get<1>(inter_pts[i])>=0) + { + all_triangle_data[i1].points.push_back(get<0>(inter_pts[i])); + CGAL_assertion(get<0>(inter_pts[i])!=t1[0]); + CGAL_assertion(get<0>(inter_pts[i])!=t1[1]); + CGAL_assertion(get<0>(inter_pts[i])!=t1[2]); + } + + if (get<2>(inter_pts[i])>=0) + { + all_triangle_data[i2].points.push_back(get<0>(inter_pts[i])); + CGAL_assertion(get<0>(inter_pts[i])!=t2[0]); + CGAL_assertion(get<0>(inter_pts[i])!=t2[1]); + CGAL_assertion(get<0>(inter_pts[i])!=t2[2]); + } ipt_ids1[i]=id1; ipt_ids2[i]=id2; } @@ -1147,6 +1197,7 @@ void autorefine_triangle_soup(PointRange& soup_points, filtered_in_triangle_ids.swap(all_triangle_data[ti].segment_input_triangle_ids); filtered_segments.swap(segments); + CGAL_assertion(points.size()-3==std::set(std::next(points.begin(),3), points.end()).size()); CGAL_assertion(points.size()==std::set(points.begin(), points.end()).size()); } }; From 81d293891007be8eea4fdaf087066f7e845d2467 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 18 Jan 2024 13:26:36 +0100 Subject: [PATCH 05/15] collect on-edge information and do not collect segments on the same edge --- .../Polygon_mesh_processing/autorefinement.h | 169 ++++++++++-------- 1 file changed, 93 insertions(+), 76 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index 50e7c5b361e5..400e8641de8a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -485,9 +485,58 @@ bool collect_intersections(const std::array& t1, //TODO: rename struct struct Triangle_data { - std::vector points; + using Point_3 = Exact_predicates_exact_constructions_kernel::Point_3; + std::vector points; + std::vector point_locations; std::vector> segments; std::vector segment_input_triangle_ids; + template + std::size_t add_point(const std::tuple& tpl) + { + if (get(tpl) < 0) + return -get(tpl)-1; + points.push_back(get<0>(tpl)); + point_locations.push_back(get(tpl)); + return points.size()-1; + } + +#ifndef NDEBUG + bool on_same_edge(std::size_t i, std::size_t j) +#else + bool on_same_edge_debug(std::size_t i, std::size_t j) +#endif + { + if (point_locations[i]==0 || point_locations[j]==0) return false; + if (point_locations[i]>0) + { + if (point_locations[j]>0) + return point_locations[i]==point_locations[j]; + return -point_locations[j]==point_locations[i] || point_locations[i]%3+1==-point_locations[j]; + } + if (point_locations[j]<0) + return true; + return -point_locations[i]==point_locations[j] || point_locations[j]%3+1==-point_locations[i]; + } +#ifdef NDEBUG + bool on_same_edge(std::size_t i, std::size_t j) + { + bool on = on_same_edge_debug(i,j); + if (!on) return false; + int eid = point_locations[i]>0 ? point_locations[i] : point_locations[j]; + if (eid<0) return true; + if (!(collinear(points[eid-1], points[(eid)%3], points[i]))) + { + std::cout << point_locations[i] << " " << point_locations[j] << " " << eid << "\n"; + } + if (!(collinear(points[eid-1], points[(eid)%3], points[j]))) + { + std::cout << point_locations[i] << " " << point_locations[j] << " " << eid << "\n"; + } + CGAL_assertion(collinear(points[eid-1], points[(eid)%3], points[i])); + CGAL_assertion(collinear(points[eid-1], points[(eid)%3], points[j])); + return true; + } +#endif }; template (inter_pts[0])>=0) - { - all_triangle_data[i1].points.push_back(get<0>(inter_pts[0])); - CGAL_assertion(get<0>(inter_pts[0])!=t1[0]); - CGAL_assertion(get<0>(inter_pts[0])!=t1[1]); - CGAL_assertion(get<0>(inter_pts[0])!=t1[2]); - } - if (get<2>(inter_pts[0])>=0) - { - all_triangle_data[i2].points.push_back(get<0>(inter_pts[0])); - CGAL_assertion(get<0>(inter_pts[0])!=t2[0]); - CGAL_assertion(get<0>(inter_pts[0])!=t2[1]); - CGAL_assertion(get<0>(inter_pts[0])!=t2[2]); - } + all_triangle_data[i1].add_point<1>(inter_pts[0]); + all_triangle_data[i2].add_point<2>(inter_pts[0]); break; case 2: { - std::size_t src_id=get<1>(inter_pts[0])<0?(-get<1>(inter_pts[0])-1):all_triangle_data[i1].points.size(); - if (get<1>(inter_pts[0])>=0) - { - all_triangle_data[i1].points.push_back(get<0>(inter_pts[0])); - CGAL_assertion(get<0>(inter_pts[0])!=t1[0]); - CGAL_assertion(get<0>(inter_pts[0])!=t1[1]); - CGAL_assertion(get<0>(inter_pts[0])!=t1[2]); - } - std::size_t tgt_id=get<1>(inter_pts[1])<0?(-get<1>(inter_pts[1])-1):all_triangle_data[i1].points.size(); - if (get<1>(inter_pts[1])>=0) + std::size_t src_id=all_triangle_data[i1].add_point<1>(inter_pts[0]), + tgt_id=all_triangle_data[i1].add_point<1>(inter_pts[1]); + if (!all_triangle_data[i1].on_same_edge(src_id, tgt_id)) { - all_triangle_data[i1].points.push_back(get<0>(inter_pts[1])); - CGAL_assertion(get<0>(inter_pts[1])!=t1[0]); - CGAL_assertion(get<0>(inter_pts[1])!=t1[1]); - CGAL_assertion(get<0>(inter_pts[1])!=t1[2]); + all_triangle_data[i1].segments.emplace_back(src_id, tgt_id); + all_triangle_data[i1].segment_input_triangle_ids.push_back(i2); } - all_triangle_data[i1].segments.emplace_back(src_id, tgt_id); - all_triangle_data[i1].segment_input_triangle_ids.push_back(i2); // - src_id=get<2>(inter_pts[0])<0?(-get<2>(inter_pts[0])-1):all_triangle_data[i2].points.size(); - if (get<2>(inter_pts[0])>=0) + src_id=all_triangle_data[i2].add_point<2>(inter_pts[0]); + tgt_id=all_triangle_data[i2].add_point<2>(inter_pts[1]); + if (!all_triangle_data[i2].on_same_edge(src_id, tgt_id)) { - all_triangle_data[i2].points.push_back(get<0>(inter_pts[0])); - CGAL_assertion(get<0>(inter_pts[0])!=t2[0]); - CGAL_assertion(get<0>(inter_pts[0])!=t2[1]); - CGAL_assertion(get<0>(inter_pts[0])!=t2[2]); + all_triangle_data[i2].segments.emplace_back(src_id, tgt_id); + all_triangle_data[i2].segment_input_triangle_ids.push_back(i1); } - tgt_id=get<2>(inter_pts[1])<0?(-get<2>(inter_pts[1])-1):all_triangle_data[i2].points.size(); - if (get<2>(inter_pts[1])>=0) - { - all_triangle_data[i2].points.push_back(get<0>(inter_pts[1])); - CGAL_assertion(get<0>(inter_pts[1])!=t2[0]); - CGAL_assertion(get<0>(inter_pts[1])!=t2[1]); - CGAL_assertion(get<0>(inter_pts[1])!=t2[2]); - } - all_triangle_data[i2].segments.emplace_back(src_id, tgt_id); - all_triangle_data[i2].segment_input_triangle_ids.push_back(i1); } break; default: { std::vector ipt_ids1(nbi+1), ipt_ids2(nbi+1); - all_triangle_data[i1].segment_input_triangle_ids.insert( - all_triangle_data[i1].segment_input_triangle_ids.end(), nbi, i2); - all_triangle_data[i2].segment_input_triangle_ids.insert( - all_triangle_data[i2].segment_input_triangle_ids.end(), nbi, i1); for (std::size_t i=0;i(inter_pts[i])<0?(-get<1>(inter_pts[i])-1):all_triangle_data[i1].points.size(); - std::size_t id2=get<2>(inter_pts[i])<0?(-get<2>(inter_pts[i])-1):all_triangle_data[i2].points.size(); - if (get<1>(inter_pts[i])>=0) - { - all_triangle_data[i1].points.push_back(get<0>(inter_pts[i])); - CGAL_assertion(get<0>(inter_pts[i])!=t1[0]); - CGAL_assertion(get<0>(inter_pts[i])!=t1[1]); - CGAL_assertion(get<0>(inter_pts[i])!=t1[2]); - } - - if (get<2>(inter_pts[i])>=0) - { - all_triangle_data[i2].points.push_back(get<0>(inter_pts[i])); - CGAL_assertion(get<0>(inter_pts[i])!=t2[0]); - CGAL_assertion(get<0>(inter_pts[i])!=t2[1]); - CGAL_assertion(get<0>(inter_pts[i])!=t2[2]); - } - ipt_ids1[i]=id1; - ipt_ids2[i]=id2; + ipt_ids1[i]=all_triangle_data[i1].add_point<1>(inter_pts[i]); + ipt_ids2[i]=all_triangle_data[i2].add_point<2>(inter_pts[i]); } ipt_ids1.back()=ipt_ids1.front(); ipt_ids2.back()=ipt_ids2.front(); for (std::size_t i=0;i tmp; + std::vector tmp_locations; tmp.reserve(unique_ids.size()+3); + tmp_locations.reserve(unique_ids.size()+3); tmp.push_back(points[0]); tmp.push_back(points[1]); tmp.push_back(points[2]); + tmp_locations.push_back(-1); + tmp_locations.push_back(-2); + tmp_locations.push_back(-3); for(std::size_t i : unique_ids) + { tmp.push_back(points[i]); + tmp_locations.push_back(all_triangle_data[ti].point_locations[i]); + } tmp.swap(points); + tmp_locations.swap(all_triangle_data[ti].point_locations); + CGAL_assertion(points.size() == all_triangle_data[ti].point_locations.size()); // now make segments unique std::size_t nbs = segments.size(); @@ -1240,7 +1257,7 @@ void autorefine_triangle_soup(PointRange& soup_points, auto refine_triangles = [&](std::size_t ti) { - if (all_triangle_data[ti].points.empty()) + if (all_triangle_data[ti].points.size()==3) new_triangles.push_back({triangles[ti], ti}); else { From 0095748cd7a5ab55eb82767b82f4e68e9e051ec5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 18 Jan 2024 14:42:47 +0100 Subject: [PATCH 06/15] insert points directly on edges --- .../Polygon_mesh_processing/autorefinement.h | 60 +++++++++++++++++-- 1 file changed, 55 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index 400e8641de8a..dcd6a0664abc 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -582,9 +582,8 @@ void generate_subtriangles(std::size_t ti, #define CGAL_AUTOREF_COUNTER_INSTRUCTION(X) #endif - #warning TODO - std::vector& points=triangle_data.points; + std::vector& point_locations=triangle_data.point_locations; std::vector>& segments=triangle_data.segments; std::vector& in_triangle_ids=triangle_data.segment_input_triangle_ids; @@ -611,7 +610,10 @@ void generate_subtriangles(std::size_t ti, { auto insert_res = point_id_map.insert(std::make_pair(pt, points.size())); if (insert_res.second) + { points.push_back(pt); + point_locations.push_back(0); + } return insert_res.first->second; }; @@ -840,8 +842,57 @@ void generate_subtriangles(std::size_t ti, vhandles[2]->set_point(t[2]); // insert points and fill vhandles - std::vector indices(triangle_data.points.size()-3); - std::iota(indices.begin(), indices.end(), 3); + + //start by points on edges + std::array, 3> indices_on_edges; + std::vector indices; + for (std::size_t i=3; i0 && point_locations[i]<4); + indices_on_edges[point_locations[i]-1].push_back(i); + } + } + + //sort points on edges and insert them + typename CDT::Vertex_handle vinf=cdt.infinite_vertex(); + for (int i=0; i<3; ++i) + { + if (indices_on_edges[i].empty()) continue; + int src_id=i, tgt_id=(i+1)%3; + //look for a sort axis + int coord = 0; + if (points[src_id].x()==points[tgt_id].x()) + { + coord=1; + if (points[src_id].y()==points[tgt_id].y()) + coord=2; + } + if (points[src_id][coord]>points[tgt_id][coord]) + std::swap(src_id, tgt_id); + + std::sort(indices_on_edges[i].begin(), indices_on_edges[i].end(), + [&](std::size_t id1, std::size_t id2) + { + return points[id1][coord]index(vinf)); + prev_id=id; + } + } + + // then points in the interior typedef typename Pointer_property_map::type Pmap; typedef Spatial_sort_traits_adapter_2 Search_traits; spatial_sort(indices.begin(), indices.end(), @@ -1067,7 +1118,6 @@ void autorefine_triangle_soup(PointRange& soup_points, std::size_t nbi = inter_pts.size(); switch(nbi) { - #warning TODO use inter pt info case 1: all_triangle_data[i1].add_point<1>(inter_pts[0]); all_triangle_data[i2].add_point<2>(inter_pts[0]); From 272a242f1b049af73bc2530e660fba4a14093142 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 19 Jan 2024 10:06:03 +0100 Subject: [PATCH 07/15] restore delaunay in CDT + autoref and deduplicate identical points in no segments case --- .../Polygon_mesh_processing/autorefinement.h | 26 +++++++++++++++++-- .../Constrained_Delaunay_triangulation_2.h | 19 ++++++++++++++ 2 files changed, 43 insertions(+), 2 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index dcd6a0664abc..8de91af1da85 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -843,13 +843,22 @@ void generate_subtriangles(std::size_t ti, // insert points and fill vhandles +#if 0 + std::vector indices(triangle_data.points.size()-3); + std::iota(indices.begin(), indices.end(), 3); +#else //start by points on edges std::array, 3> indices_on_edges; std::vector indices; for (std::size_t i=3; i0 && point_locations[i]<4); @@ -862,7 +871,7 @@ void generate_subtriangles(std::size_t ti, for (int i=0; i<3; ++i) { if (indices_on_edges[i].empty()) continue; - int src_id=i, tgt_id=(i+1)%3; + std::size_t src_id=i, tgt_id=(i+1)%3; //look for a sort axis int coord = 0; if (points[src_id].x()==points[tgt_id].x()) @@ -879,19 +888,32 @@ void generate_subtriangles(std::size_t ti, { return points[id1][coord]index(vinf)); + cdt.restore_Delaunay(vhandles[id]); // TODO maybe not each time but one global? + CGAL_assertion(cdt.is_valid()); prev_id=id; } } - +#endif // then points in the interior typedef typename Pointer_property_map::type Pmap; typedef Spatial_sort_traits_adapter_2 Search_traits; diff --git a/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h b/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h index 52c14ea54b78..714e9295b61f 100644 --- a/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h +++ b/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h @@ -342,6 +342,25 @@ class Constrained_Delaunay_triangulation_2 return number_of_vertices() - n; } + // function usable only if no constraint has been inserted + void restore_Delaunay(Vertex_handle v) + { + if(this->dimension() <= 1) return; + + Face_handle f=v->face(); + Face_handle next; + int i; + Face_handle start(f); + do { + i = f->index(v); + next = f->neighbor(ccw(i)); // turn ccw around v + propagating_flip(f,i); + f = next; + } while(next != start); + return; + } + + #ifndef CGAL_TRIANGULATION_2_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO private: //top stands for tuple-or-pair From f520602a0dcbcefe0a28a346f11fd84b15249f34 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 19 Jan 2024 10:05:30 +0000 Subject: [PATCH 08/15] Comment #warning and qualify get() with std:: --- .../CGAL/Polygon_mesh_processing/autorefinement.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index 8de91af1da85..c38d783572b4 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -468,10 +468,10 @@ bool collect_intersections(const std::array& t1, } } - #warning TODO get rid of sort and unique calls +// #warning TODO get rid of sort and unique calls // because we don't handle intersection type and can have edge-edge edge-vertex duplicates - std::sort(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return get<0>(p)(q);}); - auto last = std::unique(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return get<0>(p)==get<0>(q);}); + std::sort(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return std::get<0>(p)(q);}); + auto last = std::unique(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return std::get<0>(p)==std::get<0>(q);}); inter_pts.erase(last, inter_pts.end()); #ifdef CGAL_AUTOREF_DEBUG_DEPTH @@ -493,10 +493,10 @@ struct Triangle_data template std::size_t add_point(const std::tuple& tpl) { - if (get(tpl) < 0) - return -get(tpl)-1; - points.push_back(get<0>(tpl)); - point_locations.push_back(get(tpl)); + if (std::get(tpl) < 0) + return -std::get(tpl)-1; + points.push_back(std::get<0>(tpl)); + point_locations.push_back(std::get(tpl)); return points.size()-1; } From a783412ba77b6a0873bc62f2f4b5f09dde621548 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 9 Jan 2024 16:34:02 +0000 Subject: [PATCH 09/15] Use structural filtering for the fixed projection traits classes --- .../Kernel_23/internal/Projection_traits_3.h | 11 ++- .../include/CGAL/Triangulation_2.h | 95 +++++++++++++++---- 2 files changed, 89 insertions(+), 17 deletions(-) diff --git a/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h b/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h index dd4e24005a48..1cc09421944d 100644 --- a/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h +++ b/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h @@ -20,6 +20,7 @@ #include #include +#include namespace CGAL { @@ -1210,6 +1211,14 @@ class Projection_traits_3 }; -} } //namespace CGAL::internal + +} // namespace internal + +template +struct Triangulation_structural_filtering_traits<::CGAL::internal::Projection_traits_3 > { + typedef typename Triangulation_structural_filtering_traits::Use_structural_filtering_tag Use_structural_filtering_tag; +}; + + } //namespace CGAL #endif // CGAL_INTERNAL_PROJECTION_TRAITS_3_H diff --git a/Triangulation_2/include/CGAL/Triangulation_2.h b/Triangulation_2/include/CGAL/Triangulation_2.h index 14e8604640f1..c0f2c7ba50e0 100644 --- a/Triangulation_2/include/CGAL/Triangulation_2.h +++ b/Triangulation_2/include/CGAL/Triangulation_2.h @@ -35,6 +35,7 @@ #include #include #include +#include #include @@ -434,6 +435,81 @@ class Triangulation_2 bool has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r) const; + + +template +inline +bool +projection_traits_has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r, + const T& ) const +{ + // So that this code works well with Lazy_kernel + internal::Static_filters_predicates::Get_approx get_approx; + + const double px = to_double(get_approx(p).x()); + const double py = to_double(get_approx(p).y()); + const double qx = to_double(get_approx(q).x()); + const double qy = to_double(get_approx(q).y()); + const double rx = to_double(get_approx(r).x()); + const double ry = to_double(get_approx(r).y()); + + const double pqx = qx - px; + const double pqy = qy - py; + const double prx = rx - px; + const double pry = ry - py; + + return ( determinant(pqx, pqy, prx, pry) < 0); +} + +template +inline +bool +projection_traits_has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r, + const ::CGAL::internal::Projection_traits_3& ) const +{ // So that this code works well with Lazy_kernel + internal::Static_filters_predicates::Get_approx get_approx; + + const double px = to_double(get_approx(p).y()); + const double py = to_double(get_approx(p).z()); + const double qx = to_double(get_approx(q).y()); + const double qy = to_double(get_approx(q).z()); + const double rx = to_double(get_approx(r).y()); + const double ry = to_double(get_approx(r).z()); + + const double pqx = qx - px; + const double pqy = qy - py; + const double prx = rx - px; + const double pry = ry - py; + + return ( determinant(pqx, pqy, prx, pry) < 0); +} + + +template +inline +bool +projection_traits_has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r, + const ::CGAL::internal::Projection_traits_3& ) const +{ // So that this code works well with Lazy_kernel + internal::Static_filters_predicates::Get_approx get_approx; + + const double px = to_double(get_approx(p).x()); + const double py = to_double(get_approx(p).z()); + const double qx = to_double(get_approx(q).x()); + const double qy = to_double(get_approx(q).z()); + const double rx = to_double(get_approx(r).x()); + const double ry = to_double(get_approx(r).z()); + + const double pqx = qx - px; + const double pqy = qy - py; + const double prx = rx - px; + const double pry = ry - py; + + return ( determinant(pqx, pqy, prx, pry) < 0); +} + + + public: Face_handle inexact_locate(const Point& p, @@ -3152,6 +3228,7 @@ inexact_locate(const Point & t, Face_handle start, int n_of_turns) const return c; } + template inline bool @@ -3159,22 +3236,8 @@ Triangulation_2:: has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r) const { - // So that this code works well with Lazy_kernel - internal::Static_filters_predicates::Get_approx get_approx; - - const double px = to_double(get_approx(p).x()); - const double py = to_double(get_approx(p).y()); - const double qx = to_double(get_approx(q).x()); - const double qy = to_double(get_approx(q).y()); - const double rx = to_double(get_approx(r).x()); - const double ry = to_double(get_approx(r).y()); - - const double pqx = qx - px; - const double pqy = qy - py; - const double prx = rx - px; - const double pry = ry - py; - - return ( determinant(pqx, pqy, prx, pry) < 0); + Gt gt; + return projection_traits_has_inexact_negative_orientation(p,q,r,gt); } #endif From 6b40f5b1894e6f494435bd9221df0508edda1108 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 5 Jan 2024 17:57:31 +0100 Subject: [PATCH 10/15] restore axis aligned projection traits --- .../Polygon_mesh_processing/autorefinement.h | 44 +++++++++++++++++-- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index c38d783572b4..94b6a9c96677 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -55,6 +55,10 @@ #endif #endif +#ifdef USE_FIXED_PROJECTION_TRAITS +#include +#endif + #include //#define CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS @@ -540,6 +544,9 @@ struct Triangle_data }; template void generate_subtriangles(std::size_t ti, Triangle_data& triangle_data, @@ -813,7 +820,12 @@ void generate_subtriangles(std::size_t ti, // init CDT + insert points and constraints CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer3.start();) +#ifdef USE_FIXED_PROJECTION_TRAITS + typedef ::CGAL::internal::Projection_traits_3 P_traits; +#else typedef CGAL::Projection_traits_3 P_traits; +#endif + typedef CGAL::No_constraint_intersection_tag Itag; typedef CGAL::Constrained_Delaunay_triangulation_2 CDT_2; @@ -821,7 +833,17 @@ void generate_subtriangles(std::size_t ti, typedef CDT_2 CDT; const std::array& t = triangles[ti]; + std::vector vhandles(triangle_data.points.size()); +#ifdef USE_FIXED_PROJECTION_TRAITS + P_traits cdt_traits; + bool orientation_flipped = false; + CDT cdt(cdt_traits); + // TODO: still need to figure out why I can't make the orientation_flipped correctly + vhandles[0]=cdt.insert(t[0]); + vhandles[1]=cdt.insert(t[1]); + vhandles[2]=cdt.insert(t[2]); +#else // positive triangle normal typename EK::Vector_3 n = normal(t[0], t[1], t[2]); typename EK::Point_3 o(CGAL::ORIGIN); @@ -832,17 +854,16 @@ void generate_subtriangles(std::size_t ti, n=-n; orientation_flipped = true; } - - std::vector vhandles(triangle_data.points.size()); P_traits cdt_traits(n); + CDT cdt(cdt_traits); vhandles[0]=cdt.insert_outside_affine_hull(t[0]); vhandles[1]=cdt.insert_outside_affine_hull(t[1]); vhandles[2] = cdt.tds().insert_dim_up(cdt.infinite_vertex(), orientation_flipped); vhandles[2]->set_point(t[2]); +#endif // insert points and fill vhandles - #if 0 std::vector indices(triangle_data.points.size()-3); std::iota(indices.begin(), indices.end(), 3); @@ -1333,7 +1354,24 @@ void autorefine_triangle_soup(PointRange& soup_points, new_triangles.push_back({triangles[ti], ti}); else { +#ifdef USE_FIXED_PROJECTION_TRAITS + const std::array& t = triangles[ti]; + typename EK::Vector_3 orth = CGAL::normal(t[0], t[1], t[2]); // TODO::avoid construction? + int c = CGAL::abs(orth[0]) > CGAL::abs(orth[1]) ? 0 : 1; + c = CGAL::abs(orth[2]) > CGAL::abs(orth[c]) ? 2 : c; + + if (c == 0) { + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + } + else if (c == 1) { + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + } + else if (c == 2) { + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + } +#else autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); +#endif } #ifdef CGAL_AUTOREF_USE_PROGRESS_DISPLAY From 7aa4f37b6ba35f99cf9594744718e6fb529f1a69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 16 Feb 2024 11:39:09 +0100 Subject: [PATCH 11/15] rename macro --- .../CGAL/Polygon_mesh_processing/autorefinement.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index 94b6a9c96677..e7108e6d521f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -55,7 +55,7 @@ #endif #endif -#ifdef USE_FIXED_PROJECTION_TRAITS +#ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS #include #endif @@ -544,7 +544,7 @@ struct Triangle_data }; template @@ -820,7 +820,7 @@ void generate_subtriangles(std::size_t ti, // init CDT + insert points and constraints CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer3.start();) -#ifdef USE_FIXED_PROJECTION_TRAITS +#ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS typedef ::CGAL::internal::Projection_traits_3 P_traits; #else typedef CGAL::Projection_traits_3 P_traits; @@ -835,7 +835,7 @@ void generate_subtriangles(std::size_t ti, const std::array& t = triangles[ti]; std::vector vhandles(triangle_data.points.size()); -#ifdef USE_FIXED_PROJECTION_TRAITS +#ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS P_traits cdt_traits; bool orientation_flipped = false; CDT cdt(cdt_traits); @@ -1354,7 +1354,7 @@ void autorefine_triangle_soup(PointRange& soup_points, new_triangles.push_back({triangles[ti], ti}); else { -#ifdef USE_FIXED_PROJECTION_TRAITS +#ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS const std::array& t = triangles[ti]; typename EK::Vector_3 orth = CGAL::normal(t[0], t[1], t[2]); // TODO::avoid construction? int c = CGAL::abs(orth[0]) > CGAL::abs(orth[1]) ? 0 : 1; From 5cae5340e99570b8eb16ac2e6a5cc0b942e59365 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 16 Feb 2024 11:48:15 +0100 Subject: [PATCH 12/15] remove TODO that are not --- .../include/CGAL/Polygon_mesh_processing/autorefinement.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index e7108e6d521f..37d9a31f4f13 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -486,7 +486,6 @@ bool collect_intersections(const std::array& t1, return false; } -//TODO: rename struct struct Triangle_data { using Point_3 = Exact_predicates_exact_constructions_kernel::Point_3; @@ -950,7 +949,6 @@ void generate_subtriangles(std::size_t ti, for (const std::pair& ids : triangle_data.segments) { - //TODO remove me CGAL_assertion(ids.first < vhandles.size()); CGAL_assertion(ids.second < vhandles.size()); CGAL_assertion( vhandles[ids.first]!= typename CDT::Vertex_handle() ); From a57800ce07afdfe18429c4922460812b628f38eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 16 Feb 2024 12:12:19 +0100 Subject: [PATCH 13/15] get rid of extra container --- .../Polygon_mesh_processing/autorefinement.h | 90 +++++++++---------- 1 file changed, 41 insertions(+), 49 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index 37d9a31f4f13..ff9a8ee4d8a3 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -255,12 +255,12 @@ do_coplanar_segments_intersect(std::size_t pi, std::size_t qi, // imported from Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h template -void coplanar_intersections(const std::array& t1, - const std::array& t2, +void coplanar_intersections(const std::vector& t1, + const std::vector& t2, std::vector>& inter_pts) { - const typename K::Point_3& p1 = t1[0], q1 = t1[1], r1 = t1[2]; - const typename K::Point_3& p2 = t2[0], q2 = t2[1], r2 = t2[2]; + const typename K::Point_3& p1 = t1[0], & q1 = t1[1], & r1 = t1[2]; + const typename K::Point_3& p2 = t2[0], & q2 = t2[1], & r2 = t2[2]; std::list> l_inter_pts; l_inter_pts.push_back(Intersections::internal::Point_on_triangle(0,-1)); @@ -416,8 +416,8 @@ test_edge(const typename K::Point_3& p, const typename K::Point_3& q, } template -bool collect_intersections(const std::array& t1, - const std::array& t2, +bool collect_intersections(const std::vector& t1, + const std::vector& t2, std::vector>& inter_pts) { // test edges of t1 vs t2 @@ -548,10 +548,9 @@ template void generate_subtriangles(std::size_t ti, - Triangle_data& triangle_data, + std::vector& all_triangle_data, const std::set >& intersecting_triangles, const std::set >& coplanar_triangles, - const std::vector>& triangles, PointVector& new_triangles ) { @@ -588,6 +587,7 @@ void generate_subtriangles(std::size_t ti, #define CGAL_AUTOREF_COUNTER_INSTRUCTION(X) #endif + Triangle_data& triangle_data=all_triangle_data[ti]; std::vector& points=triangle_data.points; std::vector& point_locations=triangle_data.point_locations; std::vector>& segments=triangle_data.segments; @@ -680,9 +680,9 @@ void generate_subtriangles(std::size_t ti, { CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer6.start();) typename EK::Point_3 pt = typename EK::Construct_planes_intersection_point_3()( - triangles[in_triangle_ids[i]][0], triangles[in_triangle_ids[i]][1],triangles[in_triangle_ids[i]][2], - triangles[in_triangle_ids[j]][0], triangles[in_triangle_ids[j]][1],triangles[in_triangle_ids[j]][2], - triangles[ti][0], triangles[ti][1],triangles[ti][2]); + all_triangle_data[in_triangle_ids[i]].points[0], all_triangle_data[in_triangle_ids[i]].points[1],all_triangle_data[in_triangle_ids[i]].points[2], + all_triangle_data[in_triangle_ids[j]].points[0], all_triangle_data[in_triangle_ids[j]].points[1],all_triangle_data[in_triangle_ids[j]].points[2], + points[0], points[1],points[2]); CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer6.stop();) CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c1;) @@ -831,7 +831,6 @@ void generate_subtriangles(std::size_t ti, //typedef CGAL::Constrained_triangulation_plus_2 CDT; typedef CDT_2 CDT; - const std::array& t = triangles[ti]; std::vector vhandles(triangle_data.points.size()); #ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS @@ -839,12 +838,12 @@ void generate_subtriangles(std::size_t ti, bool orientation_flipped = false; CDT cdt(cdt_traits); // TODO: still need to figure out why I can't make the orientation_flipped correctly - vhandles[0]=cdt.insert(t[0]); - vhandles[1]=cdt.insert(t[1]); - vhandles[2]=cdt.insert(t[2]); + vhandles[0]=cdt.insert(points[0]); + vhandles[1]=cdt.insert(points[1]); + vhandles[2]=cdt.insert(points[2]); #else // positive triangle normal - typename EK::Vector_3 n = normal(t[0], t[1], t[2]); + typename EK::Vector_3 n = normal(points[0], points[1], points[2]); typename EK::Point_3 o(CGAL::ORIGIN); bool orientation_flipped = false; @@ -856,18 +855,14 @@ void generate_subtriangles(std::size_t ti, P_traits cdt_traits(n); CDT cdt(cdt_traits); - vhandles[0]=cdt.insert_outside_affine_hull(t[0]); - vhandles[1]=cdt.insert_outside_affine_hull(t[1]); + vhandles[0]=cdt.insert_outside_affine_hull(points[0]); + vhandles[1]=cdt.insert_outside_affine_hull(points[1]); vhandles[2] = cdt.tds().insert_dim_up(cdt.infinite_vertex(), orientation_flipped); - vhandles[2]->set_point(t[2]); + vhandles[2]->set_point(points[2]); #endif // insert points and fill vhandles -#if 0 - std::vector indices(triangle_data.points.size()-3); - std::iota(indices.begin(), indices.end(), 3); -#else - //start by points on edges + // start by points on edges std::array, 3> indices_on_edges; std::vector indices; for (std::size_t i=3; i::type Pmap; typedef Spatial_sort_traits_adapter_2 Search_traits; @@ -1107,22 +1102,17 @@ void autorefine_triangle_soup(PointRange& soup_points, // init the vector of triangles used for the autorefinement of triangles typedef CGAL::Exact_predicates_exact_constructions_kernel EK; - std::vector< std::array > triangles(tiid+1); // TODO get rid of triangles and use all_triangle_data // vector of data for refining triangles - std::vector all_triangle_data(triangles.size()); + std::vector all_triangle_data(tiid+1); Cartesian_converter to_exact; for(Input_TID f : intersected_faces) { std::size_t tid=tri_inter_ids[f]; - triangles[tid]= CGAL::make_array( - to_exact( get(pm, soup_points[soup_triangles[f][0]]) ), - to_exact( get(pm, soup_points[soup_triangles[f][1]]) ), - to_exact( get(pm, soup_points[soup_triangles[f][2]]) ) ); all_triangle_data[tid].points.resize(3); - all_triangle_data[tid].points[0]=triangles[tri_inter_ids[f]][0]; - all_triangle_data[tid].points[1]=triangles[tri_inter_ids[f]][1]; - all_triangle_data[tid].points[2]=triangles[tri_inter_ids[f]][2]; + all_triangle_data[tid].points[0]=to_exact( get(pm, soup_points[soup_triangles[f][0]]) ); + all_triangle_data[tid].points[1]=to_exact( get(pm, soup_points[soup_triangles[f][1]]) ); + all_triangle_data[tid].points[2]=to_exact( get(pm, soup_points[soup_triangles[f][2]]) ); all_triangle_data[tid].point_locations.resize(3); all_triangle_data[tid].point_locations[0]=-1; all_triangle_data[tid].point_locations[1]=-2; @@ -1144,8 +1134,8 @@ void autorefine_triangle_soup(PointRange& soup_points, if (i1==-1 || i2==-1) continue; //skip degenerate faces - const std::array& t1 = triangles[i1]; - const std::array& t2 = triangles[i2]; + const std::vector& t1 = all_triangle_data[i1].points; + const std::vector& t2 = all_triangle_data[i2].points; std::vector> inter_pts; bool triangles_are_coplanar = autorefine_impl::collect_intersections(t1, t2, inter_pts); @@ -1313,7 +1303,7 @@ void autorefine_triangle_soup(PointRange& soup_points, #ifdef CGAL_LINKED_WITH_TBB if (parallel_execution) { - tbb::parallel_for(tbb::blocked_range(0, triangles.size()), + tbb::parallel_for(tbb::blocked_range(0, all_triangle_data.size()), [&](const tbb::blocked_range& r) { for (size_t ti = r.begin(); ti != r.end(); ++ti) deduplicate_inserted_segments(ti); @@ -1322,7 +1312,7 @@ void autorefine_triangle_soup(PointRange& soup_points, } else #endif - for (std::size_t ti = 0; ti < triangles.size(); ++ti) { + for (std::size_t ti = 0; ti < all_triangle_data.size(); ++ti) { deduplicate_inserted_segments(ti); } @@ -1343,32 +1333,34 @@ void autorefine_triangle_soup(PointRange& soup_points, #endif #ifdef CGAL_AUTOREF_USE_PROGRESS_DISPLAY - boost::timer::progress_display pd(triangles.size()); + boost::timer::progress_display pd(all_triangle_data.size()); #endif auto refine_triangles = [&](std::size_t ti) { if (all_triangle_data[ti].points.size()==3) - new_triangles.push_back({triangles[ti], ti}); + new_triangles.push_back({CGAL::make_array(all_triangle_data[ti].points[0], + all_triangle_data[ti].points[1], + all_triangle_data[ti].points[2]), + ti}); else { #ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS - const std::array& t = triangles[ti]; - typename EK::Vector_3 orth = CGAL::normal(t[0], t[1], t[2]); // TODO::avoid construction? + typename EK::Vector_3 orth = CGAL::normal(all_triangle_data[ti].points[0], all_triangle_data[ti].points[1], all_triangle_data[ti].points[2]); // TODO::avoid construction? int c = CGAL::abs(orth[0]) > CGAL::abs(orth[1]) ? 0 : 1; c = CGAL::abs(orth[2]) > CGAL::abs(orth[c]) ? 2 : c; if (c == 0) { - autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + autorefine_impl::generate_subtriangles(ti, all_triangle_data, intersecting_triangles, coplanar_triangles, new_triangles); } else if (c == 1) { - autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + autorefine_impl::generate_subtriangles(ti, all_triangle_data, intersecting_triangles, coplanar_triangles, new_triangles); } else if (c == 2) { - autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + autorefine_impl::generate_subtriangles(ti, all_triangle_data, intersecting_triangles, coplanar_triangles, new_triangles); } #else - autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + autorefine_impl::generate_subtriangles(ti, all_triangle_data, intersecting_triangles, coplanar_triangles, new_triangles); #endif } @@ -1383,7 +1375,7 @@ void autorefine_triangle_soup(PointRange& soup_points, #ifdef CGAL_LINKED_WITH_TBB if (parallel_execution) { - tbb::parallel_for(tbb::blocked_range(0, triangles.size()), + tbb::parallel_for(tbb::blocked_range(0, all_triangle_data.size()), [&](const tbb::blocked_range& r) { for (size_t ti = r.begin(); ti != r.end(); ++ti) refine_triangles(ti); @@ -1392,7 +1384,7 @@ void autorefine_triangle_soup(PointRange& soup_points, } else #endif - for (std::size_t ti = 0; ti < triangles.size(); ++ti) { + for (std::size_t ti = 0; ti < all_triangle_data.size(); ++ti) { refine_triangles(ti); } @@ -1449,7 +1441,7 @@ void autorefine_triangle_soup(PointRange& soup_points, } // raw copy of input triangles with no intersection - std::vector tri_inter_ids_inverse(triangles.size()); + std::vector tri_inter_ids_inverse(all_triangle_data.size()); for (Input_TID f=0; f Date: Mon, 26 Feb 2024 14:18:01 +0100 Subject: [PATCH 14/15] Revert "get rid of extra container" This reverts commit a57800ce07afdfe18429c4922460812b628f38eb. points vector of Triangle_data can be updated when intersecting new intersection points. If the container is used at the same time by another thread calling generate_subtriangles, the container might be in an invalid state while resizing it. --- .../Polygon_mesh_processing/autorefinement.h | 92 ++++++++++--------- 1 file changed, 51 insertions(+), 41 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index ff9a8ee4d8a3..e9d85736f9d4 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -255,12 +255,12 @@ do_coplanar_segments_intersect(std::size_t pi, std::size_t qi, // imported from Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h template -void coplanar_intersections(const std::vector& t1, - const std::vector& t2, +void coplanar_intersections(const std::array& t1, + const std::array& t2, std::vector>& inter_pts) { - const typename K::Point_3& p1 = t1[0], & q1 = t1[1], & r1 = t1[2]; - const typename K::Point_3& p2 = t2[0], & q2 = t2[1], & r2 = t2[2]; + const typename K::Point_3& p1 = t1[0], q1 = t1[1], r1 = t1[2]; + const typename K::Point_3& p2 = t2[0], q2 = t2[1], r2 = t2[2]; std::list> l_inter_pts; l_inter_pts.push_back(Intersections::internal::Point_on_triangle(0,-1)); @@ -416,8 +416,8 @@ test_edge(const typename K::Point_3& p, const typename K::Point_3& q, } template -bool collect_intersections(const std::vector& t1, - const std::vector& t2, +bool collect_intersections(const std::array& t1, + const std::array& t2, std::vector>& inter_pts) { // test edges of t1 vs t2 @@ -548,9 +548,10 @@ template void generate_subtriangles(std::size_t ti, - std::vector& all_triangle_data, + Triangle_data& triangle_data, const std::set >& intersecting_triangles, const std::set >& coplanar_triangles, + const std::vector>& triangles, PointVector& new_triangles ) { @@ -587,7 +588,6 @@ void generate_subtriangles(std::size_t ti, #define CGAL_AUTOREF_COUNTER_INSTRUCTION(X) #endif - Triangle_data& triangle_data=all_triangle_data[ti]; std::vector& points=triangle_data.points; std::vector& point_locations=triangle_data.point_locations; std::vector>& segments=triangle_data.segments; @@ -680,9 +680,9 @@ void generate_subtriangles(std::size_t ti, { CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer6.start();) typename EK::Point_3 pt = typename EK::Construct_planes_intersection_point_3()( - all_triangle_data[in_triangle_ids[i]].points[0], all_triangle_data[in_triangle_ids[i]].points[1],all_triangle_data[in_triangle_ids[i]].points[2], - all_triangle_data[in_triangle_ids[j]].points[0], all_triangle_data[in_triangle_ids[j]].points[1],all_triangle_data[in_triangle_ids[j]].points[2], - points[0], points[1],points[2]); + triangles[in_triangle_ids[i]][0], triangles[in_triangle_ids[i]][1],triangles[in_triangle_ids[i]][2], + triangles[in_triangle_ids[j]][0], triangles[in_triangle_ids[j]][1],triangles[in_triangle_ids[j]][2], + triangles[ti][0], triangles[ti][1],triangles[ti][2]); CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer6.stop();) CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c1;) @@ -831,6 +831,7 @@ void generate_subtriangles(std::size_t ti, //typedef CGAL::Constrained_triangulation_plus_2 CDT; typedef CDT_2 CDT; + const std::array& t = triangles[ti]; std::vector vhandles(triangle_data.points.size()); #ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS @@ -838,12 +839,12 @@ void generate_subtriangles(std::size_t ti, bool orientation_flipped = false; CDT cdt(cdt_traits); // TODO: still need to figure out why I can't make the orientation_flipped correctly - vhandles[0]=cdt.insert(points[0]); - vhandles[1]=cdt.insert(points[1]); - vhandles[2]=cdt.insert(points[2]); + vhandles[0]=cdt.insert(t[0]); + vhandles[1]=cdt.insert(t[1]); + vhandles[2]=cdt.insert(t[2]); #else // positive triangle normal - typename EK::Vector_3 n = normal(points[0], points[1], points[2]); + typename EK::Vector_3 n = normal(t[0], t[1], t[2]); typename EK::Point_3 o(CGAL::ORIGIN); bool orientation_flipped = false; @@ -855,14 +856,18 @@ void generate_subtriangles(std::size_t ti, P_traits cdt_traits(n); CDT cdt(cdt_traits); - vhandles[0]=cdt.insert_outside_affine_hull(points[0]); - vhandles[1]=cdt.insert_outside_affine_hull(points[1]); + vhandles[0]=cdt.insert_outside_affine_hull(t[0]); + vhandles[1]=cdt.insert_outside_affine_hull(t[1]); vhandles[2] = cdt.tds().insert_dim_up(cdt.infinite_vertex(), orientation_flipped); - vhandles[2]->set_point(points[2]); + vhandles[2]->set_point(t[2]); #endif // insert points and fill vhandles - // start by points on edges +#if 0 + std::vector indices(triangle_data.points.size()-3); + std::iota(indices.begin(), indices.end(), 3); +#else + //start by points on edges std::array, 3> indices_on_edges; std::vector indices; for (std::size_t i=3; i::type Pmap; typedef Spatial_sort_traits_adapter_2 Search_traits; @@ -1102,17 +1107,24 @@ void autorefine_triangle_soup(PointRange& soup_points, // init the vector of triangles used for the autorefinement of triangles typedef CGAL::Exact_predicates_exact_constructions_kernel EK; + // even if the info is duplicated with Triangle_data::point, we keep this container + // so that we can use it in parallel calls of generate_subtriangles + std::vector< std::array > triangles(tiid+1); // vector of data for refining triangles - std::vector all_triangle_data(tiid+1); + std::vector all_triangle_data(triangles.size()); Cartesian_converter to_exact; for(Input_TID f : intersected_faces) { std::size_t tid=tri_inter_ids[f]; + triangles[tid]= CGAL::make_array( + to_exact( get(pm, soup_points[soup_triangles[f][0]]) ), + to_exact( get(pm, soup_points[soup_triangles[f][1]]) ), + to_exact( get(pm, soup_points[soup_triangles[f][2]]) ) ); all_triangle_data[tid].points.resize(3); - all_triangle_data[tid].points[0]=to_exact( get(pm, soup_points[soup_triangles[f][0]]) ); - all_triangle_data[tid].points[1]=to_exact( get(pm, soup_points[soup_triangles[f][1]]) ); - all_triangle_data[tid].points[2]=to_exact( get(pm, soup_points[soup_triangles[f][2]]) ); + all_triangle_data[tid].points[0]=triangles[tri_inter_ids[f]][0]; + all_triangle_data[tid].points[1]=triangles[tri_inter_ids[f]][1]; + all_triangle_data[tid].points[2]=triangles[tri_inter_ids[f]][2]; all_triangle_data[tid].point_locations.resize(3); all_triangle_data[tid].point_locations[0]=-1; all_triangle_data[tid].point_locations[1]=-2; @@ -1134,8 +1146,8 @@ void autorefine_triangle_soup(PointRange& soup_points, if (i1==-1 || i2==-1) continue; //skip degenerate faces - const std::vector& t1 = all_triangle_data[i1].points; - const std::vector& t2 = all_triangle_data[i2].points; + const std::array& t1 = triangles[i1]; + const std::array& t2 = triangles[i2]; std::vector> inter_pts; bool triangles_are_coplanar = autorefine_impl::collect_intersections(t1, t2, inter_pts); @@ -1303,7 +1315,7 @@ void autorefine_triangle_soup(PointRange& soup_points, #ifdef CGAL_LINKED_WITH_TBB if (parallel_execution) { - tbb::parallel_for(tbb::blocked_range(0, all_triangle_data.size()), + tbb::parallel_for(tbb::blocked_range(0, triangles.size()), [&](const tbb::blocked_range& r) { for (size_t ti = r.begin(); ti != r.end(); ++ti) deduplicate_inserted_segments(ti); @@ -1312,7 +1324,7 @@ void autorefine_triangle_soup(PointRange& soup_points, } else #endif - for (std::size_t ti = 0; ti < all_triangle_data.size(); ++ti) { + for (std::size_t ti = 0; ti < triangles.size(); ++ti) { deduplicate_inserted_segments(ti); } @@ -1333,34 +1345,32 @@ void autorefine_triangle_soup(PointRange& soup_points, #endif #ifdef CGAL_AUTOREF_USE_PROGRESS_DISPLAY - boost::timer::progress_display pd(all_triangle_data.size()); + boost::timer::progress_display pd(triangles.size()); #endif auto refine_triangles = [&](std::size_t ti) { if (all_triangle_data[ti].points.size()==3) - new_triangles.push_back({CGAL::make_array(all_triangle_data[ti].points[0], - all_triangle_data[ti].points[1], - all_triangle_data[ti].points[2]), - ti}); + new_triangles.push_back({triangles[ti], ti}); else { #ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS - typename EK::Vector_3 orth = CGAL::normal(all_triangle_data[ti].points[0], all_triangle_data[ti].points[1], all_triangle_data[ti].points[2]); // TODO::avoid construction? + const std::array& t = triangles[ti]; + typename EK::Vector_3 orth = CGAL::normal(t[0], t[1], t[2]); // TODO::avoid construction? int c = CGAL::abs(orth[0]) > CGAL::abs(orth[1]) ? 0 : 1; c = CGAL::abs(orth[2]) > CGAL::abs(orth[c]) ? 2 : c; if (c == 0) { - autorefine_impl::generate_subtriangles(ti, all_triangle_data, intersecting_triangles, coplanar_triangles, new_triangles); + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); } else if (c == 1) { - autorefine_impl::generate_subtriangles(ti, all_triangle_data, intersecting_triangles, coplanar_triangles, new_triangles); + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); } else if (c == 2) { - autorefine_impl::generate_subtriangles(ti, all_triangle_data, intersecting_triangles, coplanar_triangles, new_triangles); + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); } #else - autorefine_impl::generate_subtriangles(ti, all_triangle_data, intersecting_triangles, coplanar_triangles, new_triangles); + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); #endif } @@ -1375,7 +1385,7 @@ void autorefine_triangle_soup(PointRange& soup_points, #ifdef CGAL_LINKED_WITH_TBB if (parallel_execution) { - tbb::parallel_for(tbb::blocked_range(0, all_triangle_data.size()), + tbb::parallel_for(tbb::blocked_range(0, triangles.size()), [&](const tbb::blocked_range& r) { for (size_t ti = r.begin(); ti != r.end(); ++ti) refine_triangles(ti); @@ -1384,7 +1394,7 @@ void autorefine_triangle_soup(PointRange& soup_points, } else #endif - for (std::size_t ti = 0; ti < all_triangle_data.size(); ++ti) { + for (std::size_t ti = 0; ti < triangles.size(); ++ti) { refine_triangles(ti); } @@ -1441,7 +1451,7 @@ void autorefine_triangle_soup(PointRange& soup_points, } // raw copy of input triangles with no intersection - std::vector tri_inter_ids_inverse(all_triangle_data.size()); + std::vector tri_inter_ids_inverse(triangles.size()); for (Input_TID f=0; f Date: Mon, 26 Feb 2024 14:23:09 +0100 Subject: [PATCH 15/15] missing references --- .../include/CGAL/Polygon_mesh_processing/autorefinement.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index e9d85736f9d4..8de62f6b8021 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -259,8 +259,8 @@ void coplanar_intersections(const std::array& t1, const std::array& t2, std::vector>& inter_pts) { - const typename K::Point_3& p1 = t1[0], q1 = t1[1], r1 = t1[2]; - const typename K::Point_3& p2 = t2[0], q2 = t2[1], r2 = t2[2]; + const typename K::Point_3& p1 = t1[0], &q1 = t1[1], &r1 = t1[2]; + const typename K::Point_3& p2 = t2[0], &q2 = t2[1], &r2 = t2[2]; std::list> l_inter_pts; l_inter_pts.push_back(Intersections::internal::Point_on_triangle(0,-1));