From b8e3a565e63e616d742814d76dbff5e5b5afa412 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 8 May 2020 13:28:34 +0200 Subject: [PATCH 01/27] Fix incorrect NP for clipping output --- .../include/CGAL/Polygon_mesh_processing/clip.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h index d059c3eb33d..01777e9740c 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h @@ -465,10 +465,9 @@ clip(TriangleMesh& tm, parameters::choose_parameter(parameters::get_parameter(np_tm, internal_np::clip_volume), false); if(clip_volume && is_closed(tm)) - return corefine_and_compute_intersection(tm, clipper, tm, np_tm, np_c); + return corefine_and_compute_intersection(tm, clipper, tm, np_tm, np_c, np_tm); return corefine_and_compute_intersection(tm, clipper, tm, - np_tm.use_bool_op_to_clip_surface(true), - np_c); + np_tm.use_bool_op_to_clip_surface(true), np_c, np_tm); } /** From 1c5acc018a169aa7fff149abe731dde3975d9c94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 8 May 2020 13:34:22 +0200 Subject: [PATCH 02/27] Enable PMP::Corefinement to work with VPMs of different types Same value_type is still a precondition, though. --- .../Polygon_mesh_processing/corefinement.h | 99 ++- .../Corefinement/Face_graph_output_builder.h | 24 +- .../internal/Corefinement/Visitor.h | 714 +++++++++--------- .../internal/Corefinement/face_graph_utils.h | 35 +- .../intersect_triangle_and_segment_3.h | 15 +- .../Corefinement/intersection_callbacks.h | 12 +- .../internal/Corefinement/intersection_impl.h | 172 +++-- .../Corefinement/intersection_nodes.h | 86 ++- .../intersection_of_coplanar_triangles_3.h | 27 +- .../internal/Corefinement/predicates.h | 12 +- .../Polygon_mesh_processing/intersection.h | 36 +- 11 files changed, 645 insertions(+), 587 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h index e0e8c27de2f..fdc6e59fe6d 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h @@ -174,22 +174,19 @@ corefine_and_compute_boolean_operations( // Vertex point maps //for input meshes - typedef typename GetVertexPointMap::type Vpm; - typedef typename GetVertexPointMap::type Vpm2; - CGAL_USE_TYPE(Vpm2); - CGAL_assertion_code( - static const bool same_vpm = (boost::is_same::value); ) - CGAL_static_assertion(same_vpm); + typedef typename GetVertexPointMap::type VPM1; + typedef typename GetVertexPointMap::type VPM2; - Vpm vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point), - get_property_map(boost::vertex_point, tm1)); + CGAL_static_assertion((std::is_same::value_type, + typename boost::property_traits::value_type>::value)); - Vpm vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point), - get_property_map(boost::vertex_point, tm2)); + VPM1 vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point), + get_property_map(boost::vertex_point, tm1)); - typedef typename boost::property_traits::value_type Point_3; + VPM2 vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point), + get_property_map(boost::vertex_point, tm2)); + + typedef typename boost::property_traits::value_type Point_3; // for output meshes: here we have to use a trick so that if for a specific output // that is not requested, the default vpm does not have the same value type as the @@ -200,24 +197,24 @@ corefine_and_compute_boolean_operations( Corefinement::TweakedGetVertexPointMap, Corefinement::TweakedGetVertexPointMap, Corefinement::TweakedGetVertexPointMap - > Vpm_out_tuple_helper; + > VPM_out_tuple_helper; typedef std::tuple< - boost::optional< typename std::tuple_element<0, Vpm_out_tuple_helper>::type::type >, - boost::optional< typename std::tuple_element<1, Vpm_out_tuple_helper>::type::type >, - boost::optional< typename std::tuple_element<2, Vpm_out_tuple_helper>::type::type >, - boost::optional< typename std::tuple_element<3, Vpm_out_tuple_helper>::type::type > - > Vpm_out_tuple; + boost::optional< typename std::tuple_element<0, VPM_out_tuple_helper>::type::type >, + boost::optional< typename std::tuple_element<1, VPM_out_tuple_helper>::type::type >, + boost::optional< typename std::tuple_element<2, VPM_out_tuple_helper>::type::type >, + boost::optional< typename std::tuple_element<3, VPM_out_tuple_helper>::type::type > + > VPM_out_tuple; - Vpm_out_tuple vpm_out_tuple( + VPM_out_tuple vpm_out_tuple( Corefinement::get_vpm(std::get<0>(nps_out), output[0], - typename std::tuple_element<0, Vpm_out_tuple_helper>::type::Use_default_tag()), + typename std::tuple_element<0, VPM_out_tuple_helper>::type::Use_default_tag()), Corefinement::get_vpm(std::get<1>(nps_out), output[1], - typename std::tuple_element<1, Vpm_out_tuple_helper>::type::Use_default_tag()), + typename std::tuple_element<1, VPM_out_tuple_helper>::type::Use_default_tag()), Corefinement::get_vpm(std::get<2>(nps_out), output[2], - typename std::tuple_element<2, Vpm_out_tuple_helper>::type::Use_default_tag()), + typename std::tuple_element<2, VPM_out_tuple_helper>::type::Use_default_tag()), Corefinement::get_vpm(std::get<3>(nps_out), output[3], - typename std::tuple_element<3, Vpm_out_tuple_helper>::type::Use_default_tag()) + typename std::tuple_element<3, VPM_out_tuple_helper>::type::Use_default_tag()) ); if (&tm1==&tm2) @@ -348,8 +345,9 @@ corefine_and_compute_boolean_operations( // surface intersection algorithm call typedef Corefinement::Face_graph_output_builder Ob; typedef Corefinement::Surface_intersection_visitor_for_corefinement< - TriangleMesh, Vpm, Ob, Ecm_in, User_visitor> Algo_visitor; + TriangleMesh, VPM1, VPM2, Ob, Ecm_in, User_visitor> Algo_visitor; + Ecm_in ecm_in(tm1,tm2,ecm1,ecm2); Edge_mark_map_tuple ecms_out(ecm_out_0, ecm_out_1, ecm_out_2, ecm_out_3); Ob ob(tm1, tm2, vpm1, vpm2, fid_map1, fid_map2, ecm_in, vpm_out_tuple, ecms_out, uv, output); @@ -375,7 +374,7 @@ corefine_and_compute_boolean_operations( ob.setup_for_clipping_a_surface(use_compact_clipper); } - Corefinement::Intersection_of_triangle_meshes + Corefinement::Intersection_of_triangle_meshes functor(tm1, tm2, vpm1, vpm2, Algo_visitor(uv,ob,ecm_in)); functor(CGAL::Emptyset_iterator(), throw_on_self_intersection, true); @@ -645,20 +644,17 @@ corefine( TriangleMesh& tm1, choose_parameter(get_parameter(np1, internal_np::throw_on_self_intersection), false); // Vertex point maps - typedef typename GetVertexPointMap::type Vpm; - typedef typename GetVertexPointMap::type Vpm2; - CGAL_USE_TYPE(Vpm2); - CGAL_assertion_code( - static const bool same_vpm = (boost::is_same::value);) - CGAL_static_assertion(same_vpm); + typedef typename GetVertexPointMap::type VPM1; + typedef typename GetVertexPointMap::type VPM2; - Vpm vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point), - get_property_map(boost::vertex_point, tm1)); + CGAL_static_assertion((std::is_same::value_type, + typename boost::property_traits::value_type>::value)); - Vpm vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point), - get_property_map(boost::vertex_point, tm2)); + VPM1 vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point), + get_property_map(boost::vertex_point, tm1)); + + VPM2 vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point), + get_property_map(boost::vertex_point, tm2)); // Edge is-constrained maps typedef typename internal_np::Lookup_named_param_def < @@ -696,10 +692,11 @@ corefine( TriangleMesh& tm1, // surface intersection algorithm call typedef Corefinement::No_extra_output_from_corefinement Ob; typedef Corefinement::Surface_intersection_visitor_for_corefinement< - TriangleMesh, Vpm, Ob, Ecm, User_visitor> Algo_visitor; + TriangleMesh, VPM1, VPM2, Ob, Ecm, User_visitor> Algo_visitor; + Ob ob; Ecm ecm(tm1,tm2,ecm1,ecm2); - Corefinement::Intersection_of_triangle_meshes + Corefinement::Intersection_of_triangle_meshes functor(tm1, tm2, vpm1, vpm2, Algo_visitor(uv,ob,ecm)); functor(CGAL::Emptyset_iterator(), throw_on_self_intersection, true); } @@ -743,9 +740,9 @@ autorefine( TriangleMesh& tm, using parameters::get_parameter; // Vertex point maps - typedef typename GetVertexPointMap::type Vpm; + typedef typename GetVertexPointMap::type VPM; - Vpm vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), get_property_map(boost::vertex_point, tm)); // Edge is-constrained maps @@ -768,10 +765,10 @@ autorefine( TriangleMesh& tm, // surface intersection algorithm call typedef Corefinement::No_extra_output_from_corefinement Ob; typedef Corefinement::Surface_intersection_visitor_for_corefinement< - TriangleMesh, Vpm, Ob, Ecm, User_visitor,true> Algo_visitor; + TriangleMesh, VPM, VPM, Ob, Ecm, User_visitor,true> Algo_visitor; Ob ob; - Corefinement::Intersection_of_triangle_meshes + Corefinement::Intersection_of_triangle_meshes functor(tm, vpm, Algo_visitor(uv,ob,ecm) ); functor(CGAL::Emptyset_iterator(), true); @@ -817,8 +814,8 @@ autorefine_and_remove_self_intersections( TriangleMesh& tm, using parameters::get_parameter; // Vertex point maps - typedef typename GetVertexPointMap::type Vpm; - Vpm vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + typedef typename GetVertexPointMap::type VPM; + VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), get_property_map(boost::vertex_point, tm)); // Face index map @@ -844,16 +841,16 @@ autorefine_and_remove_self_intersections( TriangleMesh& tm, // surface intersection algorithm call typedef Corefinement::Output_builder_for_autorefinement Ob; typedef Corefinement::Surface_intersection_visitor_for_corefinement< - TriangleMesh, Vpm, Ob, Ecm, User_visitor,true> Algo_visitor; + TriangleMesh, VPM, VPM, Ob, Ecm, User_visitor,true> Algo_visitor; Ob ob(tm, vpm, fid_map, ecm); - Corefinement::Intersection_of_triangle_meshes + Corefinement::Intersection_of_triangle_meshes functor(tm, vpm, Algo_visitor(uv,ob,ecm) ); functor(CGAL::Emptyset_iterator(), true); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h index 325a3693d5b..6830a916731 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h @@ -58,7 +58,8 @@ namespace PMP=Polygon_mesh_processing; namespace params=PMP::parameters; template ::value_type - >::Kernel >::type Kernel; + typename boost::property_traits::value_type + >::Kernel >::type Kernel; + typedef typename Default::Get > >::type EdgeMarkMapBind; @@ -111,8 +113,8 @@ class Face_graph_output_builder //Data members TriangleMesh &tm1, &tm2; // property maps of input meshes - const VertexPointMap vpm1; - const VertexPointMap vpm2; + const VertexPointMap1& vpm1; + const VertexPointMap2& vpm2; FaceIdMap1 fids1; FaceIdMap2 fids2; EdgeMarkMapBind& marks_on_input_edges; @@ -346,8 +348,8 @@ public: Face_graph_output_builder(TriangleMesh& tm1, TriangleMesh& tm2, - const VertexPointMap vpm1, - const VertexPointMap vpm2, + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2, FaceIdMap1 fids1, FaceIdMap2 fids2, EdgeMarkMapBind& marks_on_input_edges, @@ -1065,10 +1067,6 @@ public: if (!is_tm2_closed) patch_status_not_set_tm1.reset(); - typedef Side_of_triangle_mesh Inside_poly_test; - #ifdef CGAL_COREFINEMENT_POLYHEDRA_DEBUG #warning stop using next_marked_halfedge_around_target_vertex and create lists of halfedges instead? #endif @@ -1078,7 +1076,7 @@ public: CGAL::Bounded_side in_tm2 = is_tm2_inside_out ? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE; - Inside_poly_test inside_tm2(tm2, vpm2); + Side_of_triangle_mesh inside_tm2(tm2, vpm2); for(face_descriptor f : faces(tm1)) { @@ -1140,7 +1138,7 @@ public: CGAL::Bounded_side in_tm1 = is_tm1_inside_out ? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE; - Inside_poly_test inside_tm1(tm1, vpm1); + Side_of_triangle_mesh inside_tm1(tm1, vpm1); for(face_descriptor f : faces(tm2)) { const std::size_t f_id = get(fids2, f); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h index d668782b25c..f2bb02c849a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h @@ -102,7 +102,8 @@ struct No_extra_output_from_corefinement // A visitor for Intersection_of_triangle_meshes that can be used to corefine // two meshes template< class TriangleMesh, - class VertexPointMap, + class VertexPointMap1, + class VertexPointMap2, class OutputBuilder_ = Default, class EdgeMarkMapBind_ = Default, class UserVisitor_ = Default, @@ -141,8 +142,9 @@ private: typedef boost::unordered_map Vertex_to_node_id; typedef std::map Mesh_to_vertex_to_node_id; // typedef for the CDT - typedef typename Intersection_nodes::Exact_kernel EK; + typedef Intersection_nodes INodes; + typedef typename INodes::Exact_kernel EK; typedef Triangulation_2_projection_traits_3 CDT_traits; typedef Triangulation_vertex_base_with_info_2 Vb; typedef Constrained_triangulation_face_base_2 Fb; @@ -399,16 +401,16 @@ public: //sort node ids so that we can split the hedge //consecutively - template + template void sort_vertices_along_hedge(std::vector& node_ids, halfedge_descriptor hedge, const TriangleMesh& tm, - const VertexPointMap& vpm, + const VPM& vpm, const Node_vector& nodes) { std::sort(node_ids.begin(), node_ids.end(), - Less_along_a_halfedge + Less_along_a_halfedge (hedge, tm, vpm, nodes) ); } @@ -483,6 +485,8 @@ public: }; + typedef boost::unordered_map Face_boundaries; + //update the id of input mesh vertex that are also a node void update_face_indices( std::array& f_vertices, @@ -585,22 +589,362 @@ public: return vh; } - void finalize(Intersection_nodes& nodes, - const TriangleMesh& tm1, - const TriangleMesh& tm2, - const VertexPointMap& vpm1, - const VertexPointMap& vpm2) + template + void split_halfedges(OnEdgeMapIterator it, + const VPM& vpm, + INodes& nodes, + std::map& mesh_to_face_boundaries) + { + TriangleMesh& tm=*it->first; + On_edge_map& on_edge_map=it->second; + On_face_map& on_face_map=on_face[&tm]; + Face_boundaries& face_boundaries=mesh_to_face_boundaries[&tm]; + + for(typename On_edge_map::iterator it2=on_edge_map.begin(); + it2!=on_edge_map.end(); + ++it2) + { + //the edge to be split + halfedge_descriptor hedge=halfedge(it2->first,tm); + //indices of the nodes to be inserted + Node_ids& node_ids=it2->second; + CGAL_assertion( std::set(node_ids.begin(), node_ids.end()) + .size()==node_ids.size() ); + //sort nodes along the egde to allow consecutive splits + sort_vertices_along_hedge(node_ids,hedge,tm,vpm,nodes); + + //save original face and nodes for face of hedge (1) + if ( !is_border(hedge,tm) ){ + face_descriptor f=face(hedge,tm); + typename Face_boundaries::iterator it_face = face_boundaries.find(f); + if (it_face==face_boundaries.end()) + it_face=face_boundaries.insert(std::make_pair(f,Face_boundary(hedge,tm))).first; + it_face->second.copy_node_ids(hedge,node_ids.begin(),node_ids.end()); + } + + //save original face and nodes for face of hedge->opposite (2) + typename Face_boundaries::iterator opposite_original_info=face_boundaries.end(); + halfedge_descriptor hedge_opp = opposite(hedge,tm); + if ( !is_border(hedge_opp,tm) ){ + face_descriptor f=face(hedge_opp,tm); + opposite_original_info=face_boundaries.find(f); + if (opposite_original_info==face_boundaries.end()) + opposite_original_info=face_boundaries.insert(std::make_pair(f,Face_boundary(hedge_opp,tm))).first; + opposite_original_info->second.copy_node_ids(hedge_opp,node_ids.rbegin(),node_ids.rend()); + } + + typename Mesh_to_map_node::iterator it_map=mesh_to_node_id_to_vertex.find(&tm); + CGAL_assertion(it_map!=mesh_to_node_id_to_vertex.end()); + //a map to identify the vertex in the polyhedron corresponding to an intersection point + Node_id_to_vertex& node_id_to_vertex=it_map->second; + + CGAL_assertion_code(vertex_descriptor original_vertex=source(hedge,tm);) + + //We need an edge incident to the source vertex of hedge. This is the first opposite edge created. + bool first=true; + halfedge_descriptor hedge_incident_to_src=Graph_traits::null_halfedge(); + bool hedge_is_marked = call_get(marks_on_edges,tm,edge(hedge,tm)); + //do split the edges + CGAL_assertion_code(vertex_descriptor expected_src=source(hedge,tm)); + for(std::size_t node_id : node_ids) + { + halfedge_descriptor hnew = Euler::split_edge(hedge, tm); + CGAL_assertion(expected_src==source(hnew,tm)); + vertex_descriptor vnew=target(hnew,tm); +// user_visitor.new_vertex_added(node_id, vnew, tm); // NODE_VISITOR_TAG + nodes.call_put(vpm, vnew, node_id, tm); + // register the new vertex in the output builder + output_builder.set_vertex_id(vnew, node_id, tm); + node_id_to_vertex[node_id]=vnew; + if (first){ + first=false; + hedge_incident_to_src=next(opposite(hedge,tm),tm); + } + + //update marker tags. If the edge was marked, then the resulting edges in the split must be marked + if ( hedge_is_marked ) + call_put(marks_on_edges,tm,edge(hnew,tm),true); + + CGAL_assertion_code(expected_src=vnew); + } + + CGAL_assertion(target(hedge_incident_to_src,tm)==original_vertex); + CGAL_assertion(face(hedge_incident_to_src,tm)==face(hedge_opp,tm)); + + //save original face and nodes for face of hedge->opposite (2) + if ( !is_border(hedge_opp,tm) ){ + CGAL_assertion(opposite_original_info!=face_boundaries.end()); + opposite_original_info->second.update_original_halfedge( + hedge_opp,hedge_incident_to_src,tm); + } + + //insert the two incident faces in on_face map so that they will be triangulated. + if (!is_border(hedge,tm)) on_face_map[face(hedge,tm)]; + if (!is_border(hedge_opp,tm)) on_face_map[face(hedge_opp,tm)]; + } + } + + template + void triangulate_intersected_faces(OnFaceMapIterator it, + const VPM& vpm, + INodes& nodes, + std::map& mesh_to_face_boundaries) + { + TriangleMesh& tm=*it->first; + On_face_map& on_face_map=it->second; + Face_boundaries& face_boundaries=mesh_to_face_boundaries[&tm]; + Node_id_to_vertex& node_id_to_vertex=mesh_to_node_id_to_vertex[&tm]; + Vertex_to_node_id& vertex_to_node_id=mesh_to_vertex_to_node_id[&tm]; + + const Node_id nb_nodes = nodes.size(); + + for (typename On_face_map::iterator it=on_face_map.begin(); + it!=on_face_map.end();++it) + { + face_descriptor f = it->first; //the face to be triangulated + Node_ids& node_ids = it->second; // ids of nodes in the interior of f + typename Face_boundaries::iterator it_fb=face_boundaries.find(f); + + std::map id_to_CDT_vh; + + //associate an edge of the triangulation to a halfedge in a given polyhedron + std::map,halfedge_descriptor> edge_to_hedge; + + // the vertices of f + std::array f_vertices; + // the node_id of an input vertex or a fake id (>=nb_nodes) + std::array f_indices = {{nb_nodes,nb_nodes+1,nb_nodes+2}}; + if (it_fb!=face_boundaries.end()){ //the boundary of the triangle face was refined + f_vertices[0]=it_fb->second.vertices[0]; + f_vertices[1]=it_fb->second.vertices[1]; + f_vertices[2]=it_fb->second.vertices[2]; + update_face_indices(f_vertices,f_indices,vertex_to_node_id); + if (doing_autorefinement) + it_fb->second.update_node_id_to_vertex_map(node_id_to_vertex, tm); + } + else{ + CGAL_assertion( is_triangle(halfedge(f,tm),tm) ); + halfedge_descriptor h0=halfedge(f,tm), h1=next(h0,tm), h2=next(h1,tm); + f_vertices[0]=target(h0,tm); //nb_nodes + f_vertices[1]=target(h1,tm); //nb_nodes+1 + f_vertices[2]=target(h2,tm); //nb_nodes+2 + + update_face_indices(f_vertices,f_indices,vertex_to_node_id); + edge_to_hedge[std::make_pair( f_indices[2],f_indices[0] )] = h0; + edge_to_hedge[std::make_pair( f_indices[0],f_indices[1] )] = h1; + edge_to_hedge[std::make_pair( f_indices[1],f_indices[2] )] = h2; + } + + typename EK::Point_3 p = nodes.to_exact(get(vpm,f_vertices[0])), + q = nodes.to_exact(get(vpm,f_vertices[1])), + r = nodes.to_exact(get(vpm,f_vertices[2])); +///TODO use a positive normal and remove all work around to guarantee that triangulation of coplanar patches are compatible + CDT_traits traits(typename EK::Construct_normal_3()(p,q,r)); + CDT cdt(traits); + + // insert triangle points + std::array triangle_vertices; + //we can do this to_exact because these are supposed to be input points. + triangle_vertices[0]=cdt.insert_outside_affine_hull(p); + triangle_vertices[1]=cdt.insert_outside_affine_hull(q); + triangle_vertices[2]=cdt.tds().insert_dim_up(cdt.infinite_vertex(), false); + triangle_vertices[2]->set_point(r); + + + triangle_vertices[0]->info()=f_indices[0]; + triangle_vertices[1]->info()=f_indices[1]; + triangle_vertices[2]->info()=f_indices[2]; + + node_id_to_vertex[nb_nodes ]=f_vertices[0]; + node_id_to_vertex[nb_nodes+1]=f_vertices[1]; + node_id_to_vertex[nb_nodes+2]=f_vertices[2]; + + //if one of the triangle input vertex is also a node + for (int ik=0;ik<3;++ik){ + if ( f_indices[ik]vertex(oi) ) ); + } + + // In this loop, for each original edge of the triangle, we insert + // the constrained edges and we recover the halfedge_descriptor + // corresponding to these constrained (they are already in tm) + Face_boundary& f_boundary=it_fb->second; + for (int i=0;i<3;++i){ + //handle case of halfedge starting at triangle_vertices[i] + // and ending at triangle_vertices[(i+1)%3] + + const Node_ids& ids_on_edge=f_boundary.node_ids_array[i]; + CDT_Vertex_handle previous=triangle_vertices[i]; + Node_id prev_index=f_indices[i];// node-id of the mesh vertex + halfedge_descriptor hedge = next(f_boundary.halfedges[(i+2)%3],tm); + CGAL_assertion( source(hedge,tm)==f_boundary.vertices[i] ); + if (!ids_on_edge.empty()){ //is there at least one node on this edge? + // fh must be an infinite face + // The points must be ordered from fh->vertex(cw(infinite_vertex)) to fh->vertex(ccw(infinite_vertex)) + for(Node_id id : ids_on_edge) + { + CDT_Vertex_handle vh=insert_point_on_ch_edge(cdt,infinite_faces[i],nodes.exact_node(id)); + vh->info()=id; + id_to_CDT_vh.insert(std::make_pair(id,vh)); + edge_to_hedge[std::make_pair(prev_index,id)]=hedge; + previous=vh; + hedge=next(hedge,tm); + prev_index=id; + } + } + else{ + CGAL_assertion_code(halfedge_descriptor hd=f_boundary.halfedges[i]); + CGAL_assertion( target(hd,tm) == f_boundary.vertices[(i+1)%3] ); + CGAL_assertion( source(hd,tm) == f_boundary.vertices[ i ] ); + } + CGAL_assertion(hedge==f_boundary.halfedges[i]); + edge_to_hedge[std::make_pair(prev_index,f_indices[(i+1)%3])] = + it_fb->second.halfedges[i]; + } + } + + //insert point inside face + for(Node_id node_id : node_ids) + { + CDT_Vertex_handle vh=cdt.insert(nodes.exact_node(node_id)); + vh->info()=node_id; + id_to_CDT_vh.insert(std::make_pair(node_id,vh)); + } + + std::vector > constrained_edges; + + // insert constraints that are interior to the triangle (in the case + // no edges are collinear in the meshes) + insert_constrained_edges(node_ids,cdt,id_to_CDT_vh,constrained_edges); + + // insert constraints between points that are on the boundary + // (not a contrained on the triangle boundary) + if (it_fb!=face_boundaries.end()) //is f not a triangle ? + { + for (int i=0;i<3;++i) + { + Node_ids& ids=it_fb->second.node_ids_array[i]; + insert_constrained_edges(ids,cdt,id_to_CDT_vh,constrained_edges,1); + } + } + + //insert coplanar edges for endpoints of triangles + for (int i=0;i<3;++i){ + Node_id nindex=triangle_vertices[i]->info(); + if ( nindex < nb_nodes ) + insert_constrained_edges_coplanar_case(nindex,cdt,id_to_CDT_vh); + } + + //XSL_TAG_CPL_VERT + //collect edges incident to a point that is the intersection of two + // coplanar faces. This ensure that triangulations are compatible. + if (it_fb!=face_boundaries.end()) //is f not a triangle ? + { + for (typename CDT::Finite_vertices_iterator + vit=cdt.finite_vertices_begin(), + vit_end=cdt.finite_vertices_end();vit_end!=vit;++vit) + { + //skip original vertices (that are not nodes) and non-coplanar face + // issued vertices (this is working because intersection points + // between coplanar facets are the first inserted) + if (vit->info() >= nb_nodes || + vit->info() >= number_coplanar_vertices) continue; + // \todo no need to insert constrained edges (they also are constrained + // in the other mesh)!! + typename std::map< Node_id,std::set >::iterator res = + coplanar_constraints.insert( + std::make_pair(vit->info(),std::set())).first; + //turn around the vertex and get incident edge + typename CDT::Edge_circulator start=cdt.incident_edges(vit); + typename CDT::Edge_circulator curr=start; + do{ + if (cdt.is_infinite(*curr) ) continue; + typename CDT::Edge mirror=cdt.mirror_edge(*curr); + if ( cdt.is_infinite( curr->first->vertex(curr->second) ) || + cdt.is_infinite( mirror.first->vertex(mirror.second) ) ) + continue; // skip edges that are on the boundary of the triangle + // (these are already constrained) + //insert edges in the set of constraints + CDT_Vertex_handle vh=vit; + int nindex = curr->first->vertex((curr->second+1)%3)==vh + ? (curr->second+2)%3 + : (curr->second+1)%3; + CDT_Vertex_handle vn=curr->first->vertex(nindex); + if ( vit->info() > vn->info() || vn->info()>=nb_nodes) + continue; //take only one out of the two edges + skip input + CGAL_assertion(vn->info()second.insert( vn->info() ); + }while(start!=++curr); + } + } + + // import the triangle in `cdt` in the face `f` of `tm` + triangulate_a_face(f, tm, nodes, node_ids, node_id_to_vertex, + edge_to_hedge, cdt, vpm, output_builder, user_visitor); + + // TODO Here we do the update only for internal edges. + // Update for border halfedges could be done during the split + + //3) mark halfedges that are common to two polyhedral surfaces + //recover halfedges inserted that are on the intersection + typedef std::pair Node_id_pair; + for(const Node_id_pair& node_id_pair : constrained_edges) + { + typename std::map + ::iterator it_poly_hedge=edge_to_hedge.find(node_id_pair); + //we cannot have an assertion here in case an edge or part of an edge is a constraints. + //Indeed, the graph_of_constraints report an edge 0,1 and 1,0 for example while only one of the two + //is defined as one of them defines an adjacent face + //CGAL_assertion(it_poly_hedge!=edge_to_hedge.end()); + if( it_poly_hedge!=edge_to_hedge.end() ){ + call_put(marks_on_edges,tm,edge(it_poly_hedge->second,tm),true); + output_builder.set_edge_per_polyline(tm,node_id_pair,it_poly_hedge->second); + } + else{ + //WARNING: in few case this is needed if the marked edge is on the border + //to optimize it might be better to only use sorted pair. TAG_SLXX1 + Node_id_pair opposite_pair(node_id_pair.second,node_id_pair.first); + it_poly_hedge=edge_to_hedge.find(opposite_pair); + CGAL_assertion( it_poly_hedge!=edge_to_hedge.end() ); + + call_put(marks_on_edges,tm,edge(it_poly_hedge->second,tm),true); + output_builder.set_edge_per_polyline(tm,opposite_pair,it_poly_hedge->second); + } + } + } + } + + void finalize(INodes& nodes, + const TriangleMesh& tm1, + const TriangleMesh& tm2, + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2) { nodes.all_nodes_created(); TriangleMesh* tm1_ptr = const_cast(&tm1); TriangleMesh* tm2_ptr = const_cast(&tm2); - std::map vpms; - vpms[tm1_ptr] = vpm1; - vpms[tm2_ptr] = vpm2; - vertex_descriptor null_vertex = Graph_traits::null_vertex(); const Node_id nb_nodes = nodes.size(); // we reserve nb_nodes+3 because we use the last three entries for the @@ -610,7 +954,6 @@ public: //store for each triangle face which boundary is intersected by the other surface, //original vertices (and halfedges in the refined mesh pointing on these vertices) - typedef boost::unordered_map Face_boundaries; std::map mesh_to_face_boundaries; //0) For each triangle mesh, collect original vertices that belongs to the intersection. @@ -692,94 +1035,10 @@ public: for (typename std::map::iterator it=on_edge.begin(); it!=on_edge.end(); ++it) { - TriangleMesh& tm=*it->first; - const VertexPointMap& vpm=vpms[&tm]; - On_edge_map& on_edge_map=it->second; - On_face_map& on_face_map=on_face[&tm]; - Face_boundaries& face_boundaries=mesh_to_face_boundaries[&tm]; - - for(typename On_edge_map::iterator it2=on_edge_map.begin(); - it2!=on_edge_map.end(); - ++it2) - { - //the edge to be split - halfedge_descriptor hedge=halfedge(it2->first,tm); - //indices of the nodes to be inserted - Node_ids& node_ids=it2->second; - CGAL_assertion( std::set(node_ids.begin(), node_ids.end()) - .size()==node_ids.size() ); - //sort nodes along the egde to allow consecutive splits - sort_vertices_along_hedge(node_ids,hedge,tm,vpm,nodes); - - //save original face and nodes for face of hedge (1) - if ( !is_border(hedge,tm) ){ - face_descriptor f=face(hedge,tm); - typename Face_boundaries::iterator it_face = face_boundaries.find(f); - if (it_face==face_boundaries.end()) - it_face=face_boundaries.insert(std::make_pair(f,Face_boundary(hedge,tm))).first; - it_face->second.copy_node_ids(hedge,node_ids.begin(),node_ids.end()); - } - - //save original face and nodes for face of hedge->opposite (2) - typename Face_boundaries::iterator opposite_original_info=face_boundaries.end(); - halfedge_descriptor hedge_opp = opposite(hedge,tm); - if ( !is_border(hedge_opp,tm) ){ - face_descriptor f=face(hedge_opp,tm); - opposite_original_info=face_boundaries.find(f); - if (opposite_original_info==face_boundaries.end()) - opposite_original_info=face_boundaries.insert(std::make_pair(f,Face_boundary(hedge_opp,tm))).first; - opposite_original_info->second.copy_node_ids(hedge_opp,node_ids.rbegin(),node_ids.rend()); - } - - typename Mesh_to_map_node::iterator it_map=mesh_to_node_id_to_vertex.find(&tm); - CGAL_assertion(it_map!=mesh_to_node_id_to_vertex.end()); - //a map to identify the vertex in the polyhedron corresponding to an intersection point - Node_id_to_vertex& node_id_to_vertex=it_map->second; - - CGAL_assertion_code(vertex_descriptor original_vertex=source(hedge,tm);) - - //We need an edge incident to the source vertex of hedge. This is the first opposite edge created. - bool first=true; - halfedge_descriptor hedge_incident_to_src=Graph_traits::null_halfedge(); - bool hedge_is_marked = call_get(marks_on_edges,tm,edge(hedge,tm)); - //do split the edges - CGAL_assertion_code(vertex_descriptor expected_src=source(hedge,tm)); - for(std::size_t node_id : node_ids) - { - halfedge_descriptor hnew = Euler::split_edge(hedge, tm); - CGAL_assertion(expected_src==source(hnew,tm)); - vertex_descriptor vnew=target(hnew,tm); -// user_visitor.new_vertex_added(node_id, vnew, tm); // NODE_VISITOR_TAG - nodes.call_put(vpm, vnew, node_id, tm); - // register the new vertex in the output builder - output_builder.set_vertex_id(vnew, node_id, tm); - node_id_to_vertex[node_id]=vnew; - if (first){ - first=false; - hedge_incident_to_src=next(opposite(hedge,tm),tm); - } - - //update marker tags. If the edge was marked, then the resulting edges in the split must be marked - if ( hedge_is_marked ) - call_put(marks_on_edges,tm,edge(hnew,tm),true); - - CGAL_assertion_code(expected_src=vnew); - } - - CGAL_assertion(target(hedge_incident_to_src,tm)==original_vertex); - CGAL_assertion(face(hedge_incident_to_src,tm)==face(hedge_opp,tm)); - - //save original face and nodes for face of hedge->opposite (2) - if ( !is_border(hedge_opp,tm) ){ - CGAL_assertion(opposite_original_info!=face_boundaries.end()); - opposite_original_info->second.update_original_halfedge( - hedge_opp,hedge_incident_to_src,tm); - } - - //insert the two incident faces in on_face map so that they will be triangulated. - if (!is_border(hedge,tm)) on_face_map[face(hedge,tm)]; - if (!is_border(hedge_opp,tm)) on_face_map[face(hedge_opp,tm)]; - } + if(it->first == tm1_ptr) + split_halfedges(it, vpm1, nodes, mesh_to_face_boundaries); + else + split_halfedges(it, vpm2, nodes, mesh_to_face_boundaries); } //2)triangulation of the triangle faces containing intersection point in their interior @@ -787,247 +1046,10 @@ public: for (typename std::map::iterator it=on_face.begin(); it!=on_face.end(); ++it) { - TriangleMesh& tm=*it->first; - const VertexPointMap& vpm=vpms[&tm]; - On_face_map& on_face_map=it->second; - Face_boundaries& face_boundaries=mesh_to_face_boundaries[&tm]; - Node_id_to_vertex& node_id_to_vertex=mesh_to_node_id_to_vertex[&tm]; - Vertex_to_node_id& vertex_to_node_id=mesh_to_vertex_to_node_id[&tm]; - - for (typename On_face_map::iterator it=on_face_map.begin(); - it!=on_face_map.end();++it) - { - face_descriptor f = it->first; //the face to be triangulated - Node_ids& node_ids = it->second; // ids of nodes in the interior of f - typename Face_boundaries::iterator it_fb=face_boundaries.find(f); - - std::map id_to_CDT_vh; - - //associate an edge of the triangulation to a halfedge in a given polyhedron - std::map,halfedge_descriptor> edge_to_hedge; - - // the vertices of f - std::array f_vertices; - // the node_id of an input vertex or a fake id (>=nb_nodes) - std::array f_indices = {{nb_nodes,nb_nodes+1,nb_nodes+2}}; - if (it_fb!=face_boundaries.end()){ //the boundary of the triangle face was refined - f_vertices[0]=it_fb->second.vertices[0]; - f_vertices[1]=it_fb->second.vertices[1]; - f_vertices[2]=it_fb->second.vertices[2]; - update_face_indices(f_vertices,f_indices,vertex_to_node_id); - if (doing_autorefinement) - it_fb->second.update_node_id_to_vertex_map(node_id_to_vertex, tm); - } - else{ - CGAL_assertion( is_triangle(halfedge(f,tm),tm) ); - halfedge_descriptor h0=halfedge(f,tm), h1=next(h0,tm), h2=next(h1,tm); - f_vertices[0]=target(h0,tm); //nb_nodes - f_vertices[1]=target(h1,tm); //nb_nodes+1 - f_vertices[2]=target(h2,tm); //nb_nodes+2 - - update_face_indices(f_vertices,f_indices,vertex_to_node_id); - edge_to_hedge[std::make_pair( f_indices[2],f_indices[0] )] = h0; - edge_to_hedge[std::make_pair( f_indices[0],f_indices[1] )] = h1; - edge_to_hedge[std::make_pair( f_indices[1],f_indices[2] )] = h2; - } - - typename EK::Point_3 p = nodes.to_exact(get(vpm,f_vertices[0])), - q = nodes.to_exact(get(vpm,f_vertices[1])), - r = nodes.to_exact(get(vpm,f_vertices[2])); -///TODO use a positive normal and remove all work around to guarantee that triangulation of coplanar patches are compatible - CDT_traits traits(typename EK::Construct_normal_3()(p,q,r)); - CDT cdt(traits); - - // insert triangle points - std::array triangle_vertices; - //we can do this to_exact because these are supposed to be input points. - triangle_vertices[0]=cdt.insert_outside_affine_hull(p); - triangle_vertices[1]=cdt.insert_outside_affine_hull(q); - triangle_vertices[2]=cdt.tds().insert_dim_up(cdt.infinite_vertex(), false); - triangle_vertices[2]->set_point(r); - - - triangle_vertices[0]->info()=f_indices[0]; - triangle_vertices[1]->info()=f_indices[1]; - triangle_vertices[2]->info()=f_indices[2]; - - node_id_to_vertex[nb_nodes ]=f_vertices[0]; - node_id_to_vertex[nb_nodes+1]=f_vertices[1]; - node_id_to_vertex[nb_nodes+2]=f_vertices[2]; - - //if one of the triangle input vertex is also a node - for (int ik=0;ik<3;++ik){ - if ( f_indices[ik]vertex(oi) ) ); - } - - // In this loop, for each original edge of the triangle, we insert - // the constrained edges and we recover the halfedge_descriptor - // corresponding to these constrained (they are already in tm) - Face_boundary& f_boundary=it_fb->second; - for (int i=0;i<3;++i){ - //handle case of halfedge starting at triangle_vertices[i] - // and ending at triangle_vertices[(i+1)%3] - - const Node_ids& ids_on_edge=f_boundary.node_ids_array[i]; - CDT_Vertex_handle previous=triangle_vertices[i]; - Node_id prev_index=f_indices[i];// node-id of the mesh vertex - halfedge_descriptor hedge = next(f_boundary.halfedges[(i+2)%3],tm); - CGAL_assertion( source(hedge,tm)==f_boundary.vertices[i] ); - if (!ids_on_edge.empty()){ //is there at least one node on this edge? - // fh must be an infinite face - // The points must be ordered from fh->vertex(cw(infinite_vertex)) to fh->vertex(ccw(infinite_vertex)) - for(Node_id id : ids_on_edge) - { - CDT_Vertex_handle vh=insert_point_on_ch_edge(cdt,infinite_faces[i],nodes.exact_node(id)); - vh->info()=id; - id_to_CDT_vh.insert(std::make_pair(id,vh)); - edge_to_hedge[std::make_pair(prev_index,id)]=hedge; - previous=vh; - hedge=next(hedge,tm); - prev_index=id; - } - } - else{ - CGAL_assertion_code(halfedge_descriptor hd=f_boundary.halfedges[i]); - CGAL_assertion( target(hd,tm) == f_boundary.vertices[(i+1)%3] ); - CGAL_assertion( source(hd,tm) == f_boundary.vertices[ i ] ); - } - CGAL_assertion(hedge==f_boundary.halfedges[i]); - edge_to_hedge[std::make_pair(prev_index,f_indices[(i+1)%3])] = - it_fb->second.halfedges[i]; - } - } - - //insert point inside face - for(Node_id node_id : node_ids) - { - CDT_Vertex_handle vh=cdt.insert(nodes.exact_node(node_id)); - vh->info()=node_id; - id_to_CDT_vh.insert(std::make_pair(node_id,vh)); - } - - std::vector > constrained_edges; - - // insert constraints that are interior to the triangle (in the case - // no edges are collinear in the meshes) - insert_constrained_edges(node_ids,cdt,id_to_CDT_vh,constrained_edges); - - // insert constraints between points that are on the boundary - // (not a contrained on the triangle boundary) - if (it_fb!=face_boundaries.end()) //is f not a triangle ? - { - for (int i=0;i<3;++i) - { - Node_ids& ids=it_fb->second.node_ids_array[i]; - insert_constrained_edges(ids,cdt,id_to_CDT_vh,constrained_edges,1); - } - } - - //insert coplanar edges for endpoints of triangles - for (int i=0;i<3;++i){ - Node_id nindex=triangle_vertices[i]->info(); - if ( nindex < nb_nodes ) - insert_constrained_edges_coplanar_case(nindex,cdt,id_to_CDT_vh); - } - - //XSL_TAG_CPL_VERT - //collect edges incident to a point that is the intersection of two - // coplanar faces. This ensure that triangulations are compatible. - if (it_fb!=face_boundaries.end()) //is f not a triangle ? - { - for (typename CDT::Finite_vertices_iterator - vit=cdt.finite_vertices_begin(), - vit_end=cdt.finite_vertices_end();vit_end!=vit;++vit) - { - //skip original vertices (that are not nodes) and non-coplanar face - // issued vertices (this is working because intersection points - // between coplanar facets are the first inserted) - if (vit->info() >= nb_nodes || - vit->info() >= number_coplanar_vertices) continue; - // \todo no need to insert constrained edges (they also are constrained - // in the other mesh)!! - typename std::map< Node_id,std::set >::iterator res = - coplanar_constraints.insert( - std::make_pair(vit->info(),std::set())).first; - //turn around the vertex and get incident edge - typename CDT::Edge_circulator start=cdt.incident_edges(vit); - typename CDT::Edge_circulator curr=start; - do{ - if (cdt.is_infinite(*curr) ) continue; - typename CDT::Edge mirror=cdt.mirror_edge(*curr); - if ( cdt.is_infinite( curr->first->vertex(curr->second) ) || - cdt.is_infinite( mirror.first->vertex(mirror.second) ) ) - continue; // skip edges that are on the boundary of the triangle - // (these are already constrained) - //insert edges in the set of constraints - CDT_Vertex_handle vh=vit; - int nindex = curr->first->vertex((curr->second+1)%3)==vh - ? (curr->second+2)%3 - : (curr->second+1)%3; - CDT_Vertex_handle vn=curr->first->vertex(nindex); - if ( vit->info() > vn->info() || vn->info()>=nb_nodes) - continue; //take only one out of the two edges + skip input - CGAL_assertion(vn->info()second.insert( vn->info() ); - }while(start!=++curr); - } - } - - // import the triangle in `cdt` in the face `f` of `tm` - triangulate_a_face(f, tm, nodes, node_ids, node_id_to_vertex, - edge_to_hedge, cdt, vpm, output_builder, user_visitor); - - // TODO Here we do the update only for internal edges. - // Update for border halfedges could be done during the split - - //3) mark halfedges that are common to two polyhedral surfaces - //recover halfedges inserted that are on the intersection - typedef std::pair Node_id_pair; - for(const Node_id_pair& node_id_pair : constrained_edges) - { - typename std::map - ::iterator it_poly_hedge=edge_to_hedge.find(node_id_pair); - //we cannot have an assertion here in case an edge or part of an edge is a constraints. - //Indeed, the graph_of_constraints report an edge 0,1 and 1,0 for example while only one of the two - //is defined as one of them defines an adjacent face - //CGAL_assertion(it_poly_hedge!=edge_to_hedge.end()); - if( it_poly_hedge!=edge_to_hedge.end() ){ - call_put(marks_on_edges,tm,edge(it_poly_hedge->second,tm),true); - output_builder.set_edge_per_polyline(tm,node_id_pair,it_poly_hedge->second); - } - else{ - //WARNING: in few case this is needed if the marked edge is on the border - //to optimize it might be better to only use sorted pair. TAG_SLXX1 - Node_id_pair opposite_pair(node_id_pair.second,node_id_pair.first); - it_poly_hedge=edge_to_hedge.find(opposite_pair); - CGAL_assertion( it_poly_hedge!=edge_to_hedge.end() ); - - call_put(marks_on_edges,tm,edge(it_poly_hedge->second,tm),true); - output_builder.set_edge_per_polyline(tm,opposite_pair,it_poly_hedge->second); - } - } - } + if(it->first == tm1_ptr) + triangulate_intersected_faces(it, vpm1, nodes, mesh_to_face_boundaries); + else + triangulate_intersected_faces(it, vpm2, nodes, mesh_to_face_boundaries); } nodes.finalize(); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h index c188231238c..4103dcd41ae 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h @@ -583,7 +583,8 @@ next_marked_halfedge_around_target_vertex( template void import_polyline( @@ -598,8 +599,8 @@ void import_polyline( VertexMap& pm1_to_output_vertices, const IntersectionEdgeMap& intersection_edges1, const IntersectionEdgeMap& intersection_edges2, - const VertexPointMap& vpm1, - const VertexPointMap& /*vpm2*/, + const VertexPointMap1& vpm1, + const VertexPointMap2& /*vpm2*/, const VertexPointMapOut& vpm_out, std::vector ::edge_descriptor>& output_shared_edges) @@ -1004,7 +1005,8 @@ void append_patches_to_triangle_mesh( template < class TriangleMesh, class IntersectionEdgeMap, - class VertexPointMap, + class VertexPointMap1, + class VertexPointMap2, class VertexPointMapOut, class EdgeMarkMap1, class EdgeMarkMap2, @@ -1024,8 +1026,8 @@ void fill_new_triangle_mesh( const IntersectionPolylines& polylines, const IntersectionEdgeMap& intersection_edges1, const IntersectionEdgeMap& intersection_edges2, - const VertexPointMap& vpm1, - const VertexPointMap& vpm2, + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2, const VertexPointMapOut& vpm_out, const EdgeMarkMap1& edge_mark_map1, const EdgeMarkMap2& edge_mark_map2, @@ -1261,7 +1263,8 @@ template +template std::tuple::halfedge_descriptor, bool,bool> @@ -87,16 +87,21 @@ intersection_type( typename boost::graph_traits::face_descriptor f_2, const TriangleMesh& tm1, const TriangleMesh& tm2, - const VertexPointMap& vpm1, - const VertexPointMap& vpm2) + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2) { typedef boost::graph_traits GT; typedef typename GT::halfedge_descriptor halfedge_descriptor; typedef std::tuple result_type; - typedef typename boost::property_traits::reference Point_ref; - typedef typename boost::property_traits::value_type Point_3; + typedef typename boost::property_traits::reference Point_ref; + typedef typename boost::property_traits::value_type Point_3; typedef typename Kernel_traits::Kernel Kernel; + CGAL_static_assertion((std::is_same::reference, + typename boost::property_traits::reference>::value)); + CGAL_static_assertion((std::is_same::value_type, + typename boost::property_traits::value_type>::value)); + halfedge_descriptor h_2=halfedge(f_2,tm2); Point_ref a = get(vpm2, target(h_2,tm2) ); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h index d5389e9f482..d7db498ae8c 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h @@ -67,15 +67,15 @@ public: }; template class Collect_face_bbox_per_edge_bbox_with_coplanar_handling { protected: const TriangleMesh& tm_faces; const TriangleMesh& tm_edges; - const VertexPointMap& vpmap_tmf; - const VertexPointMap& vpmap_tme; + const VertexPointMapF& vpmap_tmf; + const VertexPointMapE& vpmap_tme; EdgeToFaces& edge_to_faces; CoplanarFaceSet& coplanar_faces; @@ -83,7 +83,7 @@ protected: typedef typename Graph_traits::face_descriptor face_descriptor; typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::property_traits::reference Point; + typedef typename boost::property_traits::reference Point; typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy; typedef CGAL::Box_intersection_d::Box_with_info_d Box; @@ -92,8 +92,8 @@ public: Collect_face_bbox_per_edge_bbox_with_coplanar_handling( const TriangleMesh& tm_faces, const TriangleMesh& tm_edges, - const VertexPointMap& vpmap_tmf, - const VertexPointMap& vpmap_tme, + const VertexPointMapF& vpmap_tmf, + const VertexPointMapE& vpmap_tme, EdgeToFaces& edge_to_faces, CoplanarFaceSet& coplanar_faces) : tm_faces(tm_faces) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h index 72c94c605d7..12035e0b01b 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h @@ -83,10 +83,10 @@ struct Default_surface_intersection_visitor{ void start_new_polyline(std::size_t,std::size_t){} void add_node_to_polyline(std::size_t){} void input_have_coplanar_faces(){} - template + template void finalize(T&, const TriangleMesh&, const TriangleMesh&, - const VertexPointMap&, const VertexPointMap&) + const VPM1, const VPM2) {} void new_node_added_triple_face(std::size_t /* node_id */, face_descriptor /* f1 */, @@ -144,7 +144,7 @@ struct Node_id_set { }; template< class TriangleMesh, - class VertexPointMap, + class VertexPointMap1, class VertexPointMap2, class Node_visitor=Default_surface_intersection_visitor > class Intersection_of_triangle_meshes @@ -175,7 +175,7 @@ class Intersection_of_triangle_meshes // may contain several segments. typedef std::map< Face_pair_and_int, Node_id_set > Faces_to_nodes_map; typedef Intersection_nodes Node_vector; // data members @@ -188,11 +188,13 @@ class Intersection_of_triangle_meshes Faces_to_nodes_map f_to_node; //Associate a pair of triangles to their intersection points std::vector extra_terminal_nodes; //used only for autorefinement CGAL_assertion_code(bool doing_autorefinement;) + // member functions + template void filter_intersections(const TriangleMesh& tm_f, const TriangleMesh& tm_e, - const VertexPointMap& vpm_f, - const VertexPointMap& vpm_e, + const VPMF& vpm_f, + const VPME& vpm_e, bool throw_on_self_intersection) { std::vector face_boxes, edge_boxes; @@ -237,7 +239,7 @@ class Intersection_of_triangle_meshes Callback callback(tm_f, tm_e, edge_to_faces); #else typedef Collect_face_bbox_per_edge_bbox_with_coplanar_handling< - TriangleMesh, VertexPointMap, Edge_to_faces, Coplanar_face_set> + TriangleMesh, VPMF, VPME, Edge_to_faces, Coplanar_face_set> Callback; Callback callback(tm_f, tm_e, vpm_f, vpm_e, edge_to_faces, coplanar_faces); #endif @@ -258,8 +260,9 @@ class Intersection_of_triangle_meshes } // for autorefinement + template void filter_intersections(const TriangleMesh& tm, - const VertexPointMap& vpm) + const VPM& vpm) { std::vector face_boxes, edge_boxes; std::vector face_boxes_ptr, edge_boxes_ptr; @@ -296,7 +299,7 @@ class Intersection_of_triangle_meshes Edge_to_faces& edge_to_faces = stm_edge_to_ltm_faces; typedef Collect_face_bbox_per_edge_bbox_with_coplanar_handling_one_mesh< - TriangleMesh, VertexPointMap, Edge_to_faces, Coplanar_face_set> + TriangleMesh, VPM, Edge_to_faces, Coplanar_face_set> Callback; Callback callback(tm, vpm, edge_to_faces, coplanar_faces); @@ -527,14 +530,14 @@ class Intersection_of_triangle_meshes } } - void compute_intersection_of_coplanar_faces( - Node_id& current_node, - const TriangleMesh& tm1, - const TriangleMesh& tm2, - const VertexPointMap& vpm1, - const VertexPointMap& vpm2) + template + void compute_intersection_of_coplanar_faces(Node_id& current_node, + const TriangleMesh& tm_f, + const TriangleMesh& tm_e, + const VPMF& vpm_f, + const VPME& vpm_e) { - CGAL_assertion( &tm1 < &tm2 || &tm1==&tm2 ); + CGAL_assertion(&tm_f < &tm_e || &tm_f == &tm_e); typedef std::tuple Cpl_inter_pt; std::list inter_pts; // compute the intersection points between the two coplanar faces - intersection_coplanar_faces(f1, f2, tm1, tm2, vpm1, vpm2, inter_pts); + intersection_coplanar_faces(f1, f2, tm_f, tm_e, vpm_f, vpm_e, inter_pts); std::size_t nb_pts=inter_pts.size(); std::vector cpln_nodes; cpln_nodes.reserve(nb_pts); @@ -566,7 +569,7 @@ class Intersection_of_triangle_meshes Node_id node_id; bool is_new_node; std::tie(node_id, is_new_node) = - get_or_create_node(ipt,current_node,coplanar_node_map,tm1,tm2); + get_or_create_node(ipt,current_node,coplanar_node_map,tm_f,tm_e); cpln_nodes.push_back(node_id); switch(ipt.type_1){ @@ -574,13 +577,13 @@ class Intersection_of_triangle_meshes { switch(ipt.type_2){ case ON_VERTEX: - handle_coplanar_case_VERTEX_VERTEX(ipt.info_1,ipt.info_2,tm1,tm2,node_id,is_new_node); + handle_coplanar_case_VERTEX_VERTEX(ipt.info_1,ipt.info_2,tm_f,tm_e,node_id,is_new_node); break; case ON_EDGE: - handle_coplanar_case_VERTEX_EDGE(ipt.info_1,ipt.info_2,tm1,tm2,node_id,is_new_node); + handle_coplanar_case_VERTEX_EDGE(ipt.info_1,ipt.info_2,tm_f,tm_e,node_id,is_new_node); break; case ON_FACE: - handle_coplanar_case_VERTEX_FACE(ipt.info_1,ipt.info_2,tm1,tm2,node_id,is_new_node); + handle_coplanar_case_VERTEX_FACE(ipt.info_1,ipt.info_2,tm_f,tm_e,node_id,is_new_node); break; default: CGAL_error_msg("Should not get there!"); } @@ -590,15 +593,15 @@ class Intersection_of_triangle_meshes { switch(ipt.type_2){ case ON_VERTEX: - handle_coplanar_case_VERTEX_EDGE(ipt.info_2,ipt.info_1,tm2,tm1,node_id,is_new_node); + handle_coplanar_case_VERTEX_EDGE(ipt.info_2,ipt.info_1,tm_e,tm_f,node_id,is_new_node); break; case ON_EDGE: { if(is_new_node) - visitor.new_node_added(node_id,ON_EDGE,ipt.info_1,ipt.info_2,tm1,tm2,false,false); - typename Edge_to_faces::iterator it_ets=stm_edge_to_ltm_faces.find(edge(ipt.info_1,tm1)); + visitor.new_node_added(node_id,ON_EDGE,ipt.info_1,ipt.info_2,tm_f,tm_e,false,false); + typename Edge_to_faces::iterator it_ets=stm_edge_to_ltm_faces.find(edge(ipt.info_1,tm_f)); Face_set* fset = (it_ets!=stm_edge_to_ltm_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_edge(node_id,fset,ipt.info_1,ipt.info_2,tm1,tm2); + cip_handle_case_edge(node_id,fset,ipt.info_1,ipt.info_2,tm_f,tm_e); } break; default: CGAL_error_msg("Should not get there!"); @@ -608,7 +611,7 @@ class Intersection_of_triangle_meshes case ON_FACE: { CGAL_assertion(ipt.type_2==ON_VERTEX); - handle_coplanar_case_VERTEX_FACE(ipt.info_2,ipt.info_1,tm2,tm1,node_id,is_new_node); + handle_coplanar_case_VERTEX_FACE(ipt.info_2,ipt.info_1,tm_e,tm_f,node_id,is_new_node); } break; default: CGAL_error_msg("Should not get there!"); @@ -637,64 +640,66 @@ class Intersection_of_triangle_meshes //add a new node in the final graph. //it is the intersection of the triangle with the segment - void add_new_node(halfedge_descriptor h_1, - face_descriptor f_2, - const TriangleMesh& tm1, - const TriangleMesh& tm2, - const VertexPointMap& vpm1, - const VertexPointMap& vpm2, + template + void add_new_node(halfedge_descriptor h_f, + face_descriptor f_e, + const TriangleMesh& tm_f, + const TriangleMesh& tm_e, + const VPMF& vpm_f, + const VPME& vpm_e, std::tuple inter_res) { if ( std::get<3>(inter_res) ) // is edge target in triangle plane - nodes.add_new_node(get(vpm1, target(h_1,tm1))); + nodes.add_new_node(get(vpm_f, target(h_f, tm_f))); else{ if (std::get<2>(inter_res)) // is edge source in triangle plane - nodes.add_new_node(get(vpm1, source(h_1,tm1))); + nodes.add_new_node(get(vpm_f, source(h_f, tm_f))); else - nodes.add_new_node(h_1,f_2,tm1,tm2,vpm1,vpm2); + nodes.add_new_node(h_f, f_e, tm_f, tm_e, vpm_f, vpm_e); } } - void compute_intersection_points(Edge_to_faces& tm1_edge_to_tm2_faces, - const TriangleMesh& tm1, - const TriangleMesh& tm2, - const VertexPointMap& vpm1, - const VertexPointMap& vpm2, + template + void compute_intersection_points(Edge_to_faces& tmf_edge_to_tme_faces, + const TriangleMesh& tm_f, + const TriangleMesh& tm_e, + const VPMF& vpm_f, + const VPME& vpm_e, Node_id& current_node) { typedef std::tuple Inter_type; - for(typename Edge_to_faces::iterator it=tm1_edge_to_tm2_faces.begin(); - it!=tm1_edge_to_tm2_faces.end();++it) + for(typename Edge_to_faces::iterator it=tmf_edge_to_tme_faces.begin(); + it!=tmf_edge_to_tme_faces.end();++it) { - edge_descriptor e_1=it->first; - halfedge_descriptor h_1=halfedge(e_1,tm1); + edge_descriptor e_f=it->first; + halfedge_descriptor h_f=halfedge(e_f,tm_f); Face_set& fset=it->second; while (!fset.empty()){ - face_descriptor f_2=*fset.begin(); + face_descriptor f_e=*fset.begin(); - Inter_type res=intersection_type(h_1,f_2,tm1,tm2,vpm1,vpm2); + Inter_type res=intersection_type(h_f,f_e,tm_f,tm_e,vpm_f,vpm_e); Intersection_type type=std::get<0>(res); - //handle degenerate case: one extremity of edge belong to f_2 + //handle degenerate case: one extremity of edge belong to f_e std::vector all_edges; if ( std::get<3>(res) ) // is edge target in triangle plane - std::copy(halfedges_around_target(h_1,tm1).first, - halfedges_around_target(h_1,tm1).second, + std::copy(halfedges_around_target(h_f,tm_f).first, + halfedges_around_target(h_f,tm_f).second, std::back_inserter(all_edges)); else{ if ( std::get<2>(res) ) // is edge source in triangle plane - std::copy(halfedges_around_source(h_1,tm1).first, - halfedges_around_source(h_1,tm1).second, + std::copy(halfedges_around_source(h_f,tm_f).first, + halfedges_around_source(h_f,tm_f).second, std::back_inserter(all_edges)); else - all_edges.push_back(h_1); + all_edges.push_back(h_f); } - CGAL_precondition(all_edges[0]==h_1 || all_edges[0]==opposite(h_1,tm1)); + CGAL_precondition(all_edges[0]==h_f || all_edges[0]==opposite(h_f,tm_f)); // #ifdef USE_DETECTION_MULTIPLE_DEFINED_EDGES // check_coplanar_edges(std::next(all_edges.begin()), @@ -717,20 +722,20 @@ class Intersection_of_triangle_meshes // Case when the edge pierces the face in its interior. case ON_FACE: { - CGAL_assertion(f_2==face(std::get<1>(res),tm2)); + CGAL_assertion(f_e==face(std::get<1>(res),tm_e)); Node_id node_id=++current_node; - add_new_node(h_1,f_2,tm1,tm2,vpm1,vpm2,res); - visitor.new_node_added(node_id,ON_FACE,h_1,halfedge(f_2,tm2),tm1,tm2,std::get<3>(res),std::get<2>(res)); + add_new_node(h_f,f_e,tm_f,tm_e,vpm_f,vpm_e,res); + visitor.new_node_added(node_id,ON_FACE,h_f,halfedge(f_e,tm_e),tm_f,tm_e,std::get<3>(res),std::get<2>(res)); for (;it_edge!=all_edges.end();++it_edge){ - add_intersection_point_to_face_and_all_edge_incident_faces(f_2,*it_edge,tm2,tm1,node_id); + add_intersection_point_to_face_and_all_edge_incident_faces(f_e,*it_edge,tm_e,tm_f,node_id); //erase face from the list to test intersection with it_edge if ( it_edge==all_edges.begin() ) fset.erase(fset.begin()); else { - typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(*it_edge,tm1)); - if(it_ets!=tm1_edge_to_tm2_faces.end()) it_ets->second.erase(f_2); + typename Edge_to_faces::iterator it_ets=tmf_edge_to_tme_faces.find(edge(*it_edge,tm_f)); + if(it_ets!=tmf_edge_to_tme_faces.end()) it_ets->second.erase(f_e); } } } // end case ON_FACE @@ -740,17 +745,17 @@ class Intersection_of_triangle_meshes case ON_EDGE: { Node_id node_id=++current_node; - add_new_node(h_1,f_2,tm1,tm2,vpm1,vpm2,res); - halfedge_descriptor h_2=std::get<1>(res); - visitor.new_node_added(node_id,ON_EDGE,h_1,h_2,tm1,tm2,std::get<3>(res),std::get<2>(res)); + add_new_node(h_f,f_e,tm_f,tm_e,vpm_f,vpm_e,res); + halfedge_descriptor h_e=std::get<1>(res); + visitor.new_node_added(node_id,ON_EDGE,h_f,h_e,tm_f,tm_e,std::get<3>(res),std::get<2>(res)); for (;it_edge!=all_edges.end();++it_edge){ if ( it_edge!=all_edges.begin() ){ - typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(*it_edge,tm1)); - Face_set* fset_bis = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_edge(node_id,fset_bis,*it_edge,h_2,tm1,tm2); + typename Edge_to_faces::iterator it_ets=tmf_edge_to_tme_faces.find(edge(*it_edge,tm_f)); + Face_set* fset_bis = (it_ets!=tmf_edge_to_tme_faces.end())?&(it_ets->second):nullptr; + cip_handle_case_edge(node_id,fset_bis,*it_edge,h_e,tm_f,tm_e); } else - cip_handle_case_edge(node_id,&fset,*it_edge,h_2,tm1,tm2); + cip_handle_case_edge(node_id,&fset,*it_edge,h_e,tm_f,tm_e); } } // end case ON_EDGE break; @@ -758,18 +763,18 @@ class Intersection_of_triangle_meshes case ON_VERTEX: { Node_id node_id=++current_node; - halfedge_descriptor h_2=std::get<1>(res); - nodes.add_new_node(get(vpm2, target(h_2,tm2))); //we use the original vertex to create the node + halfedge_descriptor h_e=std::get<1>(res); + nodes.add_new_node(get(vpm_e, target(h_e,tm_e))); //we use the original vertex to create the node //before it was ON_FACE but do not remember why, probably a bug... - visitor.new_node_added(node_id,ON_VERTEX,h_1,h_2,tm1,tm2,std::get<3>(res),std::get<2>(res)); + visitor.new_node_added(node_id,ON_VERTEX,h_f,h_e,tm_f,tm_e,std::get<3>(res),std::get<2>(res)); for (;it_edge!=all_edges.end();++it_edge){ if ( it_edge!=all_edges.begin() ){ - typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(*it_edge,tm1)); - Face_set* fset_bis = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_vertex(node_id,fset_bis,*it_edge,h_2,tm1,tm2); + typename Edge_to_faces::iterator it_ets=tmf_edge_to_tme_faces.find(edge(*it_edge,tm_f)); + Face_set* fset_bis = (it_ets!=tmf_edge_to_tme_faces.end())?&(it_ets->second):nullptr; + cip_handle_case_vertex(node_id,fset_bis,*it_edge,h_e,tm_f,tm_e); } else - cip_handle_case_vertex(node_id,&fset,*it_edge,h_2,tm1,tm2); + cip_handle_case_vertex(node_id,&fset,*it_edge,h_e,tm_f,tm_e); } } // end case ON_VERTEX break; @@ -821,8 +826,9 @@ class Intersection_of_triangle_meshes } }; + template void detect_intersections_in_the_graph(const TriangleMesh& tm, - const VertexPointMap& vpm, + const VPM& vpm, Node_id& current_node) { boost::unordered_map void construct_polylines(Output_iterator out){ - typedef typename boost::property_traits::value_type Point_3; + typedef typename boost::property_traits::value_type Point_3; std::size_t nb_nodes=nodes.size(); std::vector graph(nb_nodes); //counts the number of time each node has been seen @@ -1262,8 +1268,8 @@ class Intersection_of_triangle_meshes public: Intersection_of_triangle_meshes(const TriangleMesh& tm1, const TriangleMesh& tm2, - const VertexPointMap& vpm1, - const VertexPointMap& vpm2, + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2, const Node_visitor& v=Node_visitor()) : nodes(tm1, tm2, vpm1, vpm2) , visitor(v) @@ -1275,7 +1281,7 @@ public: // for autorefinement Intersection_of_triangle_meshes(const TriangleMesh& tm, - const VertexPointMap& vpm, + const VertexPointMap1& vpm, const Node_visitor& v=Node_visitor()) : nodes(tm, tm, vpm, vpm) , visitor(v) @@ -1293,8 +1299,8 @@ public: const TriangleMesh& tm1=nodes.tm1; const TriangleMesh& tm2=nodes.tm2; - const VertexPointMap& vpm1=nodes.vpm1; - const VertexPointMap& vpm2=nodes.vpm2; + const VertexPointMap1& vpm1=nodes.vpm1; + const VertexPointMap2& vpm2=nodes.vpm2; filter_intersections(tm1, tm2, vpm1, vpm2, throw_on_self_intersection); filter_intersections(tm2, tm1, vpm2, vpm1, throw_on_self_intersection); @@ -1308,6 +1314,7 @@ public: compute_intersection_of_coplanar_faces(current_node, tm1, tm2, vpm1, vpm2); else compute_intersection_of_coplanar_faces(current_node, tm2, tm1, vpm2, vpm1); + visitor.set_number_of_intersection_points_from_coplanar_faces(current_node+1); if (!coplanar_faces.empty()) visitor.input_have_coplanar_faces(); @@ -1324,6 +1331,7 @@ public: compute_intersection_points(tm1_edge_to_tm2_faces, tm1, tm2, vpm1, vpm2, current_node); compute_intersection_points(tm2_edge_to_tm1_faces, tm2, tm1, vpm2, vpm1, current_node); + if (!build_polylines){ visitor.finalize(nodes,tm1,tm2,vpm1,vpm2); return output; @@ -1361,7 +1369,7 @@ public: CGAL_assertion(doing_autorefinement); const TriangleMesh& tm=nodes.tm1; - const VertexPointMap& vpm=nodes.vpm1; + const VertexPointMap1& vpm=nodes.vpm1; filter_intersections(tm, vpm); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h index 9d43f9cc043..84ff44e7c8b 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h @@ -28,23 +28,26 @@ namespace Corefinement { // polylines. Different specializations are available depending whether // predicates on constructions are needed. template ::value_type + typename boost::property_traits::value_type >::Kernel::FT >::value > class Intersection_nodes; //Store only the double version of the intersection points. template -class Intersection_nodes + class VertexPointMap1, class VertexPointMap2> +class Intersection_nodes { //typedefs - typedef typename boost::property_traits::value_type Point_3; + typedef typename boost::property_traits::value_type Point_3; + CGAL_static_assertion((std::is_same::value_type, + typename boost::property_traits::value_type>::value)); + typedef typename Kernel_traits::Kernel Input_kernel; typedef std::vector Nodes_vector; typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; @@ -67,12 +70,13 @@ class Intersection_nodes public: const TriangleMesh &tm1, &tm2; - VertexPointMap vpm1, vpm2; + const VertexPointMap1& vpm1; + const VertexPointMap2& vpm2; Intersection_nodes(const TriangleMesh& tm1_, const TriangleMesh& tm2_, - const VertexPointMap& vpm1_, - const VertexPointMap& vpm2_) + const VertexPointMap1& vpm1_, + const VertexPointMap2& vpm2_) : tm1(tm1_) , tm2(tm2_) , vpm1(vpm1_) @@ -101,8 +105,8 @@ public: face_descriptor f_b, const TriangleMesh& tm_a, const TriangleMesh& tm_b, - const VertexPointMap vpm_a, - const VertexPointMap& vpm_b) + const VertexPointMap1& vpm_a, + const VertexPointMap2& vpm_b) { halfedge_descriptor h_b = halfedge(f_b, tm_b); add_new_node( @@ -119,7 +123,8 @@ public: add_new_node(edge_1, face_2, tm1, tm2, vpm1, vpm2); } - void call_put(const VertexPointMap& vpm, vertex_descriptor vd, std::size_t i, TriangleMesh&) + template // VertexPointMap1 or VertexPointMap2 + void call_put(const VPM& vpm, vertex_descriptor vd, std::size_t i, TriangleMesh&) { put(vpm, vd, nodes[i]); } @@ -132,14 +137,18 @@ public: // second specializations: store an exact copy of the points so // that we can answer exactly predicates -template -class Intersection_nodes +template +class Intersection_nodes { //typedefs public: typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; + private: - typedef typename boost::property_traits::value_type Point_3; + typedef typename boost::property_traits::value_type Point_3; + CGAL_static_assertion((std::is_same::value_type, + typename boost::property_traits::value_type>::value)); + typedef typename Kernel_traits::Kernel Input_kernel; typedef Cartesian_converter Double_to_exact; @@ -160,14 +169,16 @@ private: Exact_kernel::Intersect_3 exact_intersection; std::vector tm1_vertices, tm2_vertices; const bool doing_autorefinement; + public: const TriangleMesh &tm1, &tm2; - VertexPointMap vpm1, vpm2; + const VertexPointMap1& vpm1; + const VertexPointMap2& vpm2; Intersection_nodes(const TriangleMesh& tm1_, const TriangleMesh& tm2_, - const VertexPointMap& vpm1_, - const VertexPointMap& vpm2_) + const VertexPointMap1& vpm1_, + const VertexPointMap2& vpm2_) : doing_autorefinement(&tm1_ == &tm2_) , tm1(tm1_) , tm2(tm2_) @@ -211,12 +222,13 @@ public: //add a new node in the final graph. //it is the intersection of the triangle with the segment + template // VertexPointMap1 or VertexPointMap2 void add_new_node(halfedge_descriptor h_a, face_descriptor f_b, const TriangleMesh& tm_a, const TriangleMesh& tm_b, - const VertexPointMap vpm_a, - const VertexPointMap& vpm_b) + const VPM_A vpm_a, + const VPM_B vpm_b) { halfedge_descriptor h_b = halfedge(f_b, tm_b); add_new_node( @@ -229,11 +241,12 @@ public: } // use to resolve intersection of 3 faces in autorefinement only + template void add_new_node(halfedge_descriptor h1, halfedge_descriptor h2, halfedge_descriptor h3, const TriangleMesh& tm, - const VertexPointMap& vpm) + const VPM& vpm) { // TODO Far from optimal! typedef Exact_kernel::Plane_3 Plane_3; @@ -273,7 +286,8 @@ public: tm2_vertices.resize(enodes.size(), GT::null_vertex()); } - void call_put(const VertexPointMap& vpm, vertex_descriptor vd, std::size_t i, TriangleMesh& tm) + template // VertexPointMap1 or VertexPointMap2 + void call_put(const VPM& vpm, vertex_descriptor vd, std::size_t i, TriangleMesh& tm) { put(vpm, vd, exact_to_double(enodes[i])); if (&tm1==&tm) @@ -306,11 +320,16 @@ public: //Third specialization: The kernel already has exact constructions. -template -class Intersection_nodes +template +class Intersection_nodes { //typedefs - typedef typename boost::property_traits::value_type Point_3; + typedef typename boost::property_traits::value_type Point_3; + CGAL_static_assertion((std::is_same::value_type, + typename boost::property_traits::value_type>::value)); + typedef typename Kernel_traits::Kernel Input_kernel; typedef std::vector Nodes_vector; @@ -326,12 +345,13 @@ public: typedef Input_kernel Exact_kernel; const TriangleMesh &tm1, &tm2; - VertexPointMap vpm1, vpm2; + const VertexPointMap1& vpm1; + const VertexPointMap2& vpm2; Intersection_nodes(const TriangleMesh& tm1_, const TriangleMesh& tm2_, - const VertexPointMap& vpm1_, - const VertexPointMap& vpm2_) + const VertexPointMap1& vpm1_, + const VertexPointMap2& vpm2_) : tm1(tm1_) , tm2(tm2_) , vpm1(vpm1_) @@ -346,11 +366,12 @@ public: size_t size() const {return nodes.size();} const Point_3& exact_node(std::size_t i) const {return nodes[i];} + template void add_new_node(halfedge_descriptor h1, halfedge_descriptor h2, halfedge_descriptor h3, const TriangleMesh& tm, - const VertexPointMap& vpm) + const VPM& vpm) { // TODO Far from optimal! typedef typename Exact_kernel::Plane_3 Plane_3; @@ -380,8 +401,8 @@ public: face_descriptor f_b, const TriangleMesh& tm_a, const TriangleMesh& tm_b, - const VertexPointMap vpm_a, - const VertexPointMap& vpm_b) + const VertexPointMap1& vpm_a, + const VertexPointMap2& vpm_b) { halfedge_descriptor h_b=halfedge(f_b,tm_b); @@ -407,7 +428,8 @@ public: const Point_3& to_exact(const Point_3& p) const { return p; } - void call_put(const VertexPointMap& vpm, vertex_descriptor vd, std::size_t i, TriangleMesh&) + template // VertexPointMap1 or VertexPointMap2 + void call_put(const VPM& vpm, vertex_descriptor vd, std::size_t i, TriangleMesh&) { put(vpm, vd, nodes[i]); } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_of_coplanar_triangles_3.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_of_coplanar_triangles_3.h index d4645badd45..3c4efee0abb 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_of_coplanar_triangles_3.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_of_coplanar_triangles_3.h @@ -27,10 +27,15 @@ namespace CGAL{ namespace Polygon_mesh_processing { namespace Corefinement{ -template -struct Intersect_coplanar_faces_3{ +template +struct Intersect_coplanar_faces_3 +{ // typedefs - typedef typename boost::property_traits::value_type Point; + typedef typename boost::property_traits::value_type Point; + + CGAL_static_assertion((std::is_same::value_type, + typename boost::property_traits::value_type>::value)); + typedef typename CGAL::Kernel_traits::Kernel Input_kernel; typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; @@ -40,12 +45,14 @@ struct Intersect_coplanar_faces_3{ typedef Coplanar_intersection Inter_pt_info; // data members const TriangleMesh &tm1, &tm2; - const VertexPointMap &vpm1, &vpm2; + const VertexPointMap1& vpm1; + const VertexPointMap2& vpm2; + // constructor Intersect_coplanar_faces_3(const TriangleMesh& tm1_, const TriangleMesh& tm2_, - const VertexPointMap& vpm1_, - const VertexPointMap& vpm2_) + const VertexPointMap1& vpm1_, + const VertexPointMap2& vpm2_) : tm1(tm1_), tm2(tm2_), vpm1(vpm1_), vpm2(vpm2_) {} @@ -282,14 +289,14 @@ struct Intersect_coplanar_faces_3{ } }; -template +template void intersection_coplanar_faces( typename boost::graph_traits::face_descriptor f1, typename boost::graph_traits::face_descriptor f2, const TriangleMesh& tm1, const TriangleMesh& tm2, - const VertexPointMap& vpm1, - const VertexPointMap& vpm2, + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2, std::list< Coplanar_intersection >& inter_pts) { typedef boost::graph_traits GT; @@ -297,7 +304,7 @@ void intersection_coplanar_faces( halfedge_descriptor h1=halfedge(f1,tm1), h2=halfedge(f2,tm2); - Intersect_coplanar_faces_3 + Intersect_coplanar_faces_3 intersect_cpln(tm1, tm2, vpm1, vpm2); // We will add in `inter_pts` the initial triangle of h1 diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/predicates.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/predicates.h index a296e8de90a..4b25b21e472 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/predicates.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/predicates.h @@ -122,15 +122,15 @@ bool are_triangles_coplanar_same_side( } -template +template bool are_triangles_coplanar_same_side(Node_id o_prime_index, Node_id o_index, Node_id p_index, Node_id q_index, vertex_descriptor p, vertex_descriptor q, - const Vpm& vpm_p, - const Vpm& vpm_q, + const VPMP& vpm_p, + const VPMQ& vpm_q, const Node_vector& nodes) { const Node_id NID((std::numeric_limits::max)()); @@ -142,7 +142,7 @@ bool are_triangles_coplanar_same_side(Node_id o_prime_index, ); } -template +template bool sorted_around_edge( Node_id o_prime_index, Node_id o_index, Node_id p1_index, @@ -151,8 +151,8 @@ bool sorted_around_edge( Node_id o_prime_index, vertex_descriptor p1, vertex_descriptor p2, vertex_descriptor q, - const Vpm& vpm_p, - const Vpm& vpm_q, + const VPMP& vpm_p, + const VPMQ& vpm_q, const Node_vector& nodes) { const Node_id NID((std::numeric_limits::max)()); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h index ebaf53ffa92..9dfe7c560b5 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h @@ -1666,21 +1666,18 @@ surface_intersection(const TriangleMesh& tm1, const bool throw_on_self_intersection = parameters::choose_parameter(parameters::get_parameter(np1, internal_np::throw_on_self_intersection), false); - typedef typename GetVertexPointMap::const_type Vpm; - typedef typename GetVertexPointMap::const_type Vpm2; - CGAL_USE_TYPE(Vpm2); - CGAL_assertion_code( - static const bool same_vpm = (boost::is_same::value);) - CGAL_static_assertion(same_vpm); + typedef typename GetVertexPointMap::const_type VPM1; + typedef typename GetVertexPointMap::const_type VPM2; - Vpm vpm1 = parameters::choose_parameter(parameters::get_parameter(np1, internal_np::vertex_point), - get_const_property_map(CGAL::vertex_point, tm1)); - Vpm vpm2 = parameters::choose_parameter(parameters::get_parameter(np2, internal_np::vertex_point), - get_const_property_map(CGAL::vertex_point, tm2)); + CGAL_static_assertion((std::is_same::value_type, + typename boost::property_traits::value_type>::value)); - Corefinement::Intersection_of_triangle_meshes + VPM1 vpm1 = parameters::choose_parameter(parameters::get_parameter(np1, internal_np::vertex_point), + get_const_property_map(CGAL::vertex_point, tm1)); + VPM2 vpm2 = parameters::choose_parameter(parameters::get_parameter(np2, internal_np::vertex_point), + get_const_property_map(CGAL::vertex_point, tm2)); + + Corefinement::Intersection_of_triangle_meshes functor(tm1, tm2, vpm1, vpm2); return functor(polyline_output, throw_on_self_intersection, true); } @@ -1719,17 +1716,14 @@ surface_self_intersection(const TriangleMesh& tm, const NamedParameters& np) { // Vertex point maps - typedef typename GetVertexPointMap::const_type Vpm; + typedef typename GetVertexPointMap::const_type VPM; - Vpm vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point), - get_const_property_map(CGAL::vertex_point, tm)); + VPM vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point), + get_const_property_map(CGAL::vertex_point, tm)); // surface intersection algorithm call - typedef Corefinement::Default_surface_intersection_visitor Visitor; - Corefinement::Intersection_of_triangle_meshes - functor(tm, vpm); + typedef Corefinement::Default_surface_intersection_visitor Visitor; + Corefinement::Intersection_of_triangle_meshes functor(tm, vpm); polyline_output=functor(polyline_output, true); return polyline_output; From fa2f730832995ce6ef8d8559ab942d146e2339e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 8 May 2020 16:57:54 +0200 Subject: [PATCH 03/27] Add experimental clipping with allowed self-intersections in input --- .../CGAL/Polygon_mesh_processing/clip.h | 5 +- .../internal/Corefinement/clip_experimental.h | 501 ++++++++++++++++++ 2 files changed, 505 insertions(+), 1 deletion(-) create mode 100644 Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h index 01777e9740c..ccf45d4f86d 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h @@ -859,6 +859,9 @@ void split(TriangleMesh& tm, /// \endcond -} } //end of namespace CGAL::Polygon_mesh_processing +} // namespace Polygon_mesh_processing +} // namespace CGAL + +#include #endif // CGAL_POLYGON_MESH_PROCESSING_CLIP_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h new file mode 100644 index 00000000000..53395a044a6 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -0,0 +1,501 @@ +// Copyright (c) 2019-2020 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Mael Rouxel-Labbé +// Sebastien Loriot + +#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CLIP_EXPERIMENTAL_H +#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CLIP_EXPERIMENTAL_H + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace CGAL { +namespace Polygon_mesh_processing { +namespace internal { + +template +void output_mesh(const char* filename, const PolygonMesh& pmesh) +{ + std::ofstream output(filename); + output.precision(17); + output << pmesh; + output.close(); +} + +template +bool clip_mesh_exactly(TriangleMesh& tm, + EVPM tm_evpm, + TriangleMesh clipper, // intentional copy + const NamedParameters1& np_tm, + const NamedParameters2& np_c) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + // These typedefs concern the VPM of the CLIPPER mesh + typedef typename GetK::Kernel Clipper_kernel; + typedef typename GetVertexPointMap::type Clipper_VPM; + +#if 1 + typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; + typedef typename Exact_kernel::Point_3 Exact_point; + + typedef typename CGAL::Cartesian_converter C2E; +#else + typedef typename CGAL::Exact_kernel_selector::Exact_kernel Exact_kernel; + typedef typename Exact_kernel::Point_3 Exact_point; + + typedef typename CGAL::Exact_kernel_selector::C2E C2E; +#endif + + typedef CGAL::dynamic_vertex_property_t EP_property_tag; + + // must have equal exact point types because of corefinement.h + CGAL_static_assertion((std::is_same::type>::value)); + + using parameters::get_parameter; + using parameters::choose_parameter; + + // Initialize an exact vpm for the clipper + Clipper_VPM clipper_vpm = choose_parameter(get_parameter(np_c, internal_np::vertex_point), + get_const_property_map(CGAL::vertex_point, clipper)); + + C2E to_exact; + EVPM clipper_evpm = get(EP_property_tag(), clipper); + for(vertex_descriptor vd : vertices(clipper)) + put(clipper_evpm, vd, to_exact(get(clipper_vpm, vd))); + +#ifdef CGAL_DEBUG_CLIPPING + const bool valid_input = is_valid(tm) && is_valid(clipper) && + !does_self_intersect(tm, parameters::vertex_point_map(tm_evpm)) && + !does_self_intersect(clipper, parameters::vertex_point_map(clipper_evpm)) && + does_bound_a_volume(clipper, parameters::vertex_point_map(clipper_evpm)); + if(!valid_input) + { + std::cerr << "Invalid input for clip()" << std::endl; + std::cerr << "is tm valid: " << tm.is_valid() << std::endl; + std::cerr << "is clipper valid: " << clipper.is_valid() << std::endl; + std::cerr << "does part self intersect? " << does_self_intersect(tm, parameters::vertex_point_map(tm_evpm)) << std::endl; + std::cerr << "does clipper self intersect? " << does_self_intersect(clipper, parameters::vertex_point_map(clipper_evpm)) << std::endl; + std::cerr << "clipper bounds a volume? " << does_bound_a_volume(clipper) << std::endl; + std::exit(1); + } +#endif + + bool res = clip(tm, clipper, + parameters::vertex_point_map(tm_evpm).throw_on_self_intersection(true), + parameters::vertex_point_map(clipper_evpm)); + + CGAL_postcondition(CGAL::is_valid_polygon_mesh(tm)); + CGAL_postcondition(tm.is_valid()); + CGAL_postcondition(clipper.is_valid()); + + return res; +} + +template +void fill_triangle_mesh(typename boost::graph_traits::halfedge_descriptor hd, + const TriangleMesh& tm, + VPM tm_vpm, + EVPM tm_evpm, + TriangleMesh& si_face, + EVPM si_face_evpm) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + typedef typename GetVertexPointMap::type F_VPM; + F_VPM si_face_vpm = get_property_map(CGAL::vertex_point, si_face); + + vertex_descriptor vd0 = add_vertex(si_face); + vertex_descriptor vd1 = add_vertex(si_face); + vertex_descriptor vd2 = add_vertex(si_face); + put(si_face_vpm, vd0, get(tm_vpm, source(hd, tm))); + put(si_face_vpm, vd1, get(tm_vpm, target(hd, tm))); + put(si_face_vpm, vd2, get(tm_vpm, target(next(hd, tm), tm))); + put(si_face_evpm, vd0, get(tm_evpm, source(hd, tm))); + put(si_face_evpm, vd1, get(tm_evpm, target(hd, tm))); + put(si_face_evpm, vd2, get(tm_evpm, target(next(hd, tm), tm))); + + CGAL::Euler::add_face(std::initializer_list{vd0, vd1, vd2}, si_face); +} + +template +bool clip_single_self_intersecting_cc(TriangleMesh& cc, + const FacePairRange& self_intersecting_faces, + TriangleMesh& clipper, + const NamedParameters1& np_cc, + const NamedParameters2& np_c) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename GetK::Kernel CC_kernel; + typedef typename GetVertexPointMap::type CC_VPM; + +#if 1 + typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; + typedef typename Exact_kernel::Point_3 Exact_point; + + typedef typename CGAL::Cartesian_converter C2E; + typedef typename CGAL::Cartesian_converter E2C; +#else + typedef typename CGAL::Exact_kernel_selector::Exact_kernel Exact_kernel; + typedef typename Exact_kernel::Point_3 Exact_point; + + typedef typename CGAL::Exact_kernel_selector::C2E C2E; + typedef typename CGAL::Exact_kernel_selector::E2C E2C; +#endif + + typedef CGAL::dynamic_vertex_property_t EP_property_tag; + typedef typename boost::property_map::type EVPM; + + typedef std::pair Pair_of_faces; + + using parameters::get_parameter; + using parameters::choose_parameter; + + CGAL_precondition(does_self_intersect(cc)); + +#ifdef CGAL_DEBUG_CLIPPING + std::cout << "CC self-intersects" << std::endl; +#endif + + bool res = true; + + // 1. Gather the problematic faces + // 2. Copy (one-by-one) the problematic faces into independent meshes + // 3. Remove the problematic faces from 'cc' + // 4. Clip (cc - problematic_faces) + // 2bis. Clip the problematic faces + // 5. Re-add the clipped problematic faces from step 2. + + // 1. ------------------------------ + std::set si_faces; + for(const Pair_of_faces& fp : self_intersecting_faces) + { + si_faces.insert(fp.first); + si_faces.insert(fp.second); + } + + std::size_t independent_faces_n = si_faces.size(); +#ifdef CGAL_DEBUG_CLIPPING + std::cout << independent_faces_n << " independent faces" << std::endl; +#endif + + // 2. ------------------------------ + // Switch to an exact VPM + CC_VPM cc_vpm = choose_parameter(get_parameter(np_cc, internal_np::vertex_point), + get_const_property_map(CGAL::vertex_point, cc)); + + C2E to_exact; + EVPM cc_evpm = get(EP_property_tag(), cc); + for(vertex_descriptor vd : vertices(cc)) + put(cc_evpm, vd, to_exact(get(cc_vpm, vd))); + + // forced to do that because clip() wants both meshes to be of the same type (rather than just requiring + // that the two VPMs have matching types) + std::vector clipped_si_faces; + clipped_si_faces.reserve(independent_faces_n); + std::vector si_face_evpms; + si_face_evpms.reserve(independent_faces_n); + + std::size_t counter = 0; + for(face_descriptor fd : si_faces) + { +#ifdef CGAL_DEBUG_CLIPPING + std::cout << " Build single mesh for face: " << fd << " (" << counter+1 << "/" << independent_faces_n << ")" << std::endl; +#endif + + const halfedge_descriptor hd = halfedge(fd, cc); + if(!is_degenerate_triangle_face(face(hd, cc), cc)) + { + ++counter; + + clipped_si_faces.resize(counter); + TriangleMesh& si_face = clipped_si_faces.back(); + si_face_evpms.resize(counter); + EVPM& si_face_evpm = si_face_evpms.back(); + si_face_evpm = get(EP_property_tag(), si_face); + + fill_triangle_mesh(hd, cc, cc_vpm, cc_evpm, si_face, si_face_evpm); + } + } +#ifdef CGAL_DEBUG_CLIPPING + std::cout << clipped_si_faces.size() << " problematic faces to clip" << std::endl; +#endif + + // 3. ------------------------------ + for(const face_descriptor fd : si_faces) + CGAL::Euler::remove_face(halfedge(fd, cc), cc); + + // 3bis. iteratively remove faces incident to non-manifold vertices + std::vector non_manifold_cones; + CGAL::Polygon_mesh_processing::non_manifold_vertices(cc, std::back_inserter(non_manifold_cones)); + std::set nm_vertices; + + std::unordered_set nm_faces; + for(halfedge_descriptor nm_hd : non_manifold_cones) + { + // keep the faces incident to the nm vertex for only one (the first one) of the cones + if(nm_vertices.insert(target(nm_hd, cc)).second) + continue; + + for(halfedge_descriptor hd : CGAL::halfedges_around_target(nm_hd, cc)) + if(cc.face(hd) != face_descriptor()) + nm_faces.insert(cc.face(hd)); + } + + clipped_si_faces.reserve(clipped_si_faces.size() + nm_faces.size()); + for(face_descriptor fd : nm_faces) + { + clipped_si_faces.resize(clipped_si_faces.size() + 1); + si_face_evpms.resize(si_face_evpms.size() + 1); + + TriangleMesh& si_face = clipped_si_faces.back(); + EVPM si_face_evpm = si_face_evpms.back(); + + const halfedge_descriptor hd = halfedge(fd, cc); + fill_triangle_mesh(hd, cc, cc_vpm, cc_evpm, si_face, si_face_evpm); + CGAL::Euler::remove_face(hd, cc); + } + + CGAL_postcondition(is_valid(cc)); + CGAL_postcondition(!does_self_intersect(cc)); +#ifdef CGAL_DEBUG_CLIPPING + std::cout << "Pruned mesh: " << vertices(cc).size() << " nv " << faces(cc).size() << " nf " << std::endl; + output_mesh("pruned_mesh.off", cc); +#endif + + // 4. ------------------------------ +#ifdef CGAL_DEBUG_CLIPPING + std::cout << "Clipping CC's main part (w/o self-intersecting faces)" << std::endl; +#endif + + if(!clip_mesh_exactly(cc, cc_evpm, clipper, np_cc, np_c)) + res = false; + +#ifdef CGAL_DEBUG_CLIPPING + output_mesh("pruned_mesh_clipped.off", cc); +#endif + + // 5. ------------------------------ + counter = 0; + for(std::size_t i=0, csfs=clipped_si_faces.size(); i +bool clip_single_cc(TriangleMesh& cc, + TriangleMesh& clipper, + const NamedParameters1& np_cc, + const NamedParameters2& np_c) +{ + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef std::pair Pair_of_faces; + +#ifdef CGAL_DEBUG_CLIPPING + std::cout << "Clip single CC of size " << vertices(cc).size() << " nv " << faces(cc).size() << " nf" << std::endl; +#endif + + std::vector intersecting_faces; + self_intersections(cc, std::back_inserter(intersecting_faces), np_cc); + + if(intersecting_faces.empty()) + return clip(cc, clipper, np_cc, np_c); + else + return clip_single_self_intersecting_cc(cc, intersecting_faces, clipper, np_cc, np_c); +} + +} // namespace internal + +namespace experimental { + +// Try to split the mesh in multiple connected components and clip them independently +template +void clip_self_intersecting_mesh(TriangleMesh& tm, + TriangleMesh& clipper, + const NamedParameters1& np_tm, + const NamedParameters2& np_c) +{ + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename boost::graph_traits::faces_size_type faces_size_type; + typedef typename boost::property_map >::type PatchIDMap; + + // These typedefs concern the VPM of the CLIPPED mesh + typedef typename GetVertexPointMap::type VPM; + typedef typename boost::property_traits::value_type Point; + + typedef CGAL::dynamic_vertex_property_t P_property_tag; + typedef typename boost::property_map::type DVPM; + + PatchIDMap pidmap = get(CGAL::face_patch_id_t(), tm); + faces_size_type num_cc = connected_components(tm, pidmap); + +#ifdef CGAL_DEBUG_CLIPPING + std::cout << num_cc << " connected component(s)" << std::endl; +#endif + + // order ccs by size + std::vector number_of_faces_in_cc(num_cc, 0); + for(face_descriptor f : faces(tm)) + number_of_faces_in_cc[get(pidmap, f)] += 1; + + std::vector ordered_cc_ids(num_cc); + std::iota(std::begin(ordered_cc_ids), std::end(ordered_cc_ids), 0); + std::sort(std::begin(ordered_cc_ids), std::end(ordered_cc_ids), + [&number_of_faces_in_cc](const int i, const int j) -> bool { + return (number_of_faces_in_cc[i] > number_of_faces_in_cc[j]); + }); + + for(std::size_t i=0; i ccs(num_cc - 1); // largest CC stays in 'tm' + std::vector cc_vpms(num_cc - 1); + + // Extract all but the largest CC + for(faces_size_type i=1; i tm_cc(tm, ordered_cc_ids[i], pidmap); + CGAL::copy_face_graph(tm_cc, ccs[ordered_cc_ids[i]], + np_tm, parameters::vertex_point_map(cc_vpms[i])); + } + + const CGAL::Bbox_3 clipper_bbox = CGAL::Polygon_mesh_processing::bbox(clipper, np_c); + + // Clip CC by CC + keep_connected_components(tm, std::vector(1, ordered_cc_ids[0]), pidmap); + + const CGAL::Bbox_3 tm_bbox = CGAL::Polygon_mesh_processing::bbox(tm, np_tm); + if(CGAL::do_overlap(tm_bbox, clipper_bbox)) + { + // @fixme likely don't want to forward clip_volumes to this single CC + internal::clip_single_cc(tm, clipper, np_tm, np_c); + } + + for(faces_size_type i=1; i +void generic_clip(TriangleMesh& tm, + TriangleMesh& clipper, + const NamedParameters1& np_tm, + const NamedParameters2& np_c) +{ + // @todo: do the self-intersection test only once + if(!does_self_intersect(tm, np_tm)) + clip(tm, clipper, np_tm, np_c); + else + clip_self_intersecting_mesh(tm, clipper, np_tm, np_c); +} + +template +void generic_clip(TriangleMesh& tm, TriangleMesh& clipper, const NamedParameters& np_tm) +{ + return generic_clip(tm, clipper, np_tm, parameters::all_default()); +} + +template +void generic_clip(TriangleMesh& tm, TriangleMesh& clipper) +{ + return generic_clip(tm, clipper, parameters::all_default(), parameters::all_default()); +} + +} // namespace experimental +} // namespace Polygon_mesh_processing +} // namespace CGAL + + +#endif //CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CLIP_EXPERIMENTAL_H From 091d1ec1952cc94c1362f7ae0dea0752205e1d4b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 11 May 2020 17:29:24 +0200 Subject: [PATCH 04/27] Clarify variable names to avoid confusion between f from tmf and faces --- .../internal/Corefinement/intersection_impl.h | 132 +++++++++--------- 1 file changed, 66 insertions(+), 66 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h index 12035e0b01b..ab840fe75ec 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h @@ -530,14 +530,14 @@ class Intersection_of_triangle_meshes } } - template + template void compute_intersection_of_coplanar_faces(Node_id& current_node, - const TriangleMesh& tm_f, - const TriangleMesh& tm_e, - const VPMF& vpm_f, - const VPME& vpm_e) + const TriangleMesh& tm1, + const TriangleMesh& tm2, + const VPM1& vpm1, + const VPM2& vpm2) { - CGAL_assertion(&tm_f < &tm_e || &tm_f == &tm_e); + CGAL_assertion( &tm1 < &tm2 || &tm1==&tm2 ); typedef std::tuple Cpl_inter_pt; std::list inter_pts; // compute the intersection points between the two coplanar faces - intersection_coplanar_faces(f1, f2, tm_f, tm_e, vpm_f, vpm_e, inter_pts); + intersection_coplanar_faces(f1, f2, tm1, tm2, vpm1, vpm2, inter_pts); std::size_t nb_pts=inter_pts.size(); std::vector cpln_nodes; cpln_nodes.reserve(nb_pts); @@ -569,7 +569,7 @@ class Intersection_of_triangle_meshes Node_id node_id; bool is_new_node; std::tie(node_id, is_new_node) = - get_or_create_node(ipt,current_node,coplanar_node_map,tm_f,tm_e); + get_or_create_node(ipt,current_node,coplanar_node_map,tm1,tm2); cpln_nodes.push_back(node_id); switch(ipt.type_1){ @@ -577,13 +577,13 @@ class Intersection_of_triangle_meshes { switch(ipt.type_2){ case ON_VERTEX: - handle_coplanar_case_VERTEX_VERTEX(ipt.info_1,ipt.info_2,tm_f,tm_e,node_id,is_new_node); + handle_coplanar_case_VERTEX_VERTEX(ipt.info_1,ipt.info_2,tm1,tm2,node_id,is_new_node); break; case ON_EDGE: - handle_coplanar_case_VERTEX_EDGE(ipt.info_1,ipt.info_2,tm_f,tm_e,node_id,is_new_node); + handle_coplanar_case_VERTEX_EDGE(ipt.info_1,ipt.info_2,tm1,tm2,node_id,is_new_node); break; case ON_FACE: - handle_coplanar_case_VERTEX_FACE(ipt.info_1,ipt.info_2,tm_f,tm_e,node_id,is_new_node); + handle_coplanar_case_VERTEX_FACE(ipt.info_1,ipt.info_2,tm1,tm2,node_id,is_new_node); break; default: CGAL_error_msg("Should not get there!"); } @@ -593,15 +593,15 @@ class Intersection_of_triangle_meshes { switch(ipt.type_2){ case ON_VERTEX: - handle_coplanar_case_VERTEX_EDGE(ipt.info_2,ipt.info_1,tm_e,tm_f,node_id,is_new_node); + handle_coplanar_case_VERTEX_EDGE(ipt.info_2,ipt.info_1,tm2,tm1,node_id,is_new_node); break; case ON_EDGE: { if(is_new_node) - visitor.new_node_added(node_id,ON_EDGE,ipt.info_1,ipt.info_2,tm_f,tm_e,false,false); - typename Edge_to_faces::iterator it_ets=stm_edge_to_ltm_faces.find(edge(ipt.info_1,tm_f)); + visitor.new_node_added(node_id,ON_EDGE,ipt.info_1,ipt.info_2,tm1,tm2,false,false); + typename Edge_to_faces::iterator it_ets=stm_edge_to_ltm_faces.find(edge(ipt.info_1,tm1)); Face_set* fset = (it_ets!=stm_edge_to_ltm_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_edge(node_id,fset,ipt.info_1,ipt.info_2,tm_f,tm_e); + cip_handle_case_edge(node_id,fset,ipt.info_1,ipt.info_2,tm1,tm2); } break; default: CGAL_error_msg("Should not get there!"); @@ -611,7 +611,7 @@ class Intersection_of_triangle_meshes case ON_FACE: { CGAL_assertion(ipt.type_2==ON_VERTEX); - handle_coplanar_case_VERTEX_FACE(ipt.info_2,ipt.info_1,tm_e,tm_f,node_id,is_new_node); + handle_coplanar_case_VERTEX_FACE(ipt.info_2,ipt.info_1,tm2,tm1,node_id,is_new_node); } break; default: CGAL_error_msg("Should not get there!"); @@ -640,66 +640,66 @@ class Intersection_of_triangle_meshes //add a new node in the final graph. //it is the intersection of the triangle with the segment - template - void add_new_node(halfedge_descriptor h_f, - face_descriptor f_e, - const TriangleMesh& tm_f, - const TriangleMesh& tm_e, - const VPMF& vpm_f, - const VPME& vpm_e, + template + void add_new_node(halfedge_descriptor h_1, + face_descriptor f_2, + const TriangleMesh& tm1, + const TriangleMesh& tm2, + const VPM1& vpm1, + const VPM2& vpm2, std::tuple inter_res) { if ( std::get<3>(inter_res) ) // is edge target in triangle plane - nodes.add_new_node(get(vpm_f, target(h_f, tm_f))); + nodes.add_new_node(get(vpm1, target(h_1,tm1))); else{ if (std::get<2>(inter_res)) // is edge source in triangle plane - nodes.add_new_node(get(vpm_f, source(h_f, tm_f))); + nodes.add_new_node(get(vpm1, source(h_1,tm1))); else - nodes.add_new_node(h_f, f_e, tm_f, tm_e, vpm_f, vpm_e); + nodes.add_new_node(h_1,f_2,tm1,tm2,vpm1,vpm2); } } - template - void compute_intersection_points(Edge_to_faces& tmf_edge_to_tme_faces, - const TriangleMesh& tm_f, - const TriangleMesh& tm_e, - const VPMF& vpm_f, - const VPME& vpm_e, + template + void compute_intersection_points(Edge_to_faces& tm1_edge_to_tm2_faces, + const TriangleMesh& tm1, + const TriangleMesh& tm2, + const VPM1& vpm1, + const VPM2& vpm2, Node_id& current_node) { typedef std::tuple Inter_type; - for(typename Edge_to_faces::iterator it=tmf_edge_to_tme_faces.begin(); - it!=tmf_edge_to_tme_faces.end();++it) + for(typename Edge_to_faces::iterator it=tm1_edge_to_tm2_faces.begin(); + it!=tm1_edge_to_tm2_faces.end();++it) { - edge_descriptor e_f=it->first; - halfedge_descriptor h_f=halfedge(e_f,tm_f); + edge_descriptor e_1=it->first; + halfedge_descriptor h_1=halfedge(e_1,tm1); Face_set& fset=it->second; while (!fset.empty()){ - face_descriptor f_e=*fset.begin(); + face_descriptor f_2=*fset.begin(); - Inter_type res=intersection_type(h_f,f_e,tm_f,tm_e,vpm_f,vpm_e); + Inter_type res=intersection_type(h_1,f_2,tm1,tm2,vpm1,vpm2); Intersection_type type=std::get<0>(res); - //handle degenerate case: one extremity of edge belong to f_e + //handle degenerate case: one extremity of edge belong to f_2 std::vector all_edges; if ( std::get<3>(res) ) // is edge target in triangle plane - std::copy(halfedges_around_target(h_f,tm_f).first, - halfedges_around_target(h_f,tm_f).second, + std::copy(halfedges_around_target(h_1,tm1).first, + halfedges_around_target(h_1,tm1).second, std::back_inserter(all_edges)); else{ if ( std::get<2>(res) ) // is edge source in triangle plane - std::copy(halfedges_around_source(h_f,tm_f).first, - halfedges_around_source(h_f,tm_f).second, + std::copy(halfedges_around_source(h_1,tm1).first, + halfedges_around_source(h_1,tm1).second, std::back_inserter(all_edges)); else - all_edges.push_back(h_f); + all_edges.push_back(h_1); } - CGAL_precondition(all_edges[0]==h_f || all_edges[0]==opposite(h_f,tm_f)); + CGAL_precondition(all_edges[0]==h_1 || all_edges[0]==opposite(h_1,tm1)); // #ifdef USE_DETECTION_MULTIPLE_DEFINED_EDGES // check_coplanar_edges(std::next(all_edges.begin()), @@ -722,20 +722,20 @@ class Intersection_of_triangle_meshes // Case when the edge pierces the face in its interior. case ON_FACE: { - CGAL_assertion(f_e==face(std::get<1>(res),tm_e)); + CGAL_assertion(f_2==face(std::get<1>(res),tm2)); Node_id node_id=++current_node; - add_new_node(h_f,f_e,tm_f,tm_e,vpm_f,vpm_e,res); - visitor.new_node_added(node_id,ON_FACE,h_f,halfedge(f_e,tm_e),tm_f,tm_e,std::get<3>(res),std::get<2>(res)); + add_new_node(h_1,f_2,tm1,tm2,vpm1,vpm2,res); + visitor.new_node_added(node_id,ON_FACE,h_1,halfedge(f_2,tm2),tm1,tm2,std::get<3>(res),std::get<2>(res)); for (;it_edge!=all_edges.end();++it_edge){ - add_intersection_point_to_face_and_all_edge_incident_faces(f_e,*it_edge,tm_e,tm_f,node_id); + add_intersection_point_to_face_and_all_edge_incident_faces(f_2,*it_edge,tm2,tm1,node_id); //erase face from the list to test intersection with it_edge if ( it_edge==all_edges.begin() ) fset.erase(fset.begin()); else { - typename Edge_to_faces::iterator it_ets=tmf_edge_to_tme_faces.find(edge(*it_edge,tm_f)); - if(it_ets!=tmf_edge_to_tme_faces.end()) it_ets->second.erase(f_e); + typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(*it_edge,tm1)); + if(it_ets!=tm1_edge_to_tm2_faces.end()) it_ets->second.erase(f_2); } } } // end case ON_FACE @@ -745,17 +745,17 @@ class Intersection_of_triangle_meshes case ON_EDGE: { Node_id node_id=++current_node; - add_new_node(h_f,f_e,tm_f,tm_e,vpm_f,vpm_e,res); - halfedge_descriptor h_e=std::get<1>(res); - visitor.new_node_added(node_id,ON_EDGE,h_f,h_e,tm_f,tm_e,std::get<3>(res),std::get<2>(res)); + add_new_node(h_1,f_2,tm1,tm2,vpm1,vpm2,res); + halfedge_descriptor h_2=std::get<1>(res); + visitor.new_node_added(node_id,ON_EDGE,h_1,h_2,tm1,tm2,std::get<3>(res),std::get<2>(res)); for (;it_edge!=all_edges.end();++it_edge){ if ( it_edge!=all_edges.begin() ){ - typename Edge_to_faces::iterator it_ets=tmf_edge_to_tme_faces.find(edge(*it_edge,tm_f)); - Face_set* fset_bis = (it_ets!=tmf_edge_to_tme_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_edge(node_id,fset_bis,*it_edge,h_e,tm_f,tm_e); + typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(*it_edge,tm1)); + Face_set* fset_bis = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; + cip_handle_case_edge(node_id,fset_bis,*it_edge,h_2,tm1,tm2); } else - cip_handle_case_edge(node_id,&fset,*it_edge,h_e,tm_f,tm_e); + cip_handle_case_edge(node_id,&fset,*it_edge,h_2,tm1,tm2); } } // end case ON_EDGE break; @@ -763,18 +763,18 @@ class Intersection_of_triangle_meshes case ON_VERTEX: { Node_id node_id=++current_node; - halfedge_descriptor h_e=std::get<1>(res); - nodes.add_new_node(get(vpm_e, target(h_e,tm_e))); //we use the original vertex to create the node + halfedge_descriptor h_2=std::get<1>(res); + nodes.add_new_node(get(vpm2, target(h_2,tm2))); //we use the original vertex to create the node //before it was ON_FACE but do not remember why, probably a bug... - visitor.new_node_added(node_id,ON_VERTEX,h_f,h_e,tm_f,tm_e,std::get<3>(res),std::get<2>(res)); + visitor.new_node_added(node_id,ON_VERTEX,h_1,h_2,tm1,tm2,std::get<3>(res),std::get<2>(res)); for (;it_edge!=all_edges.end();++it_edge){ if ( it_edge!=all_edges.begin() ){ - typename Edge_to_faces::iterator it_ets=tmf_edge_to_tme_faces.find(edge(*it_edge,tm_f)); - Face_set* fset_bis = (it_ets!=tmf_edge_to_tme_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_vertex(node_id,fset_bis,*it_edge,h_e,tm_f,tm_e); + typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(*it_edge,tm1)); + Face_set* fset_bis = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; + cip_handle_case_vertex(node_id,fset_bis,*it_edge,h_2,tm1,tm2); } else - cip_handle_case_vertex(node_id,&fset,*it_edge,h_e,tm_f,tm_e); + cip_handle_case_vertex(node_id,&fset,*it_edge,h_2,tm1,tm2); } } // end case ON_VERTEX break; From a620191023bb327d7df2bdb5717d872703bd6e06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 11 May 2020 17:30:37 +0200 Subject: [PATCH 05/27] Improve NP usage in clipping with self intersections --- .../internal/Corefinement/clip_experimental.h | 78 +++++++++++-------- 1 file changed, 47 insertions(+), 31 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 53395a044a6..98fd69c6bac 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -47,10 +47,10 @@ void output_mesh(const char* filename, const PolygonMesh& pmesh) template -bool clip_mesh_exactly(TriangleMesh& tm, - EVPM tm_evpm, +bool clip_mesh_exactly(TriangleMesh& cc, + EVPM cc_evpm, TriangleMesh clipper, // intentional copy - const NamedParameters1& np_tm, + const NamedParameters1& np_cc, const NamedParameters2& np_c) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; @@ -89,29 +89,38 @@ bool clip_mesh_exactly(TriangleMesh& tm, put(clipper_evpm, vd, to_exact(get(clipper_vpm, vd))); #ifdef CGAL_DEBUG_CLIPPING - const bool valid_input = is_valid(tm) && is_valid(clipper) && - !does_self_intersect(tm, parameters::vertex_point_map(tm_evpm)) && + const bool valid_input = is_valid(cc) && is_valid(clipper) && + !does_self_intersect(cc, parameters::vertex_point_map(cc_evpm)) && !does_self_intersect(clipper, parameters::vertex_point_map(clipper_evpm)) && does_bound_a_volume(clipper, parameters::vertex_point_map(clipper_evpm)); if(!valid_input) { std::cerr << "Invalid input for clip()" << std::endl; - std::cerr << "is tm valid: " << tm.is_valid() << std::endl; + std::cerr << "is cc valid: " << cc.is_valid() << std::endl; std::cerr << "is clipper valid: " << clipper.is_valid() << std::endl; - std::cerr << "does part self intersect? " << does_self_intersect(tm, parameters::vertex_point_map(tm_evpm)) << std::endl; + std::cerr << "does part self intersect? " << does_self_intersect(cc, parameters::vertex_point_map(cc_evpm)) << std::endl; std::cerr << "does clipper self intersect? " << does_self_intersect(clipper, parameters::vertex_point_map(clipper_evpm)) << std::endl; std::cerr << "clipper bounds a volume? " << does_bound_a_volume(clipper) << std::endl; std::exit(1); } #endif - bool res = clip(tm, clipper, - parameters::vertex_point_map(tm_evpm).throw_on_self_intersection(true), + bool clip_volume = choose_parameter(get_parameter(np_cc, internal_np::clip_volume), true); + bool use_compact_clipper = choose_parameter(get_parameter(np_cc, internal_np::use_compact_clipper), false); + + // @todo is it possible to forward clip_volume/use_compact_clipper/visitor? + clip_volume = false; + use_compact_clipper = false; + + bool res = clip(cc, clipper, + parameters::vertex_point_map(cc_evpm) + .clip_volume(clip_volume) + .clip_volume(use_compact_clipper) + .throw_on_self_intersection(true), parameters::vertex_point_map(clipper_evpm)); - CGAL_postcondition(CGAL::is_valid_polygon_mesh(tm)); - CGAL_postcondition(tm.is_valid()); - CGAL_postcondition(clipper.is_valid()); + CGAL_postcondition(is_valid_polygon_mesh(cc)); + CGAL_postcondition(is_valid_polygon_mesh(clipper)); return res; } @@ -386,24 +395,32 @@ namespace experimental { // Try to split the mesh in multiple connected components and clip them independently template -void clip_self_intersecting_mesh(TriangleMesh& tm, +bool clip_self_intersecting_mesh(TriangleMesh& tm, TriangleMesh& clipper, const NamedParameters1& np_tm, const NamedParameters2& np_c) { typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef typename boost::graph_traits::faces_size_type faces_size_type; - typedef typename boost::property_map >::type PatchIDMap; - - // These typedefs concern the VPM of the CLIPPED mesh typedef typename GetVertexPointMap::type VPM; typedef typename boost::property_traits::value_type Point; typedef CGAL::dynamic_vertex_property_t P_property_tag; typedef typename boost::property_map::type DVPM; + typedef typename boost::graph_traits::faces_size_type faces_size_type; + typedef typename boost::property_map >::type PatchIDMap; + + using parameters::get_parameter; + using parameters::choose_parameter; + + // @todo keep this? + if(!does_self_intersect(tm, np_tm)) + return clip(tm, clipper, np_tm, np_c); + + // Splitting connected components avoids having to treat a lot of "easy" self-intersections + // (for example from two spheres intersecting) PatchIDMap pidmap = get(CGAL::face_patch_id_t(), tm); faces_size_type num_cc = connected_components(tm, pidmap); @@ -438,6 +455,7 @@ void clip_self_intersecting_mesh(TriangleMesh& tm, np_tm, parameters::vertex_point_map(cc_vpms[i])); } + bool res = true; const CGAL::Bbox_3 clipper_bbox = CGAL::Polygon_mesh_processing::bbox(clipper, np_c); // Clip CC by CC @@ -446,8 +464,8 @@ void clip_self_intersecting_mesh(TriangleMesh& tm, const CGAL::Bbox_3 tm_bbox = CGAL::Polygon_mesh_processing::bbox(tm, np_tm); if(CGAL::do_overlap(tm_bbox, clipper_bbox)) { - // @fixme likely don't want to forward clip_volumes to this single CC - internal::clip_single_cc(tm, clipper, np_tm, np_c); + if(!internal::clip_single_cc(tm, clipper, np_tm, np_c)) + res = false; } for(faces_size_type i=1; i -void generic_clip(TriangleMesh& tm, +bool generic_clip(TriangleMesh& tm, TriangleMesh& clipper, const NamedParameters1& np_tm, const NamedParameters2& np_c) { - // @todo: do the self-intersection test only once - if(!does_self_intersect(tm, np_tm)) - clip(tm, clipper, np_tm, np_c); - else - clip_self_intersecting_mesh(tm, clipper, np_tm, np_c); + return clip_self_intersecting_mesh(tm, clipper, np_tm, np_c); } template -void generic_clip(TriangleMesh& tm, TriangleMesh& clipper, const NamedParameters& np_tm) +bool generic_clip(TriangleMesh& tm, TriangleMesh& clipper, const NamedParameters& np_tm) { return generic_clip(tm, clipper, np_tm, parameters::all_default()); } template -void generic_clip(TriangleMesh& tm, TriangleMesh& clipper) +bool generic_clip(TriangleMesh& tm, TriangleMesh& clipper) { return generic_clip(tm, clipper, parameters::all_default(), parameters::all_default()); } From d556bd0474126cdba8ec8085f7954814c7793226 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 11 May 2020 18:39:59 +0200 Subject: [PATCH 06/27] Minor bug fixes --- .../internal/Corefinement/clip_experimental.h | 22 ++++++++----------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 98fd69c6bac..b3b907fac87 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -81,7 +81,7 @@ bool clip_mesh_exactly(TriangleMesh& cc, // Initialize an exact vpm for the clipper Clipper_VPM clipper_vpm = choose_parameter(get_parameter(np_c, internal_np::vertex_point), - get_const_property_map(CGAL::vertex_point, clipper)); + get_property_map(CGAL::vertex_point, clipper)); C2E to_exact; EVPM clipper_evpm = get(EP_property_tag(), clipper); @@ -101,7 +101,7 @@ bool clip_mesh_exactly(TriangleMesh& cc, std::cerr << "does part self intersect? " << does_self_intersect(cc, parameters::vertex_point_map(cc_evpm)) << std::endl; std::cerr << "does clipper self intersect? " << does_self_intersect(clipper, parameters::vertex_point_map(clipper_evpm)) << std::endl; std::cerr << "clipper bounds a volume? " << does_bound_a_volume(clipper) << std::endl; - std::exit(1); + return false; } #endif @@ -115,7 +115,7 @@ bool clip_mesh_exactly(TriangleMesh& cc, bool res = clip(cc, clipper, parameters::vertex_point_map(cc_evpm) .clip_volume(clip_volume) - .clip_volume(use_compact_clipper) + .use_compact_clipper(use_compact_clipper) .throw_on_self_intersection(true), parameters::vertex_point_map(clipper_evpm)); @@ -135,15 +135,11 @@ void fill_triangle_mesh(typename boost::graph_traits::halfedge_des { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename GetVertexPointMap::type F_VPM; - F_VPM si_face_vpm = get_property_map(CGAL::vertex_point, si_face); - vertex_descriptor vd0 = add_vertex(si_face); vertex_descriptor vd1 = add_vertex(si_face); vertex_descriptor vd2 = add_vertex(si_face); - put(si_face_vpm, vd0, get(tm_vpm, source(hd, tm))); - put(si_face_vpm, vd1, get(tm_vpm, target(hd, tm))); - put(si_face_vpm, vd2, get(tm_vpm, target(next(hd, tm), tm))); + + // Note that we don't even fill the internal VPM put(si_face_evpm, vd0, get(tm_evpm, source(hd, tm))); put(si_face_evpm, vd1, get(tm_evpm, target(hd, tm))); put(si_face_evpm, vd2, get(tm_evpm, target(next(hd, tm), tm))); @@ -153,8 +149,8 @@ void fill_triangle_mesh(typename boost::graph_traits::halfedge_des template -bool clip_single_self_intersecting_cc(TriangleMesh& cc, - const FacePairRange& self_intersecting_faces, +bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_faces, + TriangleMesh& cc, TriangleMesh& clipper, const NamedParameters1& np_cc, const NamedParameters2& np_c) @@ -219,7 +215,7 @@ bool clip_single_self_intersecting_cc(TriangleMesh& cc, // 2. ------------------------------ // Switch to an exact VPM CC_VPM cc_vpm = choose_parameter(get_parameter(np_cc, internal_np::vertex_point), - get_const_property_map(CGAL::vertex_point, cc)); + get_property_map(CGAL::vertex_point, cc)); C2E to_exact; EVPM cc_evpm = get(EP_property_tag(), cc); @@ -385,7 +381,7 @@ bool clip_single_cc(TriangleMesh& cc, if(intersecting_faces.empty()) return clip(cc, clipper, np_cc, np_c); else - return clip_single_self_intersecting_cc(cc, intersecting_faces, clipper, np_cc, np_c); + return clip_single_self_intersecting_cc(intersecting_faces, cc, clipper, np_cc, np_c); } } // namespace internal From ef22a4d0352d46a7f142b85784101b7d5cc9a1f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 11 May 2020 19:03:35 +0200 Subject: [PATCH 07/27] Add a NP to call experimental clip_w/_self_intersections from PMP::clip --- .../CGAL/boost/graph/parameters_interface.h | 1 + .../CGAL/Polygon_mesh_processing/clip.h | 31 ++++++++++++----- .../internal/Corefinement/clip_experimental.h | 34 +++++++++++-------- 3 files changed, 44 insertions(+), 22 deletions(-) diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h index 17cf689d262..02f5755a388 100644 --- a/BGL/include/CGAL/boost/graph/parameters_interface.h +++ b/BGL/include/CGAL/boost/graph/parameters_interface.h @@ -76,6 +76,7 @@ CGAL_add_named_parameter(projection_functor_t, projection_functor, projection_fu CGAL_add_named_parameter(throw_on_self_intersection_t, throw_on_self_intersection, throw_on_self_intersection) CGAL_add_named_parameter(clip_volume_t, clip_volume, clip_volume) CGAL_add_named_parameter(use_compact_clipper_t, use_compact_clipper, use_compact_clipper) +CGAL_add_named_parameter(input_might_have_self_intersection_t, input_might_have_self_intersection, input_might_have_self_intersection) CGAL_add_named_parameter(output_iterator_t, output_iterator, output_iterator) CGAL_add_named_parameter(erase_all_duplicates_t, erase_all_duplicates, erase_all_duplicates) CGAL_add_named_parameter(require_same_orientation_t, require_same_orientation, require_same_orientation) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h index ccf45d4f86d..3bdcde94489 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h @@ -15,6 +15,8 @@ #include +#include + #include #include #include @@ -34,11 +36,9 @@ #include -namespace CGAL{ +namespace CGAL { namespace Polygon_mesh_processing { - -namespace internal -{ +namespace internal { template int @@ -402,6 +402,17 @@ void split_along_edges(TriangleMesh& tm, } // end of internal namespace +namespace experimental { + +template +bool clip_self_intersecting_mesh(TriangleMesh& tm, + TriangleMesh& clipper, + const NamedParameters1& np_tm, + const NamedParameters2& np_c); + +} // namespace experimental + /** * \ingroup PMP_corefinement_grp * clips `tm` by keeping the part that is inside the volume \link coref_def_subsec bounded \endlink @@ -461,8 +472,14 @@ clip(TriangleMesh& tm, const NamedParameters1& np_tm, const NamedParameters2& np_c) { - const bool clip_volume = - parameters::choose_parameter(parameters::get_parameter(np_tm, internal_np::clip_volume), false); + using parameters::choose_parameter; + using parameters::get_parameter; + + const bool clip_with_si = choose_parameter(get_parameter(np_tm, internal_np::input_might_have_self_intersection), false); + if(clip_with_si) + return experimental::clip_self_intersecting_mesh(tm, clipper, np_tm, np_c); + + const bool clip_volume = choose_parameter(get_parameter(np_tm, internal_np::clip_volume), false); if(clip_volume && is_closed(tm)) return corefine_and_compute_intersection(tm, clipper, tm, np_tm, np_c, np_tm); @@ -862,6 +879,4 @@ void split(TriangleMesh& tm, } // namespace Polygon_mesh_processing } // namespace CGAL -#include - #endif // CGAL_POLYGON_MESH_PROCESSING_CLIP_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index b3b907fac87..d956eb7bf7f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -16,6 +16,7 @@ #include +#include #include #include #include @@ -32,8 +33,23 @@ #include #include +#include +#include +#include +#include +#include +#include +#include + namespace CGAL { namespace Polygon_mesh_processing { + +template +bool clip(TriangleMesh& tm, + TriangleMesh& clipper, + const NamedParameters1& np_tm, + const NamedParameters2& np_c); + namespace internal { template @@ -483,26 +499,16 @@ bool clip_self_intersecting_mesh(TriangleMesh& tm, return res; } -template -bool generic_clip(TriangleMesh& tm, - TriangleMesh& clipper, - const NamedParameters1& np_tm, - const NamedParameters2& np_c) -{ - return clip_self_intersecting_mesh(tm, clipper, np_tm, np_c); -} - template -bool generic_clip(TriangleMesh& tm, TriangleMesh& clipper, const NamedParameters& np_tm) +bool clip_self_intersecting_mesh(TriangleMesh& tm, TriangleMesh& clipper, const NamedParameters& np_tm) { - return generic_clip(tm, clipper, np_tm, parameters::all_default()); + return clip_self_intersecting_mesh(tm, clipper, np_tm, parameters::all_default()); } template -bool generic_clip(TriangleMesh& tm, TriangleMesh& clipper) +bool clip_self_intersecting_mesh(TriangleMesh& tm, TriangleMesh& clipper) { - return generic_clip(tm, clipper, parameters::all_default(), parameters::all_default()); + return clip_self_intersecting_mesh(tm, clipper, parameters::all_default(), parameters::all_default()); } } // namespace experimental From 8b825eaaf9399ce49a15325d2d8cad6dabd5bf88 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 25 May 2020 13:14:01 +0200 Subject: [PATCH 08/27] Fix improper VPM forwarding --- .../internal/Corefinement/clip_experimental.h | 21 ++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index d956eb7bf7f..49a91448325 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -200,7 +200,10 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac using parameters::get_parameter; using parameters::choose_parameter; - CGAL_precondition(does_self_intersect(cc)); + CC_VPM cc_vpm = choose_parameter(get_parameter(np_cc, internal_np::vertex_point), + get_property_map(CGAL::vertex_point, cc)); + + CGAL_precondition(does_self_intersect(cc, parameters::vertex_point_map(cc_vpm))); #ifdef CGAL_DEBUG_CLIPPING std::cout << "CC self-intersects" << std::endl; @@ -230,9 +233,6 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac // 2. ------------------------------ // Switch to an exact VPM - CC_VPM cc_vpm = choose_parameter(get_parameter(np_cc, internal_np::vertex_point), - get_property_map(CGAL::vertex_point, cc)); - C2E to_exact; EVPM cc_evpm = get(EP_property_tag(), cc); for(vertex_descriptor vd : vertices(cc)) @@ -246,14 +246,14 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac si_face_evpms.reserve(independent_faces_n); std::size_t counter = 0; - for(face_descriptor fd : si_faces) + for(const face_descriptor fd : si_faces) { #ifdef CGAL_DEBUG_CLIPPING std::cout << " Build single mesh for face: " << fd << " (" << counter+1 << "/" << independent_faces_n << ")" << std::endl; #endif const halfedge_descriptor hd = halfedge(fd, cc); - if(!is_degenerate_triangle_face(face(hd, cc), cc)) + if(!is_degenerate_triangle_face(face(hd, cc), cc, np_cc)) // no need to use cc_evpm here { ++counter; @@ -265,6 +265,12 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac fill_triangle_mesh(hd, cc, cc_vpm, cc_evpm, si_face, si_face_evpm); } +#ifdef CGAL_DEBUG_CLIPPING + else + { + std::cout << "degenerate face..." << std::endl; + } +#endif } #ifdef CGAL_DEBUG_CLIPPING std::cout << clipped_si_faces.size() << " problematic faces to clip" << std::endl; @@ -306,7 +312,8 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac } CGAL_postcondition(is_valid(cc)); - CGAL_postcondition(!does_self_intersect(cc)); + CGAL_postcondition(!does_self_intersect(cc, parameters::vertex_point_map(cc_evpm))); + #ifdef CGAL_DEBUG_CLIPPING std::cout << "Pruned mesh: " << vertices(cc).size() << " nv " << faces(cc).size() << " nf " << std::endl; output_mesh("pruned_mesh.off", cc); From 24e99636f4b579ba449d1ff8e5daded073f3e278 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 25 May 2020 13:18:47 +0200 Subject: [PATCH 09/27] Fix reordered CC IDs usage --- .../internal/Corefinement/clip_experimental.h | 34 +++++++++++-------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 49a91448325..150655f024f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -452,26 +452,30 @@ bool clip_self_intersecting_mesh(TriangleMesh& tm, for(face_descriptor f : faces(tm)) number_of_faces_in_cc[get(pidmap, f)] += 1; - std::vector ordered_cc_ids(num_cc); + std::vector ordered_cc_ids(num_cc); std::iota(std::begin(ordered_cc_ids), std::end(ordered_cc_ids), 0); std::sort(std::begin(ordered_cc_ids), std::end(ordered_cc_ids), [&number_of_faces_in_cc](const int i, const int j) -> bool { return (number_of_faces_in_cc[i] > number_of_faces_in_cc[j]); }); +#ifdef CGAL_DEBUG_CLIPPING for(std::size_t i=0; i ccs(num_cc - 1); // largest CC stays in 'tm' - std::vector cc_vpms(num_cc - 1); + faces_size_type num_cc_to_clip = num_cc - 1; + std::vector ccs(num_cc_to_clip); // largest CC stays in 'tm' + std::vector cc_vpms(num_cc_to_clip); - // Extract all but the largest CC - for(faces_size_type i=1; i tm_cc(tm, ordered_cc_ids[i], pidmap); - CGAL::copy_face_graph(tm_cc, ccs[ordered_cc_ids[i]], - np_tm, parameters::vertex_point_map(cc_vpms[i])); + CGAL::Face_filtered_graph tm_cc(tm, cc_id, pidmap); + CGAL::copy_face_graph(tm_cc, ccs[i], np_tm, parameters::vertex_point_map(cc_vpms[i])); } bool res = true; @@ -487,17 +491,19 @@ bool clip_self_intersecting_mesh(TriangleMesh& tm, res = false; } - for(faces_size_type i=1; i Date: Tue, 26 May 2020 09:40:15 +0200 Subject: [PATCH 10/27] Fix add_new_node overloads not being properly templated since VPM1 might be different from VPM2 and add_new_node will be called sometimes with VPM1/VPM2 and sometimes with VPM2/VPM1 --- .../internal/Corefinement/intersection_nodes.h | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h index 84ff44e7c8b..908afd21c7c 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h @@ -101,12 +101,13 @@ public: //add a new node in the final graph. //it is the intersection of the triangle with the segment + template // VertexPointMap1 or VertexPointMap2 void add_new_node(halfedge_descriptor h_a, face_descriptor f_b, const TriangleMesh& tm_a, const TriangleMesh& tm_b, - const VertexPointMap1& vpm_a, - const VertexPointMap2& vpm_b) + const VPM_A& vpm_a, + const VPM_B& vpm_b) { halfedge_descriptor h_b = halfedge(f_b, tm_b); add_new_node( @@ -222,7 +223,7 @@ public: //add a new node in the final graph. //it is the intersection of the triangle with the segment - template // VertexPointMap1 or VertexPointMap2 + template // VertexPointMap1 or VertexPointMap2 void add_new_node(halfedge_descriptor h_a, face_descriptor f_b, const TriangleMesh& tm_a, @@ -397,12 +398,13 @@ public: //add a new node in the final graph. //it is the intersection of the triangle with the segment + template // VertexPointMap1 or VertexPointMap2 void add_new_node(halfedge_descriptor h_a, face_descriptor f_b, const TriangleMesh& tm_a, const TriangleMesh& tm_b, - const VertexPointMap1& vpm_a, - const VertexPointMap2& vpm_b) + const VPM_A& vpm_a, + const VPM_B& vpm_b) { halfedge_descriptor h_b=halfedge(f_b,tm_b); From f9ad4adfda5a7dcfd10ee2dbed8373f30f0dd412 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 26 May 2020 10:12:32 +0200 Subject: [PATCH 11/27] Use BGL API --- .../internal/Corefinement/clip_experimental.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 150655f024f..830c5d843e6 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -112,7 +112,7 @@ bool clip_mesh_exactly(TriangleMesh& cc, if(!valid_input) { std::cerr << "Invalid input for clip()" << std::endl; - std::cerr << "is cc valid: " << cc.is_valid() << std::endl; + std::cerr << "is cc valid: " << is_valid(cc) << std::endl; std::cerr << "is clipper valid: " << clipper.is_valid() << std::endl; std::cerr << "does part self intersect? " << does_self_intersect(cc, parameters::vertex_point_map(cc_evpm)) << std::endl; std::cerr << "does clipper self intersect? " << does_self_intersect(clipper, parameters::vertex_point_map(clipper_evpm)) << std::endl; @@ -293,8 +293,8 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac continue; for(halfedge_descriptor hd : CGAL::halfedges_around_target(nm_hd, cc)) - if(cc.face(hd) != face_descriptor()) - nm_faces.insert(cc.face(hd)); + if(!is_border(hd, cc)) + nm_faces.insert(face(hd, cc)); } clipped_si_faces.reserve(clipped_si_faces.size() + nm_faces.size()); From 0c1ec4b2c31cd270b18fcb30da5320e86acc5530 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 26 May 2020 10:12:54 +0200 Subject: [PATCH 12/27] Use a dynamic pmap rather than get(face_patch_id(), mesh) to keep backward compatibility A proper fix would be to have a pmap(Polyhedron (without features), face_patch_id) that creates a dynamic pmap, but... --- .../internal/Corefinement/clip_experimental.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 830c5d843e6..37e90a64bcc 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -428,8 +428,8 @@ bool clip_self_intersecting_mesh(TriangleMesh& tm, typedef typename boost::property_map::type DVPM; typedef typename boost::graph_traits::faces_size_type faces_size_type; - typedef typename boost::property_map >::type PatchIDMap; + typedef CGAL::dynamic_face_property_t Patch_ID; + typedef typename boost::property_map::type PatchIDMap; using parameters::get_parameter; using parameters::choose_parameter; @@ -440,7 +440,7 @@ bool clip_self_intersecting_mesh(TriangleMesh& tm, // Splitting connected components avoids having to treat a lot of "easy" self-intersections // (for example from two spheres intersecting) - PatchIDMap pidmap = get(CGAL::face_patch_id_t(), tm); + PatchIDMap pidmap = get(Patch_ID(), tm); faces_size_type num_cc = connected_components(tm, pidmap); #ifdef CGAL_DEBUG_CLIPPING From 808c93c0fade56ad4a033f2817f8610ee564ee7f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 26 May 2020 10:14:06 +0200 Subject: [PATCH 13/27] Fix enforcing that VPM reference types must be equal, value types are sufficient --- .../intersect_triangle_and_segment_3.h | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersect_triangle_and_segment_3.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersect_triangle_and_segment_3.h index cffced176cf..47be440cd3a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersect_triangle_and_segment_3.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersect_triangle_and_segment_3.h @@ -93,22 +93,20 @@ intersection_type( typedef boost::graph_traits GT; typedef typename GT::halfedge_descriptor halfedge_descriptor; typedef std::tuple result_type; - typedef typename boost::property_traits::reference Point_ref; + typedef typename boost::property_traits::reference Point_ref1; + typedef typename boost::property_traits::reference Point_ref2; typedef typename boost::property_traits::value_type Point_3; typedef typename Kernel_traits::Kernel Kernel; - CGAL_static_assertion((std::is_same::reference, - typename boost::property_traits::reference>::value)); - CGAL_static_assertion((std::is_same::value_type, - typename boost::property_traits::value_type>::value)); + CGAL_static_assertion((std::is_same::value_type>::value)); halfedge_descriptor h_2=halfedge(f_2,tm2); - Point_ref a = get(vpm2, target(h_2,tm2) ); - Point_ref b = get(vpm2, target(next(h_2,tm2),tm2) ); - Point_ref c = get(vpm2, source(h_2,tm2) ); - Point_ref p = get(vpm1, source(h_1,tm1) ); - Point_ref q = get(vpm1, target(h_1,tm1) ); + Point_ref2 a = get(vpm2, target(h_2,tm2) ); + Point_ref2 b = get(vpm2, target(next(h_2,tm2),tm2) ); + Point_ref2 c = get(vpm2, source(h_2,tm2) ); + Point_ref1 p = get(vpm1, source(h_1,tm1) ); + Point_ref1 q = get(vpm1, target(h_1,tm1) ); const Orientation abcp = orientation(a,b,c,p); const Orientation abcq = orientation(a,b,c,q); From b444d7ad0d2bd1c53c7d63841071f412877e7c51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 26 May 2020 15:22:43 +0200 Subject: [PATCH 14/27] remove unused overload --- .../internal/Corefinement/intersection_nodes.h | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h index 908afd21c7c..99ff0a7815f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h @@ -119,11 +119,6 @@ public: to_exact( get(vpm_a, target(h_a,tm_a)) ) ) ); } - void add_new_node(halfedge_descriptor edge_1, face_descriptor face_2) - { - add_new_node(edge_1, face_2, tm1, tm2, vpm1, vpm2); - } - template // VertexPointMap1 or VertexPointMap2 void call_put(const VPM& vpm, vertex_descriptor vd, std::size_t i, TriangleMesh&) { @@ -271,11 +266,6 @@ public: add_new_node(*pt); } - void add_new_node(halfedge_descriptor edge_1, face_descriptor face_2) - { - add_new_node(edge_1, face_2, tm1, tm2, vpm1, vpm2); - } - //the point is an input void add_new_node(const Point_3& p){ enodes.push_back(to_exact(p)); @@ -417,12 +407,6 @@ public: get(vpm_a, target(h_a,tm_a)) ) ); } - void add_new_node(halfedge_descriptor edge_1, face_descriptor face_2) - { - add_new_node(edge_1, face_2, tm1, tm2, vpm1, vpm2); - } - - void add_new_node(const Point_3& p) { nodes.push_back(p); From 7e21188463ae62fac491e5938c7a62da7ec19703 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 26 May 2020 15:23:03 +0200 Subject: [PATCH 15/27] comment unused variable @MaelRL was it expected? --- .../internal/Corefinement/clip_experimental.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 37e90a64bcc..81f7db68b41 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -144,7 +144,7 @@ bool clip_mesh_exactly(TriangleMesh& cc, template void fill_triangle_mesh(typename boost::graph_traits::halfedge_descriptor hd, const TriangleMesh& tm, - VPM tm_vpm, + VPM /* tm_vpm */, EVPM tm_evpm, TriangleMesh& si_face, EVPM si_face_evpm) From 27d4c76da35e5080e298559f81de485759d2f7a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 26 May 2020 15:52:00 +0200 Subject: [PATCH 16/27] Remove unused parameter --- .../internal/Corefinement/clip_experimental.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 81f7db68b41..55385b16903 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -141,10 +141,9 @@ bool clip_mesh_exactly(TriangleMesh& cc, return res; } -template +template void fill_triangle_mesh(typename boost::graph_traits::halfedge_descriptor hd, const TriangleMesh& tm, - VPM /* tm_vpm */, EVPM tm_evpm, TriangleMesh& si_face, EVPM si_face_evpm) @@ -263,7 +262,7 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac EVPM& si_face_evpm = si_face_evpms.back(); si_face_evpm = get(EP_property_tag(), si_face); - fill_triangle_mesh(hd, cc, cc_vpm, cc_evpm, si_face, si_face_evpm); + fill_triangle_mesh(hd, cc, cc_evpm, si_face, si_face_evpm); } #ifdef CGAL_DEBUG_CLIPPING else @@ -307,7 +306,7 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac EVPM si_face_evpm = si_face_evpms.back(); const halfedge_descriptor hd = halfedge(fd, cc); - fill_triangle_mesh(hd, cc, cc_vpm, cc_evpm, si_face, si_face_evpm); + fill_triangle_mesh(hd, cc, cc_evpm, si_face, si_face_evpm); CGAL::Euler::remove_face(hd, cc); } From 7ff74ebc891d5e1f83e06896343718f7f55a7e38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 11 Jun 2020 13:37:54 +0200 Subject: [PATCH 17/27] Fix complexity in case of append https://github.com/CGAL/cgal/pull/4778 --- .../CGAL/boost/graph/copy_face_graph.h | 38 +++++++++++-------- 1 file changed, 23 insertions(+), 15 deletions(-) diff --git a/BGL/include/CGAL/boost/graph/copy_face_graph.h b/BGL/include/CGAL/boost/graph/copy_face_graph.h index dcd8c5cedbf..d63d7f61b5a 100644 --- a/BGL/include/CGAL/boost/graph/copy_face_graph.h +++ b/BGL/include/CGAL/boost/graph/copy_face_graph.h @@ -40,7 +40,6 @@ void copy_face_graph_impl(const SourceMesh& sm, TargetMesh& tm, { typedef typename boost::graph_traits::vertex_descriptor sm_vertex_descriptor; typedef typename boost::graph_traits::vertex_descriptor tm_vertex_descriptor; - typedef typename boost::graph_traits::halfedge_iterator tm_halfedge_iterator; typedef typename boost::graph_traits::face_descriptor sm_face_descriptor; typedef typename boost::graph_traits::face_descriptor tm_face_descriptor; @@ -64,14 +63,17 @@ void copy_face_graph_impl(const SourceMesh& sm, TargetMesh& tm, const tm_face_descriptor tm_null_face = boost::graph_traits::null_face(); const tm_vertex_descriptor tm_null_vertex = boost::graph_traits::null_vertex(); - reserve(tm, static_cast::vertices_size_type>(vertices(sm).size()), - static_cast::edges_size_type>(edges(sm).size()), - static_cast::faces_size_type>(faces(sm).size()) ); + reserve(tm, static_cast::vertices_size_type>(vertices(tm).size()+vertices(sm).size()), + static_cast::edges_size_type>(edges(tm).size()+edges(sm).size()), + static_cast::faces_size_type>(faces(tm).size()+faces(sm).size()) ); //insert halfedges and create each vertex when encountering its halfedge + std::vector new_edges; + new_edges.reserve(edges(sm).size()); for(sm_edge_descriptor sm_e : edges(sm)) { tm_edge_descriptor tm_e = add_edge(tm); + new_edges.push_back(tm_e); sm_halfedge_descriptor sm_h = halfedge(sm_e, sm), sm_h_opp = opposite(sm_h, sm); tm_halfedge_descriptor tm_h = halfedge(tm_e, tm), tm_h_opp = opposite(tm_h, tm); @@ -173,9 +175,10 @@ void copy_face_graph_impl(const SourceMesh& sm, TargetMesh& tm, } // detect if there are some non-manifold umbrellas and fix missing halfedge target pointers - for (tm_halfedge_iterator it=halfedges(tm).first; it!=halfedges(tm).second; ++it) + typedef typename std::vector::iterator edge_iterator; + for (edge_iterator it=new_edges.begin(); it!=new_edges.end(); ++it) { - if (target(*it, tm) == tm_null_vertex) + if (target(*it, tm) == tm_null_vertex || source(*it, tm) == tm_null_vertex) { // create and fill a map from target halfedge to source halfedge typedef CGAL::dynamic_halfedge_property_t Dyn_th_tag; @@ -183,17 +186,22 @@ void copy_face_graph_impl(const SourceMesh& sm, TargetMesh& tm, for (sm_halfedge_descriptor hs : halfedges(sm)) put(ht_to_hs, get(hs_to_ht, hs), hs); - for(; it!=halfedges(tm).second; ++it) + for(; it!=new_edges.end(); ++it) { - if (target(*it, tm) == tm_null_vertex) + tm_halfedge_descriptor nh_t = halfedge(*it, tm); + for (int i=0; i<2; ++i) { - // we recover tm_v using the halfedge associated to the target vertex of - // the halfedge in sm corresponding to *it. This is working because we - // set the vertex halfedge pointer to the "same" halfedges. - tm_vertex_descriptor tm_v = - target( get(hs_to_ht, halfedge(target(get(ht_to_hs, *it), sm), sm)), tm); - for(tm_halfedge_descriptor ht : halfedges_around_target(*it, tm)) - set_target(ht, tm_v, tm); + if (target(nh_t, tm) == tm_null_vertex) + { + // we recover tm_v using the halfedge associated to the target vertex of + // the halfedge in sm corresponding to nh_t. This is working because we + // set the vertex halfedge pointer to the "same" halfedges. + tm_vertex_descriptor tm_v = + target( get(hs_to_ht, halfedge(target(get(ht_to_hs, nh_t), sm), sm)), tm); + for(tm_halfedge_descriptor ht : halfedges_around_target(nh_t, tm)) + set_target(ht, tm_v, tm); + } + nh_t = opposite(nh_t, tm); } } break; From 593e8c9d5a68efb83bc5eadfda84e0b2c5cd4fa6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 12 Jun 2020 09:15:37 +0200 Subject: [PATCH 18/27] Debug code (mostly tmp) --- .../internal/Corefinement/clip_experimental.h | 53 ++++++++++++++++++- 1 file changed, 52 insertions(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 55385b16903..67e3e982664 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -119,6 +119,10 @@ bool clip_mesh_exactly(TriangleMesh& cc, std::cerr << "clipper bounds a volume? " << does_bound_a_volume(clipper) << std::endl; return false; } + else + { + std::cerr << "Valid input, proceeding with clipping..." << std::endl; + } #endif bool clip_volume = choose_parameter(get_parameter(np_cc, internal_np::clip_volume), true); @@ -135,8 +139,24 @@ bool clip_mesh_exactly(TriangleMesh& cc, .throw_on_self_intersection(true), parameters::vertex_point_map(clipper_evpm)); +#ifdef CGAL_DEBUG_CLIPPING + std::cout << "Done clipping, check output" << std::endl; + + // useless, we don't care about the clipper_vpm yet, but it shows the structure is broken +// typedef typename CGAL::Cartesian_converter E2C; +// E2C from_exact; + +// for(vertex_descriptor v : vertices(cc)) +// { +// const auto& ep = from_exact(get(cc_evpm, v)); +// put(clipper_vpm, v, ep); +// } + + CGAL_postcondition(cc.is_valid()); CGAL_postcondition(is_valid_polygon_mesh(cc)); + CGAL_postcondition(clipper.is_valid()); CGAL_postcondition(is_valid_polygon_mesh(clipper)); +#endif return res; } @@ -310,7 +330,8 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac CGAL::Euler::remove_face(hd, cc); } - CGAL_postcondition(is_valid(cc)); + CGAL_postcondition(is_valid_polygon_mesh(cc)); + CGAL_assertion(cc.is_valid()); CGAL_postcondition(!does_self_intersect(cc, parameters::vertex_point_map(cc_evpm))); #ifdef CGAL_DEBUG_CLIPPING @@ -326,6 +347,9 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac if(!clip_mesh_exactly(cc, cc_evpm, clipper, np_cc, np_c)) res = false; + CGAL_assertion(cc.is_valid()); + CGAL_assertion(is_valid_polygon_mesh(cc)); + #ifdef CGAL_DEBUG_CLIPPING output_mesh("pruned_mesh_clipped.off", cc); #endif @@ -343,9 +367,15 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac std::cout << " Clipper: " << num_vertices(clipper) << " nv " << num_faces(clipper) << " nf" << std::endl; #endif + CGAL_assertion(clipped_si_face.is_valid()); + CGAL_assertion(is_valid_polygon_mesh(clipped_si_face)); + if(!clip_mesh_exactly(clipped_si_face, si_face_evpm, clipper, np_cc, np_c)) res = false; + CGAL_assertion(clipper.is_valid()); + CGAL_assertion(is_valid_polygon_mesh(clipper)); + if(is_empty(clipped_si_face)) continue; @@ -355,6 +385,11 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac std::cout << " copying..." << std::endl; #endif + CGAL_assertion(clipped_si_face.is_valid()); + CGAL_assertion(cc.is_valid()); + CGAL_assertion(is_valid_polygon_mesh(clipped_si_face)); + CGAL_assertion(is_valid_polygon_mesh(cc)); + CGAL::copy_face_graph(clipped_si_face, cc, parameters::vertex_point_map(si_face_evpm), parameters::vertex_point_map(cc_evpm)); @@ -364,8 +399,16 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac std::cout << " stitching..." << std::endl; #endif + CGAL_assertion(clipped_si_face.is_valid()); + CGAL_assertion(cc.is_valid()); + CGAL_assertion(is_valid_polygon_mesh(clipped_si_face)); + CGAL_assertion(is_valid_polygon_mesh(cc)); + stitch_borders(cc, parameters::vertex_point_map(cc_evpm)); // @todo do something local + CGAL_assertion(cc.is_valid()); + CGAL_assertion(is_valid_polygon_mesh(cc)); + #ifdef CGAL_DEBUG_CLIPPING std::cout << "\t" << vertices(cc).size() << " nv " << faces(cc).size() << " nf" << std::endl; #endif @@ -433,6 +476,10 @@ bool clip_self_intersecting_mesh(TriangleMesh& tm, using parameters::get_parameter; using parameters::choose_parameter; + CGAL_precondition(is_valid(tm)); + CGAL_precondition(is_valid(clipper)); + CGAL_precondition(!does_self_intersect(clipper)); + // @todo keep this? if(!does_self_intersect(tm, np_tm)) return clip(tm, clipper, np_tm, np_c); @@ -474,6 +521,10 @@ bool clip_self_intersecting_mesh(TriangleMesh& tm, cc_vpms[i] = get(P_property_tag(), ccs[i]); CGAL::Face_filtered_graph tm_cc(tm, cc_id, pidmap); + + CGAL_assertion(tm_cc.is_selection_valid()); + CGAL_assertion(is_valid(tm_cc)); + CGAL::copy_face_graph(tm_cc, ccs[i], np_tm, parameters::vertex_point_map(cc_vpms[i])); } From 8ce1a9c7da4cef3ab50a72558c153aa1fd25cd06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 12 Jun 2020 09:15:52 +0200 Subject: [PATCH 19/27] Duplicate possible non-manifold vertices in the input --- .../internal/Corefinement/clip_experimental.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 67e3e982664..331c5f155c9 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -484,6 +484,8 @@ bool clip_self_intersecting_mesh(TriangleMesh& tm, if(!does_self_intersect(tm, np_tm)) return clip(tm, clipper, np_tm, np_c); + duplicate_non_manifold_vertices(tm, np_tm); + // Splitting connected components avoids having to treat a lot of "easy" self-intersections // (for example from two spheres intersecting) PatchIDMap pidmap = get(Patch_ID(), tm); From e7fd86c2efd0973a80a8b8aef7c7c69e0c0af969 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 12 Jun 2020 09:16:41 +0200 Subject: [PATCH 20/27] Add example of clipping w/ SI --- .../Polygon_mesh_processing/CMakeLists.txt | 1 + .../clip_with_self_intersections.cpp | 132 ++++++++++++++++++ 2 files changed, 133 insertions(+) create mode 100644 Polygon_mesh_processing/examples/Polygon_mesh_processing/clip_with_self_intersections.cpp diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 46a7f04d886..e47fb3bbcdc 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -90,6 +90,7 @@ create_single_source_cgal_program( "repair_polygon_soup_example.cpp" ) create_single_source_cgal_program( "mesh_smoothing_example.cpp") create_single_source_cgal_program( "locate_example.cpp") create_single_source_cgal_program( "orientation_pipeline_example.cpp") +create_single_source_cgal_program( "clip_with_self_intersections.cpp") #create_single_source_cgal_program( "self_snapping_example.cpp") #create_single_source_cgal_program( "snapping_example.cpp") diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/clip_with_self_intersections.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/clip_with_self_intersections.cpp new file mode 100644 index 00000000000..6029561ea79 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/clip_with_self_intersections.cpp @@ -0,0 +1,132 @@ +#define CGAL_DEBUG_CLIPPING + +#include +#include +#include + +#include + +#include +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +// typedef CGAL::Exact_predicates_exact_constructions_kernel K; + +typedef CGAL::Surface_mesh Mesh; + +typedef K::Point_3 Point_3; +typedef K::Plane_3 Plane_3; + +namespace PMP = CGAL::Polygon_mesh_processing; + +template +bool read_mesh(const char* filename, Mesh& sm) +{ + typedef typename K::Point_3 Point; + + std::ifstream in(filename, std::ios::binary); + if(!in.good()) + { + std::cerr << "Error: can't read file: " << filename << std::endl; + std::exit(1); + } + + std::string fn(filename); + if(fn.substr(fn.find_last_of(".") + 1) == "stl") + { + std::vector points; + std::vector > faces; + if(!CGAL::read_STL(in, points, faces)) + { + std::cerr << "Error: cannot read STL mesh\n"; + return false; + } + + std::cout << "Cleaning polygon soup..." << std::endl; + PMP::repair_polygon_soup(points, faces); + + if(!CGAL::Polygon_mesh_processing::orient_polygon_soup(points, faces)) + std::cerr << "W: File does not describe a polygon mesh" << std::endl; + + CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, faces, sm); + } + else if(fn.substr(fn.find_last_of(".") + 1) == "off") + { + if(!(in >> sm)) + { + std::cerr << "Cannot read .OFF as polygon mesh, reading .OFF as a polygon soup...\n"; + + in.clear(); // clear fail and eof bits + in.seekg(0, std::ios::beg); // back at the start + + std::vector points; + std::vector > faces; + if(!CGAL::read_OFF(in, points, faces)) + { + std::cerr << "Error: cannot read OFF mesh as a polygon soup\n"; + return false; + } + + std::cout << "Cleaning polygon soup..." << std::endl; + PMP::repair_polygon_soup(points, faces); + + if(!CGAL::Polygon_mesh_processing::orient_polygon_soup(points, faces)) + std::cerr << "W: File does not describe a polygon mesh" << std::endl; + + CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, faces, sm); + } + } + else + { + std::cerr << "Unknown file type" << std::endl; + return false; + } + + if(!CGAL::is_triangle_mesh(sm)) + { + std::cerr << "Input geometry is not triangulated." << std::endl; + return false; + } + + return true; +} + +int main(int argc, char** argv) +{ + const char* filename = (argc > 1) ? argv[1] : "data/mannequin-devil.off"; + Mesh mesh; + if(!read_mesh(filename, mesh)) + return EXIT_FAILURE; + + std::cout << "Input mesh: " << vertices(mesh).size() << " nv " << faces(mesh).size() << " nf " << std::endl; + +#if 0//def CGAL_USE_CLIPPER_PLANE + Plane_3 pl(Point_3(0,0,-600), Point_3(1,0,-600), Point_3(0,1,-600)); + PMP::clip(mesh, pl, + CGAL::parameters::geom_traits(K()) + .input_might_have_self_intersection(true)); +#else + filename = (argc > 2) ? argv[2] : "data/sphere_clipper.off"; + Mesh clipper; + if(!read_mesh(filename, clipper)) + return EXIT_FAILURE; + + std::cout << "Clipper mesh: " << num_vertices(clipper) << " nv " << num_faces(clipper) << " nf" << std::endl; + + if(PMP::is_outward_oriented(clipper)) + PMP::reverse_face_orientations(clipper); + + PMP::clip(mesh, clipper, CGAL::parameters::input_might_have_self_intersection(true)); +#endif + + std::ofstream out("clipped.off"); + out.precision(17); + out << mesh; + + std::cout << "Clipped!" << std::endl; + + return EXIT_SUCCESS; +} From 00d25379751be07e285d49ea2a66f9ea6c182aee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 12 Jun 2020 10:37:34 +0200 Subject: [PATCH 21/27] Don't abuse Surface_mesh property maps --- .../internal/Corefinement/clip_experimental.h | 35 +++++++++++-------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index 331c5f155c9..b978e440b0a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -65,7 +65,7 @@ template bool clip_mesh_exactly(TriangleMesh& cc, EVPM cc_evpm, - TriangleMesh clipper, // intentional copy + TriangleMesh& clipper, const NamedParameters1& np_cc, const NamedParameters2& np_c) { @@ -99,24 +99,29 @@ bool clip_mesh_exactly(TriangleMesh& cc, Clipper_VPM clipper_vpm = choose_parameter(get_parameter(np_c, internal_np::vertex_point), get_property_map(CGAL::vertex_point, clipper)); + std::unordered_map v2v; + TriangleMesh clipper_copy; + CGAL::copy_face_graph(clipper, clipper_copy, + parameters::vertex_to_vertex_output_iterator(std::inserter(v2v, v2v.end()))); + C2E to_exact; - EVPM clipper_evpm = get(EP_property_tag(), clipper); - for(vertex_descriptor vd : vertices(clipper)) - put(clipper_evpm, vd, to_exact(get(clipper_vpm, vd))); + EVPM clipper_copy_evpm = get(EP_property_tag(), clipper_copy); + for(vertex_descriptor v : vertices(clipper)) + put(clipper_copy_evpm, v2v[v], to_exact(get(clipper_vpm, v))); #ifdef CGAL_DEBUG_CLIPPING - const bool valid_input = is_valid(cc) && is_valid(clipper) && + const bool valid_input = is_valid_polygon_mesh(cc) && is_valid_polygon_mesh(clipper_copy) && !does_self_intersect(cc, parameters::vertex_point_map(cc_evpm)) && - !does_self_intersect(clipper, parameters::vertex_point_map(clipper_evpm)) && - does_bound_a_volume(clipper, parameters::vertex_point_map(clipper_evpm)); + !does_self_intersect(clipper_copy, parameters::vertex_point_map(clipper_copy_evpm)) && + does_bound_a_volume(clipper_copy, parameters::vertex_point_map(clipper_copy_evpm)); if(!valid_input) { std::cerr << "Invalid input for clip()" << std::endl; std::cerr << "is cc valid: " << is_valid(cc) << std::endl; - std::cerr << "is clipper valid: " << clipper.is_valid() << std::endl; + std::cerr << "is clipper valid: " << clipper_copy.is_valid() << std::endl; std::cerr << "does part self intersect? " << does_self_intersect(cc, parameters::vertex_point_map(cc_evpm)) << std::endl; - std::cerr << "does clipper self intersect? " << does_self_intersect(clipper, parameters::vertex_point_map(clipper_evpm)) << std::endl; - std::cerr << "clipper bounds a volume? " << does_bound_a_volume(clipper) << std::endl; + std::cerr << "does clipper self intersect? " << does_self_intersect(clipper_copy, parameters::vertex_point_map(clipper_copy_evpm)) << std::endl; + std::cerr << "clipper bounds a volume? " << does_bound_a_volume(clipper_copy) << std::endl; return false; } else @@ -132,12 +137,12 @@ bool clip_mesh_exactly(TriangleMesh& cc, clip_volume = false; use_compact_clipper = false; - bool res = clip(cc, clipper, + bool res = clip(cc, clipper_copy, parameters::vertex_point_map(cc_evpm) .clip_volume(clip_volume) .use_compact_clipper(use_compact_clipper) .throw_on_self_intersection(true), - parameters::vertex_point_map(clipper_evpm)); + parameters::vertex_point_map(clipper_copy_evpm)); #ifdef CGAL_DEBUG_CLIPPING std::cout << "Done clipping, check output" << std::endl; @@ -154,8 +159,8 @@ bool clip_mesh_exactly(TriangleMesh& cc, CGAL_postcondition(cc.is_valid()); CGAL_postcondition(is_valid_polygon_mesh(cc)); - CGAL_postcondition(clipper.is_valid()); - CGAL_postcondition(is_valid_polygon_mesh(clipper)); + CGAL_postcondition(clipper_copy.is_valid()); + CGAL_postcondition(is_valid_polygon_mesh(clipper_copy)); #endif return res; @@ -521,7 +526,7 @@ bool clip_self_intersecting_mesh(TriangleMesh& tm, faces_size_type cc_id = ordered_cc_ids[1+i]; // extract all but the largest CC, so start at '1' CGAL_assertion(cc_id < num_cc); - cc_vpms[i] = get(P_property_tag(), ccs[i]); + cc_vpms[i] = get(P_property_tag(), ccs[i]); // ccs only contains the extracted CCs, so it starts at 0 CGAL::Face_filtered_graph tm_cc(tm, cc_id, pidmap); CGAL_assertion(tm_cc.is_selection_valid()); From 784d0cddcd9fec0622e09c30d0d8d6f891c3a037 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 12 Jun 2020 13:26:34 +0200 Subject: [PATCH 22/27] also add non-manifold border faces --- .../internal/Corefinement/clip_experimental.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h index b978e440b0a..3f5607b42e7 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h @@ -301,6 +301,8 @@ bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_fac #endif // 3. ------------------------------ + CGAL::expand_face_selection_for_removal(si_faces, cc, CGAL::make_boolean_property_map(si_faces)); + for(const face_descriptor fd : si_faces) CGAL::Euler::remove_face(halfedge(fd, cc), cc); From 7e2e9d849acd2a866ae6976152f9be1fd73bb371 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 3 Jul 2020 15:02:00 +0200 Subject: [PATCH 23/27] Remove experimental code related to clipping self-intersecting meshes See functionality added in https://github.com/CGAL/cgal/pull/4790 instead. Keeping it in a commit instead of filtering the branch to keep the code in history. --- .../clip_with_self_intersections.cpp | 132 ---- .../CGAL/Polygon_mesh_processing/clip.h | 37 +- .../internal/Corefinement/clip_experimental.h | 591 ------------------ 3 files changed, 10 insertions(+), 750 deletions(-) delete mode 100644 Polygon_mesh_processing/examples/Polygon_mesh_processing/clip_with_self_intersections.cpp delete mode 100644 Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/clip_with_self_intersections.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/clip_with_self_intersections.cpp deleted file mode 100644 index 6029561ea79..00000000000 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/clip_with_self_intersections.cpp +++ /dev/null @@ -1,132 +0,0 @@ -#define CGAL_DEBUG_CLIPPING - -#include -#include -#include - -#include - -#include -#include - -#include -#include - -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -// typedef CGAL::Exact_predicates_exact_constructions_kernel K; - -typedef CGAL::Surface_mesh Mesh; - -typedef K::Point_3 Point_3; -typedef K::Plane_3 Plane_3; - -namespace PMP = CGAL::Polygon_mesh_processing; - -template -bool read_mesh(const char* filename, Mesh& sm) -{ - typedef typename K::Point_3 Point; - - std::ifstream in(filename, std::ios::binary); - if(!in.good()) - { - std::cerr << "Error: can't read file: " << filename << std::endl; - std::exit(1); - } - - std::string fn(filename); - if(fn.substr(fn.find_last_of(".") + 1) == "stl") - { - std::vector points; - std::vector > faces; - if(!CGAL::read_STL(in, points, faces)) - { - std::cerr << "Error: cannot read STL mesh\n"; - return false; - } - - std::cout << "Cleaning polygon soup..." << std::endl; - PMP::repair_polygon_soup(points, faces); - - if(!CGAL::Polygon_mesh_processing::orient_polygon_soup(points, faces)) - std::cerr << "W: File does not describe a polygon mesh" << std::endl; - - CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, faces, sm); - } - else if(fn.substr(fn.find_last_of(".") + 1) == "off") - { - if(!(in >> sm)) - { - std::cerr << "Cannot read .OFF as polygon mesh, reading .OFF as a polygon soup...\n"; - - in.clear(); // clear fail and eof bits - in.seekg(0, std::ios::beg); // back at the start - - std::vector points; - std::vector > faces; - if(!CGAL::read_OFF(in, points, faces)) - { - std::cerr << "Error: cannot read OFF mesh as a polygon soup\n"; - return false; - } - - std::cout << "Cleaning polygon soup..." << std::endl; - PMP::repair_polygon_soup(points, faces); - - if(!CGAL::Polygon_mesh_processing::orient_polygon_soup(points, faces)) - std::cerr << "W: File does not describe a polygon mesh" << std::endl; - - CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, faces, sm); - } - } - else - { - std::cerr << "Unknown file type" << std::endl; - return false; - } - - if(!CGAL::is_triangle_mesh(sm)) - { - std::cerr << "Input geometry is not triangulated." << std::endl; - return false; - } - - return true; -} - -int main(int argc, char** argv) -{ - const char* filename = (argc > 1) ? argv[1] : "data/mannequin-devil.off"; - Mesh mesh; - if(!read_mesh(filename, mesh)) - return EXIT_FAILURE; - - std::cout << "Input mesh: " << vertices(mesh).size() << " nv " << faces(mesh).size() << " nf " << std::endl; - -#if 0//def CGAL_USE_CLIPPER_PLANE - Plane_3 pl(Point_3(0,0,-600), Point_3(1,0,-600), Point_3(0,1,-600)); - PMP::clip(mesh, pl, - CGAL::parameters::geom_traits(K()) - .input_might_have_self_intersection(true)); -#else - filename = (argc > 2) ? argv[2] : "data/sphere_clipper.off"; - Mesh clipper; - if(!read_mesh(filename, clipper)) - return EXIT_FAILURE; - - std::cout << "Clipper mesh: " << num_vertices(clipper) << " nv " << num_faces(clipper) << " nf" << std::endl; - - if(PMP::is_outward_oriented(clipper)) - PMP::reverse_face_orientations(clipper); - - PMP::clip(mesh, clipper, CGAL::parameters::input_might_have_self_intersection(true)); -#endif - - std::ofstream out("clipped.off"); - out.precision(17); - out << mesh; - - std::cout << "Clipped!" << std::endl; - - return EXIT_SUCCESS; -} diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h index 3bdcde94489..d059c3eb33d 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h @@ -15,8 +15,6 @@ #include -#include - #include #include #include @@ -36,9 +34,11 @@ #include -namespace CGAL { +namespace CGAL{ namespace Polygon_mesh_processing { -namespace internal { + +namespace internal +{ template int @@ -402,17 +402,6 @@ void split_along_edges(TriangleMesh& tm, } // end of internal namespace -namespace experimental { - -template -bool clip_self_intersecting_mesh(TriangleMesh& tm, - TriangleMesh& clipper, - const NamedParameters1& np_tm, - const NamedParameters2& np_c); - -} // namespace experimental - /** * \ingroup PMP_corefinement_grp * clips `tm` by keeping the part that is inside the volume \link coref_def_subsec bounded \endlink @@ -472,19 +461,14 @@ clip(TriangleMesh& tm, const NamedParameters1& np_tm, const NamedParameters2& np_c) { - using parameters::choose_parameter; - using parameters::get_parameter; - - const bool clip_with_si = choose_parameter(get_parameter(np_tm, internal_np::input_might_have_self_intersection), false); - if(clip_with_si) - return experimental::clip_self_intersecting_mesh(tm, clipper, np_tm, np_c); - - const bool clip_volume = choose_parameter(get_parameter(np_tm, internal_np::clip_volume), false); + const bool clip_volume = + parameters::choose_parameter(parameters::get_parameter(np_tm, internal_np::clip_volume), false); if(clip_volume && is_closed(tm)) - return corefine_and_compute_intersection(tm, clipper, tm, np_tm, np_c, np_tm); + return corefine_and_compute_intersection(tm, clipper, tm, np_tm, np_c); return corefine_and_compute_intersection(tm, clipper, tm, - np_tm.use_bool_op_to_clip_surface(true), np_c, np_tm); + np_tm.use_bool_op_to_clip_surface(true), + np_c); } /** @@ -876,7 +860,6 @@ void split(TriangleMesh& tm, /// \endcond -} // namespace Polygon_mesh_processing -} // namespace CGAL +} } //end of namespace CGAL::Polygon_mesh_processing #endif // CGAL_POLYGON_MESH_PROCESSING_CLIP_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h deleted file mode 100644 index 3f5607b42e7..00000000000 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/clip_experimental.h +++ /dev/null @@ -1,591 +0,0 @@ -// Copyright (c) 2019-2020 GeometryFactory (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// -// Author(s) : Mael Rouxel-Labbé -// Sebastien Loriot - -#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CLIP_EXPERIMENTAL_H -#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CLIP_EXPERIMENTAL_H - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -namespace CGAL { -namespace Polygon_mesh_processing { - -template -bool clip(TriangleMesh& tm, - TriangleMesh& clipper, - const NamedParameters1& np_tm, - const NamedParameters2& np_c); - -namespace internal { - -template -void output_mesh(const char* filename, const PolygonMesh& pmesh) -{ - std::ofstream output(filename); - output.precision(17); - output << pmesh; - output.close(); -} - -template -bool clip_mesh_exactly(TriangleMesh& cc, - EVPM cc_evpm, - TriangleMesh& clipper, - const NamedParameters1& np_cc, - const NamedParameters2& np_c) -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - // These typedefs concern the VPM of the CLIPPER mesh - typedef typename GetK::Kernel Clipper_kernel; - typedef typename GetVertexPointMap::type Clipper_VPM; - -#if 1 - typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; - typedef typename Exact_kernel::Point_3 Exact_point; - - typedef typename CGAL::Cartesian_converter C2E; -#else - typedef typename CGAL::Exact_kernel_selector::Exact_kernel Exact_kernel; - typedef typename Exact_kernel::Point_3 Exact_point; - - typedef typename CGAL::Exact_kernel_selector::C2E C2E; -#endif - - typedef CGAL::dynamic_vertex_property_t EP_property_tag; - - // must have equal exact point types because of corefinement.h - CGAL_static_assertion((std::is_same::type>::value)); - - using parameters::get_parameter; - using parameters::choose_parameter; - - // Initialize an exact vpm for the clipper - Clipper_VPM clipper_vpm = choose_parameter(get_parameter(np_c, internal_np::vertex_point), - get_property_map(CGAL::vertex_point, clipper)); - - std::unordered_map v2v; - TriangleMesh clipper_copy; - CGAL::copy_face_graph(clipper, clipper_copy, - parameters::vertex_to_vertex_output_iterator(std::inserter(v2v, v2v.end()))); - - C2E to_exact; - EVPM clipper_copy_evpm = get(EP_property_tag(), clipper_copy); - for(vertex_descriptor v : vertices(clipper)) - put(clipper_copy_evpm, v2v[v], to_exact(get(clipper_vpm, v))); - -#ifdef CGAL_DEBUG_CLIPPING - const bool valid_input = is_valid_polygon_mesh(cc) && is_valid_polygon_mesh(clipper_copy) && - !does_self_intersect(cc, parameters::vertex_point_map(cc_evpm)) && - !does_self_intersect(clipper_copy, parameters::vertex_point_map(clipper_copy_evpm)) && - does_bound_a_volume(clipper_copy, parameters::vertex_point_map(clipper_copy_evpm)); - if(!valid_input) - { - std::cerr << "Invalid input for clip()" << std::endl; - std::cerr << "is cc valid: " << is_valid(cc) << std::endl; - std::cerr << "is clipper valid: " << clipper_copy.is_valid() << std::endl; - std::cerr << "does part self intersect? " << does_self_intersect(cc, parameters::vertex_point_map(cc_evpm)) << std::endl; - std::cerr << "does clipper self intersect? " << does_self_intersect(clipper_copy, parameters::vertex_point_map(clipper_copy_evpm)) << std::endl; - std::cerr << "clipper bounds a volume? " << does_bound_a_volume(clipper_copy) << std::endl; - return false; - } - else - { - std::cerr << "Valid input, proceeding with clipping..." << std::endl; - } -#endif - - bool clip_volume = choose_parameter(get_parameter(np_cc, internal_np::clip_volume), true); - bool use_compact_clipper = choose_parameter(get_parameter(np_cc, internal_np::use_compact_clipper), false); - - // @todo is it possible to forward clip_volume/use_compact_clipper/visitor? - clip_volume = false; - use_compact_clipper = false; - - bool res = clip(cc, clipper_copy, - parameters::vertex_point_map(cc_evpm) - .clip_volume(clip_volume) - .use_compact_clipper(use_compact_clipper) - .throw_on_self_intersection(true), - parameters::vertex_point_map(clipper_copy_evpm)); - -#ifdef CGAL_DEBUG_CLIPPING - std::cout << "Done clipping, check output" << std::endl; - - // useless, we don't care about the clipper_vpm yet, but it shows the structure is broken -// typedef typename CGAL::Cartesian_converter E2C; -// E2C from_exact; - -// for(vertex_descriptor v : vertices(cc)) -// { -// const auto& ep = from_exact(get(cc_evpm, v)); -// put(clipper_vpm, v, ep); -// } - - CGAL_postcondition(cc.is_valid()); - CGAL_postcondition(is_valid_polygon_mesh(cc)); - CGAL_postcondition(clipper_copy.is_valid()); - CGAL_postcondition(is_valid_polygon_mesh(clipper_copy)); -#endif - - return res; -} - -template -void fill_triangle_mesh(typename boost::graph_traits::halfedge_descriptor hd, - const TriangleMesh& tm, - EVPM tm_evpm, - TriangleMesh& si_face, - EVPM si_face_evpm) -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - vertex_descriptor vd0 = add_vertex(si_face); - vertex_descriptor vd1 = add_vertex(si_face); - vertex_descriptor vd2 = add_vertex(si_face); - - // Note that we don't even fill the internal VPM - put(si_face_evpm, vd0, get(tm_evpm, source(hd, tm))); - put(si_face_evpm, vd1, get(tm_evpm, target(hd, tm))); - put(si_face_evpm, vd2, get(tm_evpm, target(next(hd, tm), tm))); - - CGAL::Euler::add_face(std::initializer_list{vd0, vd1, vd2}, si_face); -} - -template -bool clip_single_self_intersecting_cc(const FacePairRange& self_intersecting_faces, - TriangleMesh& cc, - TriangleMesh& clipper, - const NamedParameters1& np_cc, - const NamedParameters2& np_c) -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - typedef typename GetK::Kernel CC_kernel; - typedef typename GetVertexPointMap::type CC_VPM; - -#if 1 - typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; - typedef typename Exact_kernel::Point_3 Exact_point; - - typedef typename CGAL::Cartesian_converter C2E; - typedef typename CGAL::Cartesian_converter E2C; -#else - typedef typename CGAL::Exact_kernel_selector::Exact_kernel Exact_kernel; - typedef typename Exact_kernel::Point_3 Exact_point; - - typedef typename CGAL::Exact_kernel_selector::C2E C2E; - typedef typename CGAL::Exact_kernel_selector::E2C E2C; -#endif - - typedef CGAL::dynamic_vertex_property_t EP_property_tag; - typedef typename boost::property_map::type EVPM; - - typedef std::pair Pair_of_faces; - - using parameters::get_parameter; - using parameters::choose_parameter; - - CC_VPM cc_vpm = choose_parameter(get_parameter(np_cc, internal_np::vertex_point), - get_property_map(CGAL::vertex_point, cc)); - - CGAL_precondition(does_self_intersect(cc, parameters::vertex_point_map(cc_vpm))); - -#ifdef CGAL_DEBUG_CLIPPING - std::cout << "CC self-intersects" << std::endl; -#endif - - bool res = true; - - // 1. Gather the problematic faces - // 2. Copy (one-by-one) the problematic faces into independent meshes - // 3. Remove the problematic faces from 'cc' - // 4. Clip (cc - problematic_faces) - // 2bis. Clip the problematic faces - // 5. Re-add the clipped problematic faces from step 2. - - // 1. ------------------------------ - std::set si_faces; - for(const Pair_of_faces& fp : self_intersecting_faces) - { - si_faces.insert(fp.first); - si_faces.insert(fp.second); - } - - std::size_t independent_faces_n = si_faces.size(); -#ifdef CGAL_DEBUG_CLIPPING - std::cout << independent_faces_n << " independent faces" << std::endl; -#endif - - // 2. ------------------------------ - // Switch to an exact VPM - C2E to_exact; - EVPM cc_evpm = get(EP_property_tag(), cc); - for(vertex_descriptor vd : vertices(cc)) - put(cc_evpm, vd, to_exact(get(cc_vpm, vd))); - - // forced to do that because clip() wants both meshes to be of the same type (rather than just requiring - // that the two VPMs have matching types) - std::vector clipped_si_faces; - clipped_si_faces.reserve(independent_faces_n); - std::vector si_face_evpms; - si_face_evpms.reserve(independent_faces_n); - - std::size_t counter = 0; - for(const face_descriptor fd : si_faces) - { -#ifdef CGAL_DEBUG_CLIPPING - std::cout << " Build single mesh for face: " << fd << " (" << counter+1 << "/" << independent_faces_n << ")" << std::endl; -#endif - - const halfedge_descriptor hd = halfedge(fd, cc); - if(!is_degenerate_triangle_face(face(hd, cc), cc, np_cc)) // no need to use cc_evpm here - { - ++counter; - - clipped_si_faces.resize(counter); - TriangleMesh& si_face = clipped_si_faces.back(); - si_face_evpms.resize(counter); - EVPM& si_face_evpm = si_face_evpms.back(); - si_face_evpm = get(EP_property_tag(), si_face); - - fill_triangle_mesh(hd, cc, cc_evpm, si_face, si_face_evpm); - } -#ifdef CGAL_DEBUG_CLIPPING - else - { - std::cout << "degenerate face..." << std::endl; - } -#endif - } -#ifdef CGAL_DEBUG_CLIPPING - std::cout << clipped_si_faces.size() << " problematic faces to clip" << std::endl; -#endif - - // 3. ------------------------------ - CGAL::expand_face_selection_for_removal(si_faces, cc, CGAL::make_boolean_property_map(si_faces)); - - for(const face_descriptor fd : si_faces) - CGAL::Euler::remove_face(halfedge(fd, cc), cc); - - // 3bis. iteratively remove faces incident to non-manifold vertices - std::vector non_manifold_cones; - CGAL::Polygon_mesh_processing::non_manifold_vertices(cc, std::back_inserter(non_manifold_cones)); - std::set nm_vertices; - - std::unordered_set nm_faces; - for(halfedge_descriptor nm_hd : non_manifold_cones) - { - // keep the faces incident to the nm vertex for only one (the first one) of the cones - if(nm_vertices.insert(target(nm_hd, cc)).second) - continue; - - for(halfedge_descriptor hd : CGAL::halfedges_around_target(nm_hd, cc)) - if(!is_border(hd, cc)) - nm_faces.insert(face(hd, cc)); - } - - clipped_si_faces.reserve(clipped_si_faces.size() + nm_faces.size()); - for(face_descriptor fd : nm_faces) - { - clipped_si_faces.resize(clipped_si_faces.size() + 1); - si_face_evpms.resize(si_face_evpms.size() + 1); - - TriangleMesh& si_face = clipped_si_faces.back(); - EVPM si_face_evpm = si_face_evpms.back(); - - const halfedge_descriptor hd = halfedge(fd, cc); - fill_triangle_mesh(hd, cc, cc_evpm, si_face, si_face_evpm); - CGAL::Euler::remove_face(hd, cc); - } - - CGAL_postcondition(is_valid_polygon_mesh(cc)); - CGAL_assertion(cc.is_valid()); - CGAL_postcondition(!does_self_intersect(cc, parameters::vertex_point_map(cc_evpm))); - -#ifdef CGAL_DEBUG_CLIPPING - std::cout << "Pruned mesh: " << vertices(cc).size() << " nv " << faces(cc).size() << " nf " << std::endl; - output_mesh("pruned_mesh.off", cc); -#endif - - // 4. ------------------------------ -#ifdef CGAL_DEBUG_CLIPPING - std::cout << "Clipping CC's main part (w/o self-intersecting faces)" << std::endl; -#endif - - if(!clip_mesh_exactly(cc, cc_evpm, clipper, np_cc, np_c)) - res = false; - - CGAL_assertion(cc.is_valid()); - CGAL_assertion(is_valid_polygon_mesh(cc)); - -#ifdef CGAL_DEBUG_CLIPPING - output_mesh("pruned_mesh_clipped.off", cc); -#endif - - // 5. ------------------------------ - counter = 0; - for(std::size_t i=0, csfs=clipped_si_faces.size(); i -bool clip_single_cc(TriangleMesh& cc, - TriangleMesh& clipper, - const NamedParameters1& np_cc, - const NamedParameters2& np_c) -{ - typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef std::pair Pair_of_faces; - -#ifdef CGAL_DEBUG_CLIPPING - std::cout << "Clip single CC of size " << vertices(cc).size() << " nv " << faces(cc).size() << " nf" << std::endl; -#endif - - std::vector intersecting_faces; - self_intersections(cc, std::back_inserter(intersecting_faces), np_cc); - - if(intersecting_faces.empty()) - return clip(cc, clipper, np_cc, np_c); - else - return clip_single_self_intersecting_cc(intersecting_faces, cc, clipper, np_cc, np_c); -} - -} // namespace internal - -namespace experimental { - -// Try to split the mesh in multiple connected components and clip them independently -template -bool clip_self_intersecting_mesh(TriangleMesh& tm, - TriangleMesh& clipper, - const NamedParameters1& np_tm, - const NamedParameters2& np_c) -{ - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - typedef typename GetVertexPointMap::type VPM; - typedef typename boost::property_traits::value_type Point; - - typedef CGAL::dynamic_vertex_property_t P_property_tag; - typedef typename boost::property_map::type DVPM; - - typedef typename boost::graph_traits::faces_size_type faces_size_type; - typedef CGAL::dynamic_face_property_t Patch_ID; - typedef typename boost::property_map::type PatchIDMap; - - using parameters::get_parameter; - using parameters::choose_parameter; - - CGAL_precondition(is_valid(tm)); - CGAL_precondition(is_valid(clipper)); - CGAL_precondition(!does_self_intersect(clipper)); - - // @todo keep this? - if(!does_self_intersect(tm, np_tm)) - return clip(tm, clipper, np_tm, np_c); - - duplicate_non_manifold_vertices(tm, np_tm); - - // Splitting connected components avoids having to treat a lot of "easy" self-intersections - // (for example from two spheres intersecting) - PatchIDMap pidmap = get(Patch_ID(), tm); - faces_size_type num_cc = connected_components(tm, pidmap); - -#ifdef CGAL_DEBUG_CLIPPING - std::cout << num_cc << " connected component(s)" << std::endl; -#endif - - // order ccs by size - std::vector number_of_faces_in_cc(num_cc, 0); - for(face_descriptor f : faces(tm)) - number_of_faces_in_cc[get(pidmap, f)] += 1; - - std::vector ordered_cc_ids(num_cc); - std::iota(std::begin(ordered_cc_ids), std::end(ordered_cc_ids), 0); - std::sort(std::begin(ordered_cc_ids), std::end(ordered_cc_ids), - [&number_of_faces_in_cc](const int i, const int j) -> bool { - return (number_of_faces_in_cc[i] > number_of_faces_in_cc[j]); - }); - -#ifdef CGAL_DEBUG_CLIPPING - for(std::size_t i=0; i ccs(num_cc_to_clip); // largest CC stays in 'tm' - std::vector cc_vpms(num_cc_to_clip); - - for(faces_size_type i=0; i tm_cc(tm, cc_id, pidmap); - - CGAL_assertion(tm_cc.is_selection_valid()); - CGAL_assertion(is_valid(tm_cc)); - - CGAL::copy_face_graph(tm_cc, ccs[i], np_tm, parameters::vertex_point_map(cc_vpms[i])); - } - - bool res = true; - const CGAL::Bbox_3 clipper_bbox = CGAL::Polygon_mesh_processing::bbox(clipper, np_c); - - // Clip CC by CC - keep_connected_components(tm, std::vector(1, ordered_cc_ids[0]), pidmap); - - const CGAL::Bbox_3 tm_bbox = CGAL::Polygon_mesh_processing::bbox(tm, np_tm); - if(CGAL::do_overlap(tm_bbox, clipper_bbox)) - { - if(!internal::clip_single_cc(tm, clipper, np_tm, np_c)) - res = false; - } - - for(faces_size_type i=0; i -bool clip_self_intersecting_mesh(TriangleMesh& tm, TriangleMesh& clipper, const NamedParameters& np_tm) -{ - return clip_self_intersecting_mesh(tm, clipper, np_tm, parameters::all_default()); -} - -template -bool clip_self_intersecting_mesh(TriangleMesh& tm, TriangleMesh& clipper) -{ - return clip_self_intersecting_mesh(tm, clipper, parameters::all_default(), parameters::all_default()); -} - -} // namespace experimental -} // namespace Polygon_mesh_processing -} // namespace CGAL - - -#endif //CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CLIP_EXPERIMENTAL_H From bbedb5e8b40a9bc4a780c1f8ff7259c714b84290 Mon Sep 17 00:00:00 2001 From: Mael Date: Mon, 13 Jul 2020 10:19:06 +0200 Subject: [PATCH 24/27] Fix name of property maps --- .../corefinement_consecutive_bool_op.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp index 7b13f109a6d..6ade126528c 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp @@ -83,10 +83,10 @@ int main(int argc, char* argv[]) } Exact_point_map mesh1_exact_points = - mesh1.add_property_map("e:exact_point").first; + mesh1.add_property_map("v:exact_point").first; Exact_point_map mesh2_exact_points = - mesh2.add_property_map("e:exact_point").first; + mesh2.add_property_map("v:exact_point").first; Exact_vertex_point_map mesh1_vpm(mesh1_exact_points, mesh1); Exact_vertex_point_map mesh2_vpm(mesh2_exact_points, mesh2); From 6553b0dfaa24f5ab4471c20ee7b8f26af37c70c7 Mon Sep 17 00:00:00 2001 From: Mael Date: Thu, 16 Jul 2020 09:27:55 +0200 Subject: [PATCH 25/27] Remove deleted example from CMakeLists.txt --- .../examples/Polygon_mesh_processing/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index e47fb3bbcdc..46a7f04d886 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -90,7 +90,6 @@ create_single_source_cgal_program( "repair_polygon_soup_example.cpp" ) create_single_source_cgal_program( "mesh_smoothing_example.cpp") create_single_source_cgal_program( "locate_example.cpp") create_single_source_cgal_program( "orientation_pipeline_example.cpp") -create_single_source_cgal_program( "clip_with_self_intersections.cpp") #create_single_source_cgal_program( "self_snapping_example.cpp") #create_single_source_cgal_program( "snapping_example.cpp") From cec2440297595e407e419e2ae1cd9717cc8a77d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 22 Sep 2020 11:27:18 +0200 Subject: [PATCH 26/27] Remove obsolete named parameter --- BGL/include/CGAL/boost/graph/parameters_interface.h | 1 - 1 file changed, 1 deletion(-) diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h index 32f4b8b8b16..86ce60eca70 100644 --- a/BGL/include/CGAL/boost/graph/parameters_interface.h +++ b/BGL/include/CGAL/boost/graph/parameters_interface.h @@ -76,7 +76,6 @@ CGAL_add_named_parameter(projection_functor_t, projection_functor, projection_fu CGAL_add_named_parameter(throw_on_self_intersection_t, throw_on_self_intersection, throw_on_self_intersection) CGAL_add_named_parameter(clip_volume_t, clip_volume, clip_volume) CGAL_add_named_parameter(use_compact_clipper_t, use_compact_clipper, use_compact_clipper) -CGAL_add_named_parameter(input_might_have_self_intersection_t, input_might_have_self_intersection, input_might_have_self_intersection) CGAL_add_named_parameter(output_iterator_t, output_iterator, output_iterator) CGAL_add_named_parameter(erase_all_duplicates_t, erase_all_duplicates, erase_all_duplicates) CGAL_add_named_parameter(require_same_orientation_t, require_same_orientation, require_same_orientation) From 4d4eef94a40cdbe1371bce570a96f3f34d723fd7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 28 Sep 2020 15:44:29 +0200 Subject: [PATCH 27/27] Enable different VPMs in PMP::clip --- .../CGAL/Polygon_mesh_processing/clip.h | 20 ++++++++++++------- .../Generic_clip_output_builder.h | 15 +++++++------- 2 files changed, 21 insertions(+), 14 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h index ae94faba8e2..311528101f8 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h @@ -33,7 +33,14 @@ #include +#include +#include +#include +#include +#include #include +#include +#include namespace CGAL{ namespace Polygon_mesh_processing { @@ -420,10 +427,9 @@ generic_clip_impl( NamedParameters1>::type Vpm; typedef typename GetVertexPointMap::type Vpm2; - CGAL_USE_TYPE(Vpm2); - CGAL_assertion_code( - static const bool same_vpm = (boost::is_same::value); ) - CGAL_static_assertion(same_vpm); + + CGAL_static_assertion((std::is_same::value_type, + typename boost::property_traits::value_type>::value)); Vpm vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point), get_property_map(boost::vertex_point, tm1)); @@ -477,17 +483,17 @@ generic_clip_impl( // surface intersection algorithm call typedef Corefinement::Generic_clip_output_builder Ob; typedef Corefinement::Surface_intersection_visitor_for_corefinement< - TriangleMesh, Vpm, Ob, Ecm_in, User_visitor> Algo_visitor; + TriangleMesh, Vpm, Vpm2, Ob, Ecm_in, User_visitor> Algo_visitor; Ecm_in ecm_in(tm1,tm2,ecm1,ecm2); Ob ob(tm1, tm2, vpm1, vpm2, algo_ecm1, fid_map1, use_compact_clipper); - Corefinement::Intersection_of_triangle_meshes + Corefinement::Intersection_of_triangle_meshes functor(tm1, tm2, vpm1, vpm2, Algo_visitor(uv,ob,ecm_in,&tm2)); functor(CGAL::Emptyset_iterator(), false, true); } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Generic_clip_output_builder.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Generic_clip_output_builder.h index bcc8e75469d..d84788ff868 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Generic_clip_output_builder.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Generic_clip_output_builder.h @@ -38,7 +38,8 @@ namespace PMP=Polygon_mesh_processing; namespace params=PMP::parameters; template @@ -48,7 +49,7 @@ class Generic_clip_output_builder typedef typename Default::Get< Kernel_, typename Kernel_traits< - typename boost::property_traits::value_type + typename boost::property_traits::value_type >::Kernel >::type Kernel; // graph_traits typedefs @@ -76,8 +77,8 @@ class Generic_clip_output_builder //Data members TriangleMesh &tm1, &tm2; // property maps of input meshes - const VertexPointMap vpm1; - const VertexPointMap vpm2; + const VertexPointMap1 vpm1; + const VertexPointMap2 vpm2; Ecm1 ecm1; FaceIdMap1 fids1; bool use_compact_clipper; @@ -105,8 +106,8 @@ public: Generic_clip_output_builder(TriangleMesh& tm1, TriangleMesh& tm2, - const VertexPointMap vpm1, - const VertexPointMap vpm2, + const VertexPointMap1 vpm1, + const VertexPointMap2 vpm2, const Ecm1& ecm1, FaceIdMap1 fids1, bool use_compact_clipper) @@ -167,7 +168,7 @@ public: typedef Side_of_triangle_mesh Inside_poly_test; + VertexPointMap2> Inside_poly_test; CGAL::Bounded_side in_tm2 = is_tm2_inside_out ? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE;