diff --git a/Box_intersection_d/include/CGAL/box_intersection_d.h b/Box_intersection_d/include/CGAL/box_intersection_d.h index 2e8d103e2c9c..b689ca9f56cb 100644 --- a/Box_intersection_d/include/CGAL/box_intersection_d.h +++ b/Box_intersection_d/include/CGAL/box_intersection_d.h @@ -65,6 +65,10 @@ void box_intersection_segment_tree_d( static_assert (!std::is_convertible::value, "Parallel_tag is enabled but TBB is unavailable."); #else // CGAL_LINKED_WITH_TBB + // if(std::is_convertible::value) + // { + // std::nth_element() + // } if(std::is_convertible::value) { // Here is an illustration for n=2. diff --git a/PMP_Boolean_operations/benchmark/CMakeLists.txt b/PMP_Boolean_operations/benchmark/CMakeLists.txt index 1cd65d7416b7..566cbfa4705e 100644 --- a/PMP_Boolean_operations/benchmark/CMakeLists.txt +++ b/PMP_Boolean_operations/benchmark/CMakeLists.txt @@ -9,6 +9,7 @@ find_package(CGAL REQUIRED) create_single_source_cgal_program("rotated_cubes_autorefinement.cpp") create_single_source_cgal_program("coplanar_cubes_autorefinement.cpp") +create_single_source_cgal_program("corefinement_mesh_union_bench.cpp") create_single_source_cgal_program("Performance/performance_snap_polygon_soup.cpp") create_single_source_cgal_program("Robustness/robustness_snap_polygon_soup.cpp") @@ -22,6 +23,7 @@ include(CGAL_TBB_support) if(TARGET CGAL::TBB_support) target_link_libraries(rotated_cubes_autorefinement PRIVATE CGAL::TBB_support) target_link_libraries(coplanar_cubes_autorefinement PRIVATE CGAL::TBB_support) + target_link_libraries(corefinement_mesh_union_bench PRIVATE CGAL::TBB_support) else() message(STATUS "NOTICE: Intel TBB was not found. Sequential code will be used.") endif() diff --git a/PMP_Boolean_operations/benchmark/corefinement_mesh_union_bench.cpp b/PMP_Boolean_operations/benchmark/corefinement_mesh_union_bench.cpp new file mode 100644 index 000000000000..c6f0381b1ec9 --- /dev/null +++ b/PMP_Boolean_operations/benchmark/corefinement_mesh_union_bench.cpp @@ -0,0 +1,235 @@ +#define CGAL_OUTPUT_BUILDER_RUNNING_TIME + +#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; + +namespace PMP = CGAL::Polygon_mesh_processing; + +struct Visitor_rep{ + + Visitor_rep(double normalize = 4) + : normalize(normalize) + { + t.start(); + } + + void progress_filtering_intersections(double d) + { + // d /= normalize; + // total += d; + // if(total > bound){ + // // std::cout << std::setprecision(3) << total*100 << " % in " << std::setprecision(5) << t.time() << " sec." << std::endl; + // bound += 0.1; + // } + } + + void start_triangulating_faces(std::size_t tf) + { + // tfaces = tf; + // bound_faces = tf/10; + } + + void face_triangulation(std::size_t i) + { + // if(i> bound_faces){ + // // std::cout << double(i)/double(tfaces) * 100 << " %" << std::endl; + // bound_faces += tfaces/10; + // } + } + + void start_coplanar_faces(std::size_t tc) + { + // std::cout << "Visitor::start_coplanar_faces() at " << t.time() << " sec." << std::endl; + tcoplanar= tc; + count_coplanar = 0; + bound_coplanar = tcoplanar/10; + } + + void intersection_of_coplanar_faces_step() + { + ++count_coplanar; + if(count_coplanar> bound_coplanar){ + // std::cout << "Visitor::coplanar_faces: " << double(count_coplanar)/double(tcoplanar) * 100 << " % " << std::endl; + bound_coplanar += tcoplanar/10; + } + } + + void start_intersection_points(std::size_t ti) + { + // std::cout << "Visitor::start_intersection_points() at " << t.time() << " sec." << std::endl; + tintersection= ti; + count_intersection = 0; + bound_intersection = tintersection/10; + } + + void edge_face_intersections_step() + { + ++count_intersection; + if(count_intersection> bound_intersection){ + // std::cout << "Visitor::intersection_points: " << double(count_intersection)/double(tintersection) * 100 << " % " << std::endl; + bound_intersection += tintersection/10; + } + } + + double time() const + { + return t.time(); + } + + void reset_timer(){ t.stop(); t.reset(); t.start(); } + + double normalize; + double bound = 0.1; + double total = 0; + std::size_t count = 0; + + std::size_t bound_faces = 0; + std::size_t tfaces = 0; + + std::size_t bound_coplanar = 0; + std::size_t tcoplanar = 0; + std::size_t count_coplanar = 0; + + std::size_t bound_intersection = 0; + std::size_t tintersection = 0; + std::size_t count_intersection = 0; + CGAL::Real_timer t; + CGAL::Real_timer local; +}; + + +struct Visitor : + public PMP::Corefinement::Default_visitor +{ + std::shared_ptr sptr; + mutable std::size_t tf_counter = 0; + + Visitor() + : sptr(std::make_shared()) + {} + + void progress_filtering_intersections(double d) + { + sptr->progress_filtering_intersections(d); + } + + void start_filtering_intersections() const + { + // std::cout << "Visitor::start_filtering_intersections() at " << sptr->time() << " sec." << std::endl; + sptr->reset_timer(); + } + void end_filtering_intersections() const + { + std::cout << "Filtering_intersections: " << sptr->time() << std::endl; + } + + void start_triangulating_faces(std::size_t tf) const + { + std::cout << "Triangulate " << tf << " faces in "; + sptr->start_triangulating_faces(tf); + tf_counter = 0; + sptr->reset_timer(); + } + + void triangulating_faces_step() const + { + // sptr->face_triangulation(tf_counter++); + } + + void end_triangulating_faces()const + { + std::cout << sptr->time() << " sec." << std::endl; + } + + void start_handling_intersection_of_coplanar_faces(std::size_t i) const + { + sptr->start_coplanar_faces(i); + sptr->reset_timer(); + } + + void intersection_of_coplanar_faces_step() const + { + sptr->intersection_of_coplanar_faces_step(); + } + + void end_handling_intersection_of_coplanar_faces() const + { + std::cout << "Handling coplanar in " << sptr->time() << " sec." << std::endl; + } + + void start_handling_edge_face_intersections(std::size_t i) const + { + sptr->start_intersection_points(i); + sptr->reset_timer(); + } + + void edge_face_intersections_step() const + { + sptr->edge_face_intersections_step(); + } + + void end_handling_edge_face_intersections() const + { + std::cout << "Compute intersection points in " << sptr->time() << " sec." << std::endl; + } + + void start_building_output() const + { + sptr->reset_timer(); + } + + void end_building_output() const + { + std::cout << "Build the output in " << sptr->time() << " sec." << std::endl; + } +}; + + +int main(int argc, char* argv[]) +{ + const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); + const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off"); + + Mesh mesh1, mesh2; + if(!CGAL::IO::read_polygon_mesh(filename1, mesh1) || !CGAL::IO::read_polygon_mesh(filename2, mesh2)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + CGAL::Real_timer rt; + CGAL::Timer t; + rt.start(); t.start(); + Mesh out; + Visitor visitor; + + bool valid_union = PMP::corefine_and_compute_union (mesh1, mesh2, out, CGAL::parameters::visitor(visitor).concurrency_tag(CGAL::Parallel_tag())); + // bool valid_union = PMP::corefine_and_compute_difference (mesh1, mesh2, out, CGAL::parameters::visitor(visitor).concurrency_tag(CGAL::Parallel_tag())); + + std::cout << "Global timer = " << rt.time() << " sec." << " ( " << t.time() << " cpu time)" << std::endl; + + + if(valid_union) + { + std::cout << "Union was successfully computed\n"; + CGAL::IO::write_polygon_mesh("union.off", out, CGAL::parameters::stream_precision(17)); + return 0; + } + + std::cout << "Union could not be computed\n"; + + return 1; +} diff --git a/PMP_Boolean_operations/examples/PMP_Boolean_operations/CMakeLists.txt b/PMP_Boolean_operations/examples/PMP_Boolean_operations/CMakeLists.txt index ae035b0858de..2eb75498e21c 100644 --- a/PMP_Boolean_operations/examples/PMP_Boolean_operations/CMakeLists.txt +++ b/PMP_Boolean_operations/examples/PMP_Boolean_operations/CMakeLists.txt @@ -53,6 +53,9 @@ if(TARGET CGAL::TBB_support) create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp") target_link_libraries(corefinement_parallel_union_meshes PRIVATE CGAL::TBB_support) + + target_link_libraries(corefinement_mesh_union_progress PRIVATE CGAL::TBB_support) + target_link_libraries(corefinement_mesh_union PRIVATE CGAL::TBB_support) else() message(STATUS "NOTICE: Intel TBB was not found. Sequential code will be used.") endif() diff --git a/PMP_Boolean_operations/examples/PMP_Boolean_operations/corefinement_mesh_union.cpp b/PMP_Boolean_operations/examples/PMP_Boolean_operations/corefinement_mesh_union.cpp index bf226b0dc0c0..f6c501116007 100644 --- a/PMP_Boolean_operations/examples/PMP_Boolean_operations/corefinement_mesh_union.cpp +++ b/PMP_Boolean_operations/examples/PMP_Boolean_operations/corefinement_mesh_union.cpp @@ -1,6 +1,8 @@ #include #include +#include + #include #include @@ -25,7 +27,11 @@ int main(int argc, char* argv[]) } Mesh out; - bool valid_union = PMP::corefine_and_compute_union(mesh1,mesh2, out); + CGAL::Real_timer rt; CGAL::Timer t; + rt.start(); t.start(); + // bool valid_union = PMP::corefine_and_compute_union(mesh1,mesh2, out); + bool valid_union = PMP::corefine_and_compute_union(mesh1,mesh2, out, CGAL::parameters::concurrency_tag(CGAL::Parallel_tag())); + std::cout << "run in " << rt.time() << " (" << t.time() << "s)" << std::endl; if(valid_union) { diff --git a/PMP_Boolean_operations/examples/PMP_Boolean_operations/corefinement_mesh_union_progress.cpp b/PMP_Boolean_operations/examples/PMP_Boolean_operations/corefinement_mesh_union_progress.cpp index 67c2c3615fe5..6d7bf3abfcac 100644 --- a/PMP_Boolean_operations/examples/PMP_Boolean_operations/corefinement_mesh_union_progress.cpp +++ b/PMP_Boolean_operations/examples/PMP_Boolean_operations/corefinement_mesh_union_progress.cpp @@ -4,6 +4,8 @@ #include #include +#include + #include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; @@ -17,6 +19,7 @@ struct Visitor_rep{ : normalize(normalize) { t.start(); + rt.start(); } void progress_filtering_intersections(double d) @@ -24,7 +27,7 @@ struct Visitor_rep{ d /= normalize; total += d; if(total > bound){ - std::cout << std::setprecision(3) << total*100 << " % in " << std::setprecision(5) << t.time() << " sec." << std::endl; + std::cout << std::setprecision(3) << total*100 << " % in " << std::setprecision(5) << rt.time() << " sec (" << t.time() << "s)." << std::endl; bound += 0.1; } } @@ -45,7 +48,7 @@ struct Visitor_rep{ void start_coplanar_faces(std::size_t tc) { - std::cout << "Visitor::start_coplanar_faces() at " << t.time() << " sec." << std::endl; + std::cout << "Visitor::start_coplanar_faces() at " << rt.time() << " sec (" << t.time() << "s)." << std::endl; tcoplanar= tc; count_coplanar = 0; bound_coplanar = tcoplanar/10; @@ -62,7 +65,7 @@ struct Visitor_rep{ void start_intersection_points(std::size_t ti) { - std::cout << "Visitor::start_intersection_points() at " << t.time() << " sec." << std::endl; + std::cout << "Visitor::start_intersection_points() at " << rt.time() << " sec (" << t.time() << "s)." << std::endl; tintersection= ti; count_intersection = 0; bound_intersection = tintersection/10; @@ -82,6 +85,11 @@ struct Visitor_rep{ return t.time(); } + double real_time() const + { + return rt.time(); + } + double normalize; double bound = 0.1; double total = 0; @@ -98,6 +106,7 @@ struct Visitor_rep{ std::size_t tintersection = 0; std::size_t count_intersection = 0; CGAL::Timer t; + CGAL::Real_timer rt; }; @@ -118,16 +127,16 @@ struct Visitor : void start_filtering_intersections() const { - std::cout << "Visitor::start_filtering_intersections() at " << sptr->time() << " sec." << std::endl; + std::cout << "Visitor::start_filtering_intersections() at " << sptr->real_time() << " sec (" << sptr->time() << "s)." << std::endl; } void end_filtering_intersections() const { - std::cout << "Visitor::end_filtering_intersections() at " << sptr->time() << " sec." << std::endl; + std::cout << "Visitor::end_filtering_intersections() at " << sptr->real_time() << " sec (" << sptr->time() << "s)." << std::endl; } void start_triangulating_faces(std::size_t tf) const { - std::cout << "Visitor::start_triangulation() with " << tf << " faces at " << sptr->time() << " sec." << std::endl; + std::cout << "Visitor::start_triangulation() with " << tf << " faces at " << sptr->real_time() << " sec (" << sptr->time() << "s)." << std::endl; sptr->start_triangulating_faces(tf); tf_counter = 0; } @@ -139,7 +148,7 @@ struct Visitor : void end_triangulating_faces()const { - std::cout << "Visitor::end_triangulating_faces() at " << sptr->time() << " sec." << std::endl; + std::cout << "Visitor::end_triangulating_faces() at " << sptr->real_time() << " sec (" << sptr->time() << "s)." << std::endl; } void start_handling_intersection_of_coplanar_faces(std::size_t i) const @@ -154,7 +163,7 @@ struct Visitor : void end_handling_intersection_of_coplanar_faces() const { - std::cout << "Visitor::end_coplanar_faces() at " << sptr->time() << " sec." << std::endl; + std::cout << "Visitor::end_coplanar_faces() at " << sptr->real_time() << " sec (" << sptr->time() << "s)." << std::endl; } void start_handling_edge_face_intersections(std::size_t i) const @@ -169,25 +178,27 @@ struct Visitor : void end_handling_edge_face_intersections() const { - std::cout << "Visitor::end_intersection_points() at " << sptr->time() << " sec." << std::endl; + std::cout << "Visitor::end_intersection_points() at " << sptr->real_time() << " sec (" << sptr->time() << "s)." << std::endl; } void start_building_output() const { - std::cout << "Visitor::start_building_output() at " << sptr->time() << " sec."<< std::endl; + std::cout << "Visitor::start_building_output() at " << sptr->real_time() << " sec (" << sptr->time() << "s)." << std::endl; } void end_building_output() const { - std::cout << "Visitor::end_building_output() at " << sptr->time() << " sec." << std::endl; + std::cout << "Visitor::end_building_output() at " << sptr->real_time() << " sec (" << sptr->time() << "s)." << std::endl; } }; int main(int argc, char* argv[]) { - const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); - const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off"); + // const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); + // const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off"); + const std::string filename1 = (argc > 1) ? argv[1] : "/home/oem/Data/data/iphigenia.off"; + const std::string filename2 = (argc > 2) ? argv[2] : "/home/oem/Data/data/iphigenia_transformed.off"; Mesh mesh1, mesh2; if(!CGAL::IO::read_polygon_mesh(filename1, mesh1) || !CGAL::IO::read_polygon_mesh(filename2, mesh2)) @@ -196,14 +207,16 @@ int main(int argc, char* argv[]) return 1; } + CGAL::Real_timer rt; CGAL::Timer t; - t.start(); + rt.start(); t.start(); Mesh out; Visitor visitor; - bool valid_union = PMP::corefine_and_compute_union (mesh1,mesh2, out, CGAL::parameters::visitor(visitor)); + bool valid_union = PMP::corefine_and_compute_union (mesh1, mesh2, out, CGAL::parameters::visitor(visitor).concurrency_tag(CGAL::Parallel_tag())); + // bool valid_union = PMP::corefine_and_compute_difference (mesh1, mesh2, out, CGAL::parameters::visitor(visitor).concurrency_tag(CGAL::Parallel_tag())); - std::cout << "Global timer = " << t.time() << " sec." << std::endl; + std::cout << "Global timer = " << rt.time() << " sec." << " ( " << t.time() << " cpu time)" << std::endl; if(valid_union) diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/corefinement.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/corefinement.h index 4176b99fbe4e..5cd10b8f1636 100644 --- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/corefinement.h +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/corefinement.h @@ -284,6 +284,12 @@ corefine_and_compute_boolean_operations( using parameters::choose_parameter; using parameters::get_parameter; + using Concurrency_tag = typename internal_np::Lookup_named_param_def < + internal_np::concurrency_tag_t, + NPIn1, + Sequential_tag + > ::type; + const bool throw_on_self_intersection = choose_parameter(get_parameter(np1, internal_np::throw_on_self_intersection), false); @@ -489,7 +495,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); @@ -762,6 +768,12 @@ corefine( TriangleMesh& tm1, using parameters::choose_parameter; using parameters::get_parameter; + using Concurrency_tag = typename internal_np::Lookup_named_param_def < + internal_np::concurrency_tag_t, + NamedParameters1, + Sequential_tag + > ::type; + TriangleMesh* const_mesh_ptr=nullptr; if (choose_parameter(get_parameter(np1, internal_np::do_not_modify), false)) { @@ -835,7 +847,7 @@ corefine( TriangleMesh& tm1, 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,const_mesh_ptr), const_mesh_ptr); // Fill non-manifold feature maps if provided diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h index 85a73ae8628b..4f484448f713 100644 --- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h @@ -22,6 +22,17 @@ #include #include +#include + +namespace CGAL::internal +{ +template +struct is_surface_mesh : std::false_type {}; + +template +struct is_surface_mesh> : std::true_type {}; +} + #include #include @@ -1077,32 +1088,41 @@ class Face_graph_output_builder // (1) Assign a patch id to each facet indicating in which connected // component limited by intersection edges of the surface they are. - // ... for tm1 + std::vector tm1_patch_sizes; + std::vector tm2_patch_sizes; + std::size_t nb_patches_tm1, nb_patches_tm2; std::vector tm1_patch_ids( num_faces(tm1),NID ); Border_edge_map is_marked_1(intersection_edges1, tm1); - std::size_t nb_patches_tm1 = + auto connected_components_1=[&]() + { + // ... for tm1 + nb_patches_tm1 = connected_components(tm1, make_compose_property_map(fids1,make_property_map(&tm1_patch_ids[0])), parameters::edge_is_constrained_map(is_marked_1) .face_index_map(fids1)); - - std::vector tm1_patch_sizes(nb_patches_tm1, 0); + tm1_patch_sizes.resize(nb_patches_tm1, 0); for(std::size_t i : tm1_patch_ids) if(i!=NID) ++tm1_patch_sizes[i]; - // ... for tm2 + }; + std::vector tm2_patch_ids( num_faces(tm2),NID ); Border_edge_map is_marked_2(intersection_edges2, tm2); - std::size_t nb_patches_tm2 = + auto connected_components_2=[&]() + { + // ... for tm2 + nb_patches_tm2 = connected_components(tm2, make_compose_property_map(fids2,make_property_map(&tm2_patch_ids[0])), parameters::edge_is_constrained_map(is_marked_2) .face_index_map(fids2)); - - std::vector tm2_patch_sizes(nb_patches_tm2, 0); + tm2_patch_sizes.resize(nb_patches_tm2, 0); for(Node_id i : tm2_patch_ids) if(i!=NID) ++tm2_patch_sizes[i]; + }; + tbb::parallel_invoke(connected_components_1, connected_components_2); #ifdef CGAL_COREFINEMENT_DEBUG std::cout << "nb_patches_tm1 = " << nb_patches_tm1 << "\n"; @@ -2361,24 +2381,68 @@ class Face_graph_output_builder ); std::vector shared_edges; - - #define CGAL_COREF_FUNCTION_CALL_DEF(BO_type) \ - fill_new_triangle_mesh( \ - output, \ - patches_of_tm1_used[BO_type], patches_of_tm2_used[BO_type], \ - patches_of_tm1, patches_of_tm2, \ - BO_type == TM2_MINUS_TM1, BO_type == TM1_MINUS_TM2, \ - polylines, \ - intersection_edges1, intersection_edges2, \ - vpm1, vpm2, *std::get(output_vpms), \ - marks_on_input_edges.ecm1, \ - marks_on_input_edges.ecm2, \ - std::get(out_edge_mark_maps), \ - shared_edges, \ - user_visitor \ - ) - CGAL_COREF_FUNCTION_CALL(operation) - #undef CGAL_COREF_FUNCTION_CALL_DEF + if constexpr(::CGAL::internal::is_surface_mesh::value){ + #define CGAL_COREF_FUNCTION_CALL_DEF(BO_type) \ + output.template fill_new_triangle_mesh( \ + patches_of_tm1_used[BO_type], patches_of_tm2_used[BO_type], \ + patches_of_tm1, patches_of_tm2, \ + polylines, \ + intersection_edges1, intersection_edges2, \ + vpm1, vpm2, *std::get(output_vpms), \ + marks_on_input_edges.ecm1, \ + marks_on_input_edges.ecm2, \ + std::get(out_edge_mark_maps), \ + shared_edges, \ + user_visitor \ + ) + CGAL_COREF_FUNCTION_CALL(operation) + #undef CGAL_COREF_FUNCTION_CALL_DEF + } + else + { +#if 0 + std::vector< typename VertexPointMap1::value_type > pts; + std::vector< boost::container::small_vector > faces; + #define CGAL_COREF_FUNCTION_CALL_DEF(BO_type) \ + fill_triangle_soup( \ + pts, \ + faces, \ + patches_of_tm1_used[BO_type], patches_of_tm2_used[BO_type], \ + patches_of_tm1, patches_of_tm2, \ + BO_type == TM2_MINUS_TM1, BO_type == TM1_MINUS_TM2, \ + polylines, \ + intersection_edges1, intersection_edges2, \ + vpm1, vpm2, \ + marks_on_input_edges.ecm1, \ + marks_on_input_edges.ecm2, \ + user_visitor \ + ) + CGAL_COREF_FUNCTION_CALL(operation) + #undef CGAL_COREF_FUNCTION_CALL_DEF + #define CGAL_COREF_FUNCTION_CALL_DEF(BO_type) \ + polygon_soup_to_polygon_mesh(pts, faces, output, parameters::default_values(), parameters::vertex_point_map(*std::get(output_vpms))) + CGAL_COREF_FUNCTION_CALL(operation) + #undef CGAL_COREF_FUNCTION_CALL_DEF +#else + #define CGAL_COREF_FUNCTION_CALL_DEF(BO_type) \ + fill_new_triangle_mesh( \ + output, \ + patches_of_tm1_used[BO_type], patches_of_tm2_used[BO_type], \ + patches_of_tm1, patches_of_tm2, \ + BO_type == TM2_MINUS_TM1, BO_type == TM1_MINUS_TM2, \ + polylines, \ + intersection_edges1, intersection_edges2, \ + vpm1, vpm2, *std::get(output_vpms), \ + marks_on_input_edges.ecm1, \ + marks_on_input_edges.ecm2, \ + std::get(out_edge_mark_maps), \ + shared_edges, \ + user_visitor \ + ) + CGAL_COREF_FUNCTION_CALL(operation) + #undef CGAL_COREF_FUNCTION_CALL_DEF +#endif + } mark_edges(out_edge_mark_maps, shared_edges, operation); } @@ -2568,10 +2632,10 @@ class Face_graph_output_builder for(halfedge_descriptor h : patches_of_tm1[i].interior_edges) { - if ( !border_vertices.count( target(h,tm1) ) ) - patches_of_tm1[i].interior_vertices.insert( target(h,tm1) ); - if ( !border_vertices.count( source(h,tm1) ) ) - patches_of_tm1[i].interior_vertices.insert( source(h,tm1) ); + if( halfedge(target(h, tm1),tm1) == h && !border_vertices.count( target(h,tm1) )) // We insert a vertex only once + patches_of_tm1[i].interior_vertices.push_back( target(h,tm1) ); + if( halfedge(source(h, tm1),tm1) == opposite(h, tm1) && !border_vertices.count( source(h,tm1) )) // We insert a vertex only once + patches_of_tm1[i].interior_vertices.push_back( source(h,tm1) ); } } } @@ -2597,8 +2661,10 @@ class Face_graph_output_builder if (all_removed) id_p_rm.erase(id_p_rm.begin()); // remove the vertex from the interior vertices of patches to be removed - for(std::size_t pid : id_p_rm) - patches_of_tm1[pid].interior_vertices.erase(vd); + for(std::size_t pid : id_p_rm){ + auto &interior_vertices = patches_of_tm1[pid].interior_vertices; + interior_vertices.erase(std::find(interior_vertices.begin(), interior_vertices.end(), vd)); + } // we now need to update the next/prev relationship induced by the future removal of patches // that will not be updated after patch removal diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h index 3794336dad27..4d0e8ebd9462 100644 --- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h @@ -38,6 +38,10 @@ namespace CGAL { namespace Polygon_mesh_processing { namespace Corefinement{ + + CGAL_MUTEX vertices_mutex; +CGAL_MUTEX triangles_mutex; + // TODO option to ignore internal edges for patches of coplanar faces //binds two edge constrained pmaps @@ -1256,12 +1260,15 @@ class Surface_intersection_visitor_for_corefinement{ //check if one of the triangle input vertex is also a node for (int ik=0;ik<3;++ik) - if ( f_indices[ik]second; for (int i=0;i<3;++i) std::copy(f_boundary.node_ids_array[i].begin(), @@ -1291,6 +1298,7 @@ class Surface_intersection_visitor_for_corefinement{ if (is_face_border) { call_put(marks_on_edges,tm,edge(h,tm),true); + CGAL_SCOPED_LOCK(triangles_mutex); output_builder.set_edge_per_polyline(tm,std::make_pair(id, id_n),h); } else @@ -1302,6 +1310,7 @@ class Surface_intersection_visitor_for_corefinement{ halfedge_descriptor hn=halfedge(vn, tm); while(face(hn, tm) != f) hn=opposite(next(hn, tm), tm); + CGAL_SCOPED_LOCK(triangles_mutex); constraints.emplace_back(make_array(std::make_pair(hi,id),std::make_pair(hn, id_n))); } } @@ -1317,6 +1326,7 @@ class Surface_intersection_visitor_for_corefinement{ std::vector new_faces; for (const std::array, 2>& a : constraints) { + CGAL_SCOPED_LOCK(triangles_mutex); halfedge_descriptor nh = Euler::split_face(a[0].first, a[1].first, tm); new_faces.push_back(face(opposite(nh, tm), tm)); @@ -1327,6 +1337,7 @@ class Surface_intersection_visitor_for_corefinement{ // now triangulate new faces if (!new_faces.empty()) { + CGAL_SCOPED_LOCK(triangles_mutex); new_faces.push_back(f); for(face_descriptor nf : new_faces) { diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h index 4276641ce7c7..cb95d6e0e7ca 100644 --- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h @@ -15,6 +15,7 @@ #include +#include #include #include @@ -23,6 +24,14 @@ #include #include +#include +#include +#include + +#ifdef CGAL_OUTPUT_BUILDER_RUNNING_TIME +#include +#endif + namespace CGAL { namespace Polygon_mesh_processing { namespace Corefinement { @@ -244,6 +253,98 @@ struct Compute_bbox { BPM bpm; }; +template +struct AABB_tree_build_helper +{ + using Primitive = AABB_face_graph_triangle_primitive; + using Traits = AABB_traits_3; + using Tree = AABB_tree; + + using Graph_traits = boost::graph_traits; + using vertex_descriptor = typename Graph_traits::vertex_descriptor; + using face_descriptor = typename Graph_traits::face_descriptor; + + template + struct Split_primitives + { + Split_primitives(const RPM &rpm): rpm(rpm){} + + template + void operator()(PrimitiveIterator first, + PrimitiveIterator beyond, + const CGAL::Bbox_3& bbox) const + { + auto longest_axis=[](const CGAL::Bbox_3& bbox){ + const double dx = bbox.x_span(); + const double dy = bbox.y_span(); + const double dz = bbox.z_span(); + return (dx>=dy) ? ((dx>=dz) ? 0 : 2) : ((dy>=dz) ? 1 : 2); + }; + + PrimitiveIterator middle = first + (beyond - first)/2; + const int crd=longest_axis(bbox); + std::nth_element(first, middle, beyond, [this, crd](const Primitive& p1, const Primitive& p2){ return get(rpm, p1.id())[crd] < get(rpm, p2.id())[crd];}); + } + const RPM &rpm; + }; + + // For exact side_of_triangle_mesh + template + struct Compute_bbox { + Compute_bbox(const BPM& bpm): bpm(bpm){} + + template + CGAL::Bbox_3 operator()(ConstPrimitiveIterator first, + ConstPrimitiveIterator beyond) const + { + CGAL::Bbox_3 bbox = get(bpm, first->id()); + for(++first; first != beyond; ++first) + bbox += get(bpm, first->id()); + return bbox; + } + BPM bpm; + }; + + template + void build(Tree& tree, const TriangleMesh& tm, VertexPointMap& vpm){ + using Face_bbox_tag = typename CGAL::dynamic_face_property_t; + using Face_ref_point_tag = typename CGAL::dynamic_face_property_t; + using Bbox_map = typename boost::property_map::const_type; + using Ref_point_map = typename boost::property_map::const_type; + + using VPM_kernel = typename Kernel_traits::value_type>::Kernel; + CGAL::Cartesian_converter to_input; + +#ifdef CGAL_OUTPUT_BUILDER_RUNNING_TIME + CGAL::Real_timer t; + t.start(); +#endif + Bbox_map bb_map = get(Face_bbox_tag(), tm); + Ref_point_map rp_map = get(Face_ref_point_tag(), tm); + +#ifdef CGAL_LINKED_WITH_TBB + if constexpr(std::is_same_v) + { + tbb::parallel_for(std::size_t(0), faces(tm).size(), [&](std::size_t i){ + face_descriptor f(i); + put(bb_map, f, face_bbox(f, tm)); + put(rp_map, f, to_input(get(vpm, target(halfedge(f, tm), tm))) ); + }); + } + else +#endif + { + for(face_descriptor f : faces(tm)){ + put(bb_map, f, face_bbox(f, tm)); + put(rp_map, f, to_input(get(vpm, target(halfedge(f, tm), tm))) ); + } + } + Compute_bbox compute_bbox(bb_map); + Split_primitives split_primitives(rp_map); + tree.template custom_build(compute_bbox, split_primitives); + } +}; + // For exact side_of_triangle_mesh template faces; - std::set interior_vertices; + std::vector interior_vertices; std::vector interior_edges; std::vector shared_edges; bool is_initialized; @@ -682,7 +783,7 @@ template void extract_patch_simplices( PolygonMesh& pm, std::vector::face_descriptor>& patch_faces, - std::set::vertex_descriptor>& interior_vertices, + std::vector::vertex_descriptor>& interior_vertices, std::vector::halfedge_descriptor>& interior_edges, std::vector::halfedge_descriptor>& shared_edges, const IsIntersectionEdge& is_intersection_edge) @@ -717,10 +818,10 @@ void extract_patch_simplices( for(halfedge_descriptor h : interior_edges) { - if ( !border_vertices.count( target(h,pm) ) ) - interior_vertices.insert( target(h,pm) ); - if ( !border_vertices.count( source(h,pm) ) ) - interior_vertices.insert( source(h,pm) ); + if( halfedge(target(h, pm),pm) == h && !border_vertices.count( target(h,pm) )) // We insert a vertex only once + interior_vertices.push_back( target(h,pm) ); + if( halfedge(source(h, pm),pm) == opposite(h, pm) && !border_vertices.count( source(h,pm) )) // We insert a vertex only once + interior_vertices.push_back( source(h,pm) ); } } @@ -875,10 +976,11 @@ void import_polyline( EdgeMap& pm1_to_output_edges, EdgeMap& pm2_to_output_edges, VertexMap& pm1_to_output_vertices, + VertexMap& pm2_to_output_vertices, const IntersectionEdgeMap& intersection_edges1, const IntersectionEdgeMap& intersection_edges2, const VertexPointMap1& vpm1, - const VertexPointMap2& /*vpm2*/, + const VertexPointMap2& vpm2, const VertexPointMapOut& vpm_out, std::vector ::edge_descriptor>& output_shared_edges, @@ -891,6 +993,15 @@ void import_polyline( output_shared_edges.push_back(add_edge(output)); halfedge_descriptor h_out = halfedge(output_shared_edges.back(),output); + auto set_output_vertex = [&](vertex_descriptor v, halfedge_descriptor h_out){ + user_visitor.before_vertex_copy(v, pm1, output); + vertex_descriptor new_v = add_vertex(output); + set_halfedge(new_v, opposite(h_out, output),output); + put(vpm_out, new_v, get(vpm1, v)); + user_visitor.after_vertex_copy(v, pm1, new_v, output); + return new_v; + }; + //make sure the first vertex does not already exist vertex_descriptor src = GT::null_vertex(); std::pair< typename VertexMap::iterator, bool > insert_res= @@ -898,12 +1009,9 @@ void import_polyline( if( insert_res.second ) { - user_visitor.before_vertex_copy(source(h1,pm1), pm1, output); - src = add_vertex(output); - set_halfedge(src, opposite(h_out,output),output); - put(vpm_out, src, get(vpm1, source(h1,pm1))); - user_visitor.after_vertex_copy(source(h1,pm1), pm1, src, output); + src = set_output_vertex(source(h1, pm1), h_out); insert_res.first->second = src; + pm2_to_output_vertices.insert( std::make_pair( source(h2,pm2), src ) ); } else src = insert_res.first->second; @@ -916,22 +1024,16 @@ void import_polyline( pm1_to_output_vertices.insert( std::make_pair( target(h1,pm1), tgt ) ); if( insert_res.second ) { - user_visitor.before_vertex_copy(target(h1,pm1), pm1, output); - tgt = add_vertex(output); - set_halfedge(tgt, h_out, output); - put(vpm_out, tgt, get(vpm1, target(h1,pm1))); - user_visitor.after_vertex_copy(target(h1,pm1), pm1, tgt, output); + tgt = set_output_vertex(target(h1, pm1), h_out); insert_res.first->second = tgt; + pm2_to_output_vertices.insert( std::make_pair( target(h2,pm2), tgt ) ); } else tgt = insert_res.first->second; } else{ - user_visitor.before_vertex_copy(target(h1,pm1), pm1, output); - tgt = add_vertex(output); - set_halfedge(tgt, h_out, output); - put(vpm_out, tgt, get(vpm1, target(h1,pm1))); - user_visitor.after_vertex_copy(target(h1,pm1), pm1, tgt, output); + tgt = set_output_vertex(target(h1, pm1), h_out); + pm2_to_output_vertices.insert( std::make_pair( target(h2,pm2), tgt ) ); } //update source and target vertex of the edge created @@ -943,10 +1045,8 @@ void import_polyline( halfedge_descriptor prev2=h2; //set the correspondence - pm1_to_output_edges.insert( - std::make_pair(edge(prev1, pm1), edge(prev_out, output)) ); - pm2_to_output_edges.insert( - std::make_pair(edge(prev2, pm2), edge(prev_out, output)) ); + pm1_to_output_edges.insert( std::make_pair(edge(prev1, pm1), edge(prev_out, output)) ); + pm2_to_output_edges.insert( std::make_pair(edge(prev2, pm2), edge(prev_out, output)) ); user_visitor.intersection_edge_copy(prev1, pm1, prev2, pm2, h_out, output); @@ -965,23 +1065,14 @@ void import_polyline( //if this is the final segment, only create a target vertex if it does not exist if (i+1!=nb_segments) { - user_visitor.before_vertex_copy(target(h1,pm1), pm1, output); - tgt=add_vertex(output); - set_halfedge(tgt, h_out, output); - put(vpm_out, tgt, get(vpm1, target(h1,pm1))); - user_visitor.after_vertex_copy(target(h1,pm1), pm1, tgt, output); + tgt = set_output_vertex(target(h1, pm1), h_out); } else{ std::pair< typename VertexMap::iterator, bool > insert_res = pm1_to_output_vertices.insert(std::make_pair(target(h1,pm1), tgt)); if (insert_res.second) { - user_visitor.before_vertex_copy(target(h1,pm1), pm1, output); - tgt=add_vertex(output); - set_halfedge(tgt, h_out, output); - put(vpm_out, tgt, get(vpm1, target(h1,pm1))); - insert_res.first->second = tgt; - user_visitor.after_vertex_copy(target(h1,pm1), pm1, tgt, output); + tgt = set_output_vertex(target(h1, pm1), h_out); } else tgt = insert_res.first->second; @@ -995,10 +1086,8 @@ void import_polyline( prev2 = h2; src = tgt; - pm1_to_output_edges.insert( - std::make_pair(edge(prev1, pm1), edge(prev_out, output)) ); - pm2_to_output_edges.insert( - std::make_pair(edge(prev2, pm2), edge(prev_out, output)) ); + pm1_to_output_edges.insert( std::make_pair(edge(prev1, pm1), edge(prev_out, output)) ); + pm2_to_output_edges.insert( std::make_pair(edge(prev2, pm2), edge(prev_out, output)) ); } } @@ -1091,6 +1180,7 @@ struct Triangle_mesh_extension_helper template < bool reverse_patch_orientation, class TriangleMesh, class PatchContainer, + // class VertexToPointMap, class VertexPointMap, class VertexPointMapOut, class EdgeMarkMapOut, @@ -1100,6 +1190,7 @@ void append_patches_to_triangle_mesh( TriangleMesh& output, const boost::dynamic_bitset<>& patches_to_append, PatchContainer& patches, + // const VertexToPointMap &tm_to_output_vertices, const VertexPointMapOut& vpm_out, const VertexPointMap& vpm_tm, EdgeMarkMapOut& edge_mark_map_out, @@ -1162,29 +1253,31 @@ void append_patches_to_triangle_mesh( CGAL_assertion(is_border(opposite(new_h,output), output)); //create a copy of interior vertices only once - if ( halfedge(target(h,tm),tm)==h && - patch.interior_vertices.count(target(h, tm)) ) - { - user_visitor.before_vertex_copy(target(h,tm), tm, output); - vertex_descriptor v = add_vertex(output); - set_halfedge(v, new_h, output); - set_target(new_h, v, output); - put(vpm_out, v, get(vpm_tm, target(h, tm) ) ); - user_visitor.after_vertex_copy(target(h,tm), tm, v, output); - interior_vertex_halfedges.push_back( new_h ); - } - if ( halfedge(source(h,tm),tm)==opposite(h,tm) && - patch.interior_vertices.count(source(h,tm)) ) - { - user_visitor.before_vertex_copy(source(h,tm), tm, output); - vertex_descriptor v = add_vertex(output); - halfedge_descriptor new_h_opp = opposite(new_h, output); - set_halfedge(v, new_h_opp, output); - set_target(new_h_opp, v, output); - put(vpm_out, v, get(vpm_tm, source(h, tm) ) ); - user_visitor.after_vertex_copy(source(h,tm), tm, v, output); - interior_vertex_halfedges.push_back( new_h_opp ); - } + // if ( halfedge(target(h,tm),tm)==h && + // tm_to_output_vertices[target(h, tm)] == std::numeric_limits::max()) + // { + // user_visitor.before_vertex_copy(target(h,tm), tm, output); + // vertex_descriptor v = add_vertex(output); + // set_halfedge(v, new_h, output); + // set_target(new_h, v, output); + // put(vpm_out, v, get(vpm_tm, target(h, tm) ) ); + // user_visitor.after_vertex_copy(target(h,tm), tm, v, output); + // interior_vertex_halfedges.push_back( new_h ); + // tm_to_output_vertices[target(h, tm)] = v; + // } + // if ( halfedge(source(h,tm),tm)==opposite(h,tm) && + // tm_to_output_vertices[source(h, tm)] == std::numeric_limits::max() ) + // { + // user_visitor.before_vertex_copy(source(h,tm), tm, output); + // vertex_descriptor v = add_vertex(output); + // halfedge_descriptor new_h_opp = opposite(new_h, output); + // set_halfedge(v, new_h_opp, output); + // set_target(new_h_opp, v, output); + // put(vpm_out, v, get(vpm_tm, source(h, tm) ) ); + // user_visitor.after_vertex_copy(source(h,tm), tm, v, output); + // interior_vertex_halfedges.push_back( new_h_opp ); + // tm_to_output_vertices[source(h, tm)] = v; + // } } } @@ -1236,12 +1329,12 @@ void append_patches_to_triangle_mesh( vertex_descriptor tgt = target(h, tm); if (reverse_patch_orientation) std::swap(src, tgt); - if ( !patch.interior_vertices.count(src) ) - border_halfedges_source_to_link.push_back(helper.get_hedge(h)); - if ( !patch.interior_vertices.count(tgt) ){ - border_halfedges_target_to_link.push_back(helper.get_hedge(h)); - continue; // since the next halfedge should not be in the same patch - } + // if ( !patch.interior_vertices.count(src) ) + // border_halfedges_source_to_link.push_back(helper.get_hedge(h)); + // if ( !patch.interior_vertices.count(tgt) ){ + // border_halfedges_target_to_link.push_back(helper.get_hedge(h)); + // continue; // since the next halfedge should not be in the same patch + // } CGAL_assertion( is_border(h, tm) && is_border(prev(h, tm),tm) && is_border(next(h, tm),tm)); @@ -1324,6 +1417,349 @@ void append_patches_to_triangle_mesh( } } +template +void import_vertices_from_polyline( + typename boost::graph_traits::halfedge_descriptor h1, + typename boost::graph_traits::halfedge_descriptor h2, + const PolygonMesh& pm1, + const PolygonMesh& pm2, + std::size_t nb_segments, + VertexMap& pm1_to_output_points, + VertexMap& pm2_to_output_points, + const IntersectionEdgeMap& intersection_edges1, + const IntersectionEdgeMap& intersection_edges2, + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2, + PointVector& pts, + UserVisitor& user_visitor) +{ + typedef boost::graph_traits GT; + typedef typename GT::halfedge_descriptor halfedge_descriptor; + typedef typename GT::vertex_descriptor vertex_descriptor; + + auto insert_new_point = [&](vertex_descriptor v1, vertex_descriptor v2){ + std::size_t idx = pts.size(); + pts.push_back(get(vpm1, v1)); + // pm1_to_output_points[v1] = idx; + // pm2_to_output_points[v2] = idx; + put( pm1_to_output_points, v1, idx); + put( pm2_to_output_points, v2, idx); + }; + + //make sure the first vertex does not already exist + vertex_descriptor v1 = source(h1, pm1); + vertex_descriptor v2 = source(h2, pm2); + // if( pm1_to_output_points.find( v1 ) == pm1_to_output_points.end()){ + // CGAL_assertion(pm2_to_output_points.find( v2 ) == pm2_to_output_points.end()); + if( get(pm1_to_output_points, v1) == std::numeric_limits::max()){ + insert_new_point(v1, v2); + } + + //make sure the target vertex does not already exist if it is a polyline endpoint + v1 = target(h1, pm1); + v2 = target(h2, pm2); + // if ( nb_segments!=1 || (pm1_to_output_points.find( v1 ) == pm1_to_output_points.end())){ + // CGAL_assertion(pm2_to_output_points.find( v2 ) == pm2_to_output_points.end()); + if(nb_segments!=1 || get(pm1_to_output_points, v1) == std::numeric_limits::max()){ + insert_new_point(v1, v2); + } + + for (std::size_t i=1; i::max()){ + // CGAL_assertion(pm2_to_output_points.find( v2 ) == pm2_to_output_points.end()); + insert_new_point(v1, v2); + } + } +} + +template < class PointRange, + class FaceRange, + class IntersectionEdgeMap, + class VertexPointMap1, + class VertexPointMap2, + class EdgeMarkMap1, + class EdgeMarkMap2, + class IntersectionPolylines, + class PatchContainer1, + class PatchContainer2, + class UserVisitor> +void fill_triangle_soup( + PointRange &pts, + FaceRange &faces, + const boost::dynamic_bitset<>& patches_of_tm1_to_import, + const boost::dynamic_bitset<>& patches_of_tm2_to_import, + PatchContainer1& patches_of_tm1, + PatchContainer2& patches_of_tm2, + bool reverse_orientation_of_patches_from_tm1, + bool reverse_orientation_of_patches_from_tm2, + const IntersectionPolylines& polylines, + const IntersectionEdgeMap& intersection_edges1, + const IntersectionEdgeMap& intersection_edges2, + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2, + const EdgeMarkMap1& edge_mark_map1, + const EdgeMarkMap2& edge_mark_map2, + UserVisitor& user_visitor) +{ + using TriangleMesh = std::remove_reference_t; + using Point = typename VertexPointMap1::value_type; + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using face_descriptor = typename boost::graph_traits::face_descriptor; + + using V2V_tag = typename CGAL::dynamic_vertex_property_t; + using Vertex_to_point_map = typename boost::property_map::const_type; + + using Concurrency_tag = Parallel_if_available_tag; + const std::size_t parallel_cutoff = 10000; + const bool parallel_execution = std::is_same_v; + +#ifdef CGAL_OUTPUT_BUILDER_RUNNING_TIME + Real_timer t; + t.start(); +#endif + + const TriangleMesh &tm1 = patches_of_tm1.pm; + const TriangleMesh &tm2 = patches_of_tm2.pm; + Vertex_to_point_map v2p_tm1 = get(V2V_tag(), tm1, std::numeric_limits::max()); + Vertex_to_point_map v2p_tm2 = get(V2V_tag(), tm2, std::numeric_limits::max()); + + std::size_t nb_polylines = polylines.lengths.size(); + for (std::size_t i=0; i < nb_polylines; ++i) + if (!polylines.to_skip.test(i)) + import_vertices_from_polyline(polylines.tm1[i], polylines.tm2[i], + tm1, tm2, + polylines.lengths[i], + v2p_tm1, v2p_tm2, + intersection_edges1, intersection_edges2, + vpm1, vpm2, pts, + user_visitor); + +#ifdef CGAL_OUTPUT_BUILDER_RUNNING_TIME + std::cout << "import vertices from polylines time: " << t.time() << std::endl; +#endif + + // Get ids of patch to append + std::vector ids_of_patches_to_append_from_tm1; + std::vector ids_of_patches_to_append_from_tm2; + ids_of_patches_to_append_from_tm1.reserve(patches_of_tm1_to_import.count()); + for (std::size_t i= patches_of_tm1_to_import.find_first(); + i < patches_of_tm1_to_import.npos; + i = patches_of_tm1_to_import.find_next(i)){ + ids_of_patches_to_append_from_tm1.push_back(i); + } + ids_of_patches_to_append_from_tm2.reserve(patches_of_tm2_to_import.count()); + for (std::size_t i= patches_of_tm2_to_import.find_first(); + i < patches_of_tm2_to_import.npos; + i = patches_of_tm2_to_import.find_next(i)){ + ids_of_patches_to_append_from_tm2.push_back(i); + } + + // Compute final sizes + std::size_t pn = pts.size(); + std::size_t fpn = pts.size(); + for (std::size_t i : ids_of_patches_to_append_from_tm1) + fpn += patches_of_tm1[i].interior_vertices.size(); + for (std::size_t i : ids_of_patches_to_append_from_tm2) + fpn += patches_of_tm2[i].interior_vertices.size(); + pts.resize(fpn); + + std::size_t fn = 0; + std::size_t ffn = 0; + for (std::size_t i : ids_of_patches_to_append_from_tm1) + ffn += patches_of_tm1[i].faces.size(); + for (std::size_t i : ids_of_patches_to_append_from_tm2) + ffn += patches_of_tm2[i].faces.size(); + faces.resize(ffn); + + // Append pts + for (std::size_t i : ids_of_patches_to_append_from_tm1){ + std::vector< vertex_descriptor > interior_vertices(patches_of_tm1[i].interior_vertices.begin(), patches_of_tm1[i].interior_vertices.end()); +#ifdef CGAL_LINKED_WITH_TBB + if constexpr(parallel_execution){ + if( interior_vertices.size() > parallel_cutoff) + tbb::parallel_for( + tbb::blocked_range(0, interior_vertices.size()), + [&](const tbb::blocked_range& r){ + for (std::size_t i = r.begin(); i != r.end(); ++i){ + pts[i+pn] = get(vpm1, interior_vertices[i]); + put(v2p_tm1, interior_vertices[i], i+pn); + } + } + ); + else + for(std::size_t i = 0; i interior_vertices(patches_of_tm2[i].interior_vertices.begin(), patches_of_tm2[i].interior_vertices.end()); +#ifdef CGAL_LINKED_WITH_TBB + if constexpr(parallel_execution){ + if( interior_vertices.size() > parallel_cutoff) + tbb::parallel_for( + tbb::blocked_range(0, interior_vertices.size()), + [&](const tbb::blocked_range& r){ + for (std::size_t i = r.begin(); i != r.end(); ++i){ + pts[i+pn] = get(vpm2, interior_vertices[i]); + put(v2p_tm2, interior_vertices[i], i+pn); + } + } + ); + else + for(std::size_t i = 0; i parallel_cutoff) + tbb::parallel_for( + tbb::blocked_range(0, tm1_faces.size()), + [&](const tbb::blocked_range& r){ + for (std::size_t i = r.begin(); i != r.end(); ++i){ + auto h = halfedge(tm1_faces[i], tm1); + + std::size_t id1 = get(v2p_tm1, source(h, tm1)); + std::size_t id2 = get(v2p_tm1, target(h, tm1)); + std::size_t id3 = get(v2p_tm1, target(next(h, tm1), tm1)); + + if (reverse_orientation_of_patches_from_tm1) + faces[fn + i] = {id1, id3, id2}; + else + faces[fn + i] = {id1, id2, id3}; + } + } + ); + else + for(std::size_t i = 0; i parallel_cutoff) + tbb::parallel_for( + tbb::blocked_range(0, tm2_faces.size()), + [&](const tbb::blocked_range& r){ + for (std::size_t i = r.begin(); i != r.end(); ++i){ + auto h = halfedge(tm2_faces[i], tm2); + + std::size_t id1 = get(v2p_tm2, source(h, tm2)); + std::size_t id2 = get(v2p_tm2, target(h, tm2)); + std::size_t id3 = get(v2p_tm2, target(next(h, tm2), tm2)); + + if (reverse_orientation_of_patches_from_tm2) + faces[fn + i] = {id1, id3, id2}; + else + faces[fn + i] = {id1, id2, id3}; + } + } + ); + else + for(std::size_t i = 0; i tm1_to_output_vertices; - std::unordered_map tm1_to_output_edges, - tm2_to_output_edges; + std::unordered_map tm1_to_output_vertices, tm2_to_output_vertices; + std::unordered_map tm1_to_output_edges, tm2_to_output_edges; for (std::size_t i=0; i < nb_polylines; ++i) if (!polylines.to_skip.test(i)) @@ -1380,7 +1814,7 @@ void fill_new_triangle_mesh( patches_of_tm2.pm, polylines.lengths[i], tm1_to_output_edges, tm2_to_output_edges, - tm1_to_output_vertices, + tm1_to_output_vertices, tm2_to_output_vertices, intersection_edges1, intersection_edges2, vpm1, vpm2, vpm_out, output_shared_edges, @@ -1391,6 +1825,7 @@ void fill_new_triangle_mesh( append_patches_to_triangle_mesh(output, patches_of_tm1_to_import, patches_of_tm1, + // tm1_to_output_vertices, vpm_out, vpm1, edge_mark_map_out, @@ -1401,6 +1836,7 @@ void fill_new_triangle_mesh( append_patches_to_triangle_mesh(output, patches_of_tm1_to_import, patches_of_tm1, + // tm1_to_output_vertices, vpm_out, vpm1, edge_mark_map_out, @@ -1413,6 +1849,7 @@ void fill_new_triangle_mesh( append_patches_to_triangle_mesh(output, patches_of_tm2_to_import, patches_of_tm2, + // tm2_to_output_vertices, vpm_out, vpm2, edge_mark_map_out, @@ -1423,6 +1860,7 @@ void fill_new_triangle_mesh( append_patches_to_triangle_mesh(output, patches_of_tm2_to_import, patches_of_tm2, + // tm2_to_output_vertices, vpm_out, vpm2, edge_mark_map_out, @@ -1721,19 +2159,19 @@ remove_patches(TriangleMesh& tm, if(is_border_edge(h, tm)) { if (is_border(h, tm)) h=opposite(h, tm); - if ( !patch.interior_vertices.count(target(h, tm)) ) - { - // look for the halfedge belonging to shared_edges - // having the prev pointer not correctly set - halfedge_descriptor nxt=next(h, tm); - while(!is_border(nxt, tm)) - nxt=next(opposite(nxt, tm), tm); - CGAL_assertion( is_border(nxt, tm) );//we marked it above! - // now update the prev pointer - halfedge_descriptor prv=prev(opposite(h, tm), tm); - set_next(prv, nxt, tm); - set_halfedge(target(prv, tm), prv, tm); - } + // if ( !patch.interior_vertices.count(target(h, tm)) ) + // { + // // look for the halfedge belonging to shared_edges + // // having the prev pointer not correctly set + // halfedge_descriptor nxt=next(h, tm); + // while(!is_border(nxt, tm)) + // nxt=next(opposite(nxt, tm), tm); + // CGAL_assertion( is_border(nxt, tm) );//we marked it above! + // // now update the prev pointer + // halfedge_descriptor prv=prev(opposite(h, tm), tm); + // set_next(prv, nxt, tm); + // set_halfedge(target(prv, tm), prv, tm); + // } } //now remove the simplices diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h index 6e0a2f38f9ce..ee64c0e54cef 100644 --- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h @@ -29,6 +29,10 @@ namespace CGAL { namespace Polygon_mesh_processing { namespace Corefinement { +#ifdef CGAL_LINKED_WITH_TBB +CGAL_MUTEX corefinement_intersection_callback_mutex; +#endif + template class Collect_face_bbox_per_edge_bbox { protected: @@ -57,7 +61,9 @@ class Collect_face_bbox_per_edge_bbox { { halfedge_descriptor fh = face_box.info(); halfedge_descriptor eh = edge_box.info(); - +#ifdef CGAL_LINKED_WITH_TBB + CGAL_SCOPED_LOCK(corefinement_intersection_callback_mutex); +#endif edge_to_faces[edge(eh,tm_edges)].insert(face(fh, tm_faces)); } @@ -86,7 +92,7 @@ class Collect_face_bbox_per_edge_bbox_with_coplanar_handling { 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::value_type Point; typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy; typedef CGAL::Box_intersection_d::Box_with_info_d Box; @@ -109,19 +115,18 @@ class Collect_face_bbox_per_edge_bbox_with_coplanar_handling { , visitor(visitor) {} - void operator()( const Box& face_box, const Box& edge_box) const { - halfedge_descriptor fh = face_box.info(); - halfedge_descriptor eh = edge_box.info(); + void operator()( halfedge_descriptor fh, halfedge_descriptor eh ) const { if(is_border(eh,tm_edges)) eh = opposite(eh, tm_edges); //check if the segment intersects the plane of the facet or if it is included in the plane - Point a = get(vpmap_tmf, source(fh, tm_faces)); - Point b = get(vpmap_tmf, target(fh, tm_faces)); - Point c = get(vpmap_tmf, target(next(fh, tm_faces), tm_faces)); + Point& a = get(vpmap_tmf, source(fh, tm_faces)); + Point& b = get(vpmap_tmf, target(fh, tm_faces)); + Point& c = get(vpmap_tmf, target(next(fh, tm_faces), tm_faces)); /// SHOULD_USE_TRAITS_TAG - const Orientation abcp = orientation(a,b,c, get(vpmap_tme, target(eh, tm_edges))); - const Orientation abcq = orientation(a,b,c, get(vpmap_tme, source(eh, tm_edges))); + using K = typename Kernel_traits::Kernel; + auto orientation = K().orientation_3_object(); + const auto [abcp,abcq] = orientation(a,b,c, get(vpmap_tme, target(eh, tm_edges)), get(vpmap_tme, source(eh, tm_edges))); if (abcp==abcq){ if (abcp!=COPLANAR){ return; //no intersection @@ -149,39 +154,45 @@ class Collect_face_bbox_per_edge_bbox_with_coplanar_handling { return; } // non-coplanar case +#ifdef CGAL_LINKED_WITH_TBB + CGAL_SCOPED_LOCK(corefinement_intersection_callback_mutex); +#endif edge_to_faces[edge(eh,tm_edges)].insert(face(fh, tm_faces)); } bool is_face_degenerated(halfedge_descriptor fh) const { - Point a = get(vpmap_tmf, source(fh, tm_faces)); - Point b = get(vpmap_tmf, target(fh, tm_faces)); - Point c = get(vpmap_tmf, target(next(fh, tm_faces), tm_faces)); + Point& a = get(vpmap_tmf, source(fh, tm_faces)); + Point& b = get(vpmap_tmf, target(fh, tm_faces)); + Point& c = get(vpmap_tmf, target(next(fh, tm_faces), tm_faces)); return collinear(a, b, c); } bool are_edge_faces_degenerated(halfedge_descriptor eh) const { - Point a = get(vpmap_tme, source(eh, tm_edges)); - Point b = get(vpmap_tme, target(eh, tm_edges)); + Point& a = get(vpmap_tme, source(eh, tm_edges)); + Point& b = get(vpmap_tme, target(eh, tm_edges)); if(!is_border(eh,tm_edges)) { - Point c = get(vpmap_tme, target(next(eh, tm_edges), tm_edges)); + Point& c = get(vpmap_tme, target(next(eh, tm_edges), tm_edges)); if (collinear(a, b, c)) return true; } eh = opposite(eh, tm_edges); if(!is_border(eh,tm_edges)) { - Point c = get(vpmap_tme, target(next(eh, tm_edges), tm_edges)); + Point& c = get(vpmap_tme, target(next(eh, tm_edges), tm_edges)); if (collinear(a, b, c)) return true; } return false; } + void operator()( const Box& face_box, const Box& edge_box) const { + operator()(face_box.info(), edge_box.info()); + } void operator()(const Box* face_box_ptr, const Box* edge_box_ptr) const { @@ -353,6 +364,9 @@ class Collect_face_bbox_per_edge_bbox_with_coplanar_handling_one_mesh { } // non-coplanar case +#ifdef CGAL_LINKED_WITH_TBB + CGAL_SCOPED_LOCK(corefinement_intersection_callback_mutex); +#endif edge_to_faces[edge(eh,tm)].insert(face(fh, tm)); } diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h index 2ad06651e85b..74bee6c3a1c5 100644 --- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h @@ -15,9 +15,7 @@ #include - #include -#include #include #include #include @@ -25,6 +23,7 @@ #include #include #include +#include #include #include @@ -51,8 +50,7 @@ struct Triple_intersection_exception : // The algorithm works as follow: // From each triangle mesh, we extract a set of segments from the edges and a // set of triangles from the faces. -// We use Box_intersection_d to filter intersection between the segments from -// one mesh with the triangles from the other one. +// We use AABB_trees to filter pair of candidates intersecting faces and extract the intersection of between an edge and a face. // From this filtered set, for each pair (segment,triangle), we look at the // intersection type. If not empty, we can have three different cases // 1)the segment intersect the interior of the triangle: @@ -195,7 +193,8 @@ struct Node_id_set { template< class TriangleMesh, class VertexPointMap1, class VertexPointMap2, - class Node_visitor=Default_surface_intersection_visitor + class Node_visitor=Default_surface_intersection_visitor, + class Concurrency_tag=Sequential_tag > class Intersection_of_triangle_meshes { @@ -205,9 +204,6 @@ class Intersection_of_triangle_meshes typedef typename graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename graph_traits::vertex_descriptor vertex_descriptor; - typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy; - typedef CGAL::Box_intersection_d::Box_with_info_d Box; - typedef std::unordered_set Face_set; typedef std::unordered_map Edge_to_faces; @@ -243,218 +239,127 @@ class Intersection_of_triangle_meshes CGAL_assertion_code(bool doing_autorefinement;) // member functions - template - void filter_intersections(const TriangleMesh& tm_f, - const TriangleMesh& tm_e, - const VPMF& vpm_f, - const VPME& vpm_e, - const Non_manifold_feature_map& non_manifold_feature_map, + + // template + template + void filter_intersections(const TriangleMesh& tm1, + const TriangleMesh& tm2, + const VPMA& vpm1, + const VPMB& vpm2, + const Non_manifold_feature_map& non_manifold_feature_map_1, + const Non_manifold_feature_map& non_manifold_feature_map_2, bool throw_on_self_intersection, - std::set& tm_f_faces, - std::set& tm_e_faces, - Bbox_3 g_bb, - IFSM is_shared_map_f, - IESM is_shared_map_e, - bool run_check) + std::set& tm1_faces, + std::set& tm2_faces, + Bbox_3 tm_1_bb, + Bbox_3 tm_2_bb) { - std::vector face_boxes, edge_boxes; - std::vector face_boxes_ptr, edge_boxes_ptr; - - face_boxes.reserve(num_faces(tm_f)); - face_boxes_ptr.reserve(num_faces(tm_f)); - for(face_descriptor fd : faces(tm_f)) - { - if (get(is_shared_map_f,fd)) continue; - halfedge_descriptor h=halfedge(fd,tm_f); - Bbox_3 bb = get(vpm_f,source(h,tm_f)).bbox() + - get(vpm_f,target(h,tm_f)).bbox() + - get(vpm_f,target(next(h,tm_f),tm_f)).bbox(); - if (do_overlap(bb, g_bb)) - { - face_boxes.emplace_back(bb, h); - face_boxes_ptr.push_back( &face_boxes.back() ); - } - } - - edge_boxes.reserve(num_edges(tm_e)); - edge_boxes_ptr.reserve(num_edges(tm_e)); - if (non_manifold_feature_map.non_manifold_edges.empty()) - // general manifold case - for(edge_descriptor ed : edges(tm_e)) - { - if (get(is_shared_map_e,ed)) continue; - halfedge_descriptor h=halfedge(ed,tm_e); - Bbox_3 bb = get(vpm_e,source(h,tm_e)).bbox() + - get(vpm_e,target(h,tm_e)).bbox(); - - if (do_overlap(bb, g_bb)) - { - edge_boxes.emplace_back(bb,h); - edge_boxes_ptr.push_back( &edge_boxes.back() ); - } - } - else - // non-manifold case - for(edge_descriptor ed : edges(tm_e)) - { - std::size_t eid=get(non_manifold_feature_map.e_nm_id, ed); - halfedge_descriptor h=halfedge(ed,tm_e); - // insert only one copy of a non-manifold edge - if (eid!=NM_NID) - { - if (non_manifold_feature_map.non_manifold_edges[eid].front()!=ed) - continue; - else - // make sure the halfedge used is consistent with stored one - h = halfedge(non_manifold_feature_map.non_manifold_edges[eid].front(), tm_e); - } - if (get(is_shared_map_e,edge(h,tm_e))) continue; - Bbox_3 bb = get(vpm_e,source(h,tm_e)).bbox() + - get(vpm_e,target(h,tm_e)).bbox(); - - if (do_overlap(bb, g_bb)) - { - edge_boxes.emplace_back(bb,h); - edge_boxes_ptr.push_back( &edge_boxes.back() ); - } - } +#ifdef CGAL_OUTPUT_BUILDER_RUNNING_TIME + CGAL::Real_timer t; + t.start(); +#endif - /// \todo experiments different cutoff values - std::ptrdiff_t cutoff = 2 * std::ptrdiff_t( - std::sqrt(face_boxes.size()+edge_boxes.size()) ); + using GT = typename GetGeomTraits::type; + using AABB_tree_helper = AABB_tree_build_helper; + using Tree = typename AABB_tree_helper::Tree; + AABB_tree_helper helper; - Edge_to_faces& edge_to_faces = &tm_e < &tm_f - ? stm_edge_to_ltm_faces - : ltm_edge_to_stm_faces; + Tree tree1(faces(tm1).begin(), faces(tm1).end(), tm1); + Tree tree2(faces(tm2).begin(), faces(tm2).end(), tm2); + // std::cout << "Computing Bbox and reference point " << t.time() << std::endl; #ifdef DO_NOT_HANDLE_COPLANAR_FACES typedef Collect_face_bbox_per_edge_bbox Callback; - Callback callback(tm_f, tm_e, edge_to_faces); + Callback callback12(tm1, tm2, ltm_edge_to_stm_faces); + Callback callback21(tm2, tm1, stm_edge_to_ltm_faces); #else typedef Collect_face_bbox_per_edge_bbox_with_coplanar_handling< - TriangleMesh, VPMF, VPME, Edge_to_faces, Coplanar_face_set, Node_visitor> + TriangleMesh, VPMA, VPMB, Edge_to_faces, Coplanar_face_set, Node_visitor> Callback; - Callback callback(tm_f, tm_e, vpm_f, vpm_e, edge_to_faces, coplanar_faces, visitor); + Callback callback12(tm1, tm2, vpm1, vpm2, ltm_edge_to_stm_faces, coplanar_faces, visitor); + Callback callback21(tm2, tm1, vpm2, vpm1, stm_edge_to_ltm_faces, coplanar_faces, visitor); #endif - //using pointers in box_intersection_d is about 10% faster - if (throw_on_self_intersection){ - Callback_with_self_intersection_report callback_si(callback, tm_f_faces, tm_e_faces); - CGAL::box_intersection_d(face_boxes_ptr.begin(), face_boxes_ptr.end(), - edge_boxes_ptr.begin(), edge_boxes_ptr.end(), - callback_si, cutoff); - if (run_check && callback_si.self_intersections_found()) - throw Self_intersection_exception(); - } - else { - if (const_mesh_ptr==&tm_e) - { - // tm_f might feature degenerate faces - auto filtered_callback = [&callback](const Box* fb, const Box* eb) - { - if (!callback.is_face_degenerated(fb->info())) - callback(fb, eb); - }; - CGAL::box_intersection_d( face_boxes_ptr.begin(), face_boxes_ptr.end(), - edge_boxes_ptr.begin(), edge_boxes_ptr.end(), - filtered_callback, cutoff ); - } - else - { - if (const_mesh_ptr==&tm_f) - { - // tm_e might feature degenerate edges - auto filtered_callback = [&,this](const Box* fb, const Box* eb) - { - if (get(vpm_e, source(eb->info(), tm_e)) != get(vpm_e, target(eb->info(), tm_e))) - callback(fb, eb); - else - { - halfedge_descriptor hf = fb->info(); - halfedge_descriptor he = eb->info(); - for (int i=0; i<2; ++i) - { - if (!is_border(he, tm_e)) - { - if ( get(vpm_e, target(next(he, tm_e), tm_e))==get(vpm_e, target(he, tm_e)) && - coplanar(get(vpm_f, source(hf, tm_f)), - get(vpm_f, target(hf, tm_f)), - get(vpm_f, target(next(hf, tm_f), tm_f)), - get(vpm_e, target(he, tm_e))) ) - { - coplanar_faces.insert( - &tm_e < &tm_f - ? std::make_pair(face(he, tm_e), face(hf, tm_f)) - : std::make_pair(face(hf, tm_f), face(he, tm_e)) - ); - } - } - he=opposite(he, tm_e); - } - } - }; - CGAL::box_intersection_d( face_boxes_ptr.begin(), face_boxes_ptr.end(), - edge_boxes_ptr.begin(), edge_boxes_ptr.end(), - filtered_callback, cutoff ); - } - else - CGAL::box_intersection_d( face_boxes_ptr.begin(), face_boxes_ptr.end(), - edge_boxes_ptr.begin(), edge_boxes_ptr.end(), - callback, cutoff ); - } - } - } - // for autorefinement - template - void filter_intersections(const TriangleMesh& tm, - const VPM& vpm) - { - std::vector face_boxes, edge_boxes; - std::vector face_boxes_ptr, edge_boxes_ptr; - - face_boxes.reserve(num_faces(tm)); - face_boxes_ptr.reserve(num_faces(tm)); - for(face_descriptor fd : faces(tm)) + #ifdef CGAL_LINKED_WITH_TBB + if constexpr(std::is_same_v) { - halfedge_descriptor h=halfedge(fd,tm); - face_boxes.push_back( Box( - get(vpm,source(h,tm)).bbox() + - get(vpm,target(h,tm)).bbox() + - get(vpm,target(next(h,tm),tm)).bbox(), - h ) ); - face_boxes_ptr.push_back( &face_boxes.back() ); + oneapi::tbb::task_group tg; + tg.run([&]{ helper.template build(tree1, tm1, vpm1); }); + helper.template build(tree2, tm2, vpm2); + tg.wait(); + + // std::cout << "Building AABB tree " << t.time() << std::endl; + + tbb::concurrent_vector> inter; + CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter), parameters::concurrency_tag(Concurrency_tag())); + + // std::cout << "Compute " << inter.size() << " candidates " << t.time() << std::endl; + // for(const auto& [f_1, f_2]: inter){ + tbb::parallel_for(std::size_t(0), inter.size(), [&](std::size_t i){ + const auto& [f_1, f_2] = inter[i]; + + halfedge_descriptor hf1_0 = halfedge(f_1, tm1); + halfedge_descriptor hf1_1 = next(hf1_0, tm1); + halfedge_descriptor hf1_2 = next(hf1_1, tm1); + + halfedge_descriptor hf2_0 = halfedge(f_2, tm2); + halfedge_descriptor hf2_1 = next(hf2_0, tm2); + halfedge_descriptor hf2_2 = next(hf2_1, tm2); + + if (is_border(hf2_0, tm2) || hf2_0 < opposite(hf2_0, tm2)) + callback12(hf1_0, hf2_0); + if (is_border(hf2_1, tm2) || hf2_1 < opposite(hf2_1, tm2)) + callback12(hf1_0, hf2_1); + if (is_border(hf2_2, tm2) || hf2_2 < opposite(hf2_2, tm2)) + callback12(hf1_0, hf2_2); + + if (is_border(hf1_0, tm1) || hf1_0 < opposite(hf1_0, tm1)) + callback21(hf2_0, hf1_0); + if (is_border(hf1_1, tm1) || hf1_1 < opposite(hf1_1, tm1)) + callback21(hf2_0, hf1_1); + if (is_border(hf1_2, tm1) || hf1_2 < opposite(hf1_2, tm1)) + callback21(hf2_0, hf1_2); + }); + // } + + // std::cout << "process candidates " << t.time() << std::endl; } - - edge_boxes.reserve(num_edges(tm)); - edge_boxes_ptr.reserve(num_edges(tm)); - for(edge_descriptor ed : edges(tm)) + else + #endif { - halfedge_descriptor h=halfedge(ed,tm); - edge_boxes.push_back( Box( - get(vpm,source(h,tm)).bbox() + - get(vpm,target(h,tm)).bbox(), - h ) ); - edge_boxes_ptr.push_back( &edge_boxes.back() ); + helper.template build(tree1, tm1, vpm1); + helper.template build(tree2, tm2, vpm2); + // using InternOutputIterator= std::back_insert_iterator>>; + std::vector> inter; + CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter)); + for(const auto& [f_1, f_2]: inter){ + halfedge_descriptor hf1_0 = halfedge(f_1, tm1); + halfedge_descriptor hf1_1 = next(hf1_0, tm1); + halfedge_descriptor hf1_2 = next(hf1_1, tm1); + + halfedge_descriptor hf2_0 = halfedge(f_2, tm2); + halfedge_descriptor hf2_1 = next(hf2_0, tm2); + halfedge_descriptor hf2_2 = next(hf2_1, tm2); + + if (is_border(hf2_0, tm2) || hf2_0 < opposite(hf2_0, tm2)) + callback12(hf1_0, hf2_0); + if (is_border(hf2_1, tm2) || hf2_1 < opposite(hf2_1, tm2)) + callback12(hf1_0, hf2_1); + if (is_border(hf2_2, tm2) || hf2_2 < opposite(hf2_2, tm2)) + callback12(hf1_0, hf2_2); + + if (is_border(hf1_0, tm1) || hf1_0 < opposite(hf1_0, tm1)) + callback21(hf2_0, hf1_0); + if (is_border(hf1_1, tm1) || hf1_1 < opposite(hf1_1, tm1)) + callback21(hf2_0, hf1_1); + if (is_border(hf1_2, tm1) || hf1_2 < opposite(hf1_2, tm1)) + callback21(hf2_0, hf1_2); + } } - - /// \todo experiments different cutoff values - std::ptrdiff_t cutoff = 2 * std::ptrdiff_t( - std::sqrt(face_boxes.size()+edge_boxes.size()) ); - - Edge_to_faces& edge_to_faces = stm_edge_to_ltm_faces; - - typedef Collect_face_bbox_per_edge_bbox_with_coplanar_handling_one_mesh< - TriangleMesh, VPM, Edge_to_faces, Coplanar_face_set> - Callback; - Callback callback(tm, vpm, edge_to_faces, coplanar_faces); - - //using pointers in box_intersection_d is about 10% faster - CGAL::box_intersection_d( face_boxes_ptr.begin(), face_boxes_ptr.end(), - edge_boxes_ptr.begin(), edge_boxes_ptr.end(), - callback, cutoff ); } + template std::pair get_or_create_node(const Cpl_inter_pt& ipt, @@ -943,8 +848,12 @@ class Intersection_of_triangle_meshes } } + CGAL_MUTEX m; + template - void compute_intersection_points(Edge_to_faces& tm1_edge_to_tm2_faces, + void compute_intersection_points(edge_descriptor e_1, + Face_set& fset, + Edge_to_faces& tm1_edge_to_tm2_faces, const TriangleMesh& tm1, const TriangleMesh& tm2, const VPM1& vpm1, @@ -955,35 +864,34 @@ class Intersection_of_triangle_meshes { typedef std::tuple Inter_type; - visitor.start_handling_edge_face_intersections(tm1_edge_to_tm2_faces.size()); + visitor.edge_face_intersections_step(); - for(typename Edge_to_faces::iterator it=tm1_edge_to_tm2_faces.begin(); - it!=tm1_edge_to_tm2_faces.end();++it) - { + CGAL_MUTEX m; + auto loop_step = [&](auto &pair){ visitor.edge_face_intersections_step(); - edge_descriptor e_1=it->first; + edge_descriptor e_1=pair.first; halfedge_descriptor h_1=halfedge(e_1,tm1); - Face_set& fset=it->second; + Face_set& fset=pair.second; while (!fset.empty()){ face_descriptor f_2=*fset.begin(); 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_2 + //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 { if (!nm_features_map_1.non_manifold_edges.empty()) { - std::size_t vid1 = get(nm_features_map_1.v_nm_id, target(h_1, tm1)); + std::size_t vid1 = get(nm_features_map_1.v_nm_id, source(h_1, tm1)); if (vid1 != NM_NID) { for (vertex_descriptor vd : nm_features_map_1.non_manifold_vertices[vid1]) { - std::copy(halfedges_around_target(vd,tm1).first, - halfedges_around_target(vd,tm1).second, + std::copy(halfedges_around_source(vd,tm1).first, + halfedges_around_source(vd,tm1).second, std::back_inserter(all_edges)); } if (all_edges.front()!=h_1) @@ -996,188 +904,167 @@ class Intersection_of_triangle_meshes } } else - std::copy(halfedges_around_target(h_1,tm1).first, - halfedges_around_target(h_1,tm1).second, + std::copy(halfedges_around_source(h_1,tm1).first, + halfedges_around_source(h_1,tm1).second, std::back_inserter(all_edges)); } else - std::copy(halfedges_around_target(h_1,tm1).first, - halfedges_around_target(h_1,tm1).second, + std::copy(halfedges_around_source(h_1,tm1).first, + halfedges_around_source(h_1,tm1).second, std::back_inserter(all_edges)); } - else{ - if ( std::get<2>(res) ) // is edge source in triangle plane + else + { + all_edges.push_back(h_1); + edge_descriptor e_1 = edge(h_1, tm1); + if (!nm_features_map_1.non_manifold_edges.empty()) { - if (!nm_features_map_1.non_manifold_edges.empty()) + std::size_t eid1 = get(nm_features_map_1.e_nm_id, e_1); + if (eid1 != NM_NID) { - std::size_t vid1 = get(nm_features_map_1.v_nm_id, source(h_1, tm1)); - if (vid1 != NM_NID) + CGAL_assertion( nm_features_map_1.non_manifold_edges[eid1][0]==e_1 ); + for (std::size_t k=1; + k::iterator pos = - std::find(all_edges.begin(), all_edges.end(), h_1); - CGAL_assertion(pos!=all_edges.end()); - std::swap(*pos, all_edges.front()); - } + edge_descriptor e_1b = nm_features_map_1.non_manifold_edges[eid1][k]; + // note that the orientation of the halfedge pushed back is + // not relevant for how it is used in the following + all_edges.push_back(halfedge(e_1b, tm1)); } - else - std::copy(halfedges_around_source(h_1,tm1).first, - halfedges_around_source(h_1,tm1).second, - std::back_inserter(all_edges)); } - else - 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_1); - edge_descriptor e_1 = edge(h_1, tm1); - if (!nm_features_map_1.non_manifold_edges.empty()) + } + } + 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()), + // all_edges.end(),std::get<1>(res),type); + // #endif + + CGAL_SCOPED_LOCK(m); + typename std::vector::iterator it_edge=all_edges.begin(); + switch(type){ + case COPLANAR_TRIANGLES: + #ifndef DO_NOT_HANDLE_COPLANAR_FACES + CGAL_error_msg("COPLANAR_TRIANGLES : this point should never be reached!"); + #else + //nothing needs to be done, cf. comments at the beginning of the file + #endif + break; + case EMPTY: + fset.erase(fset.begin()); + break; + + // Case when the edge pierces the face in its interior. + case ON_FACE: + { + CGAL_assertion(f_2==face(std::get<1>(res),tm2)); + + 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)); + 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); + // erase face from the list to test intersection with it_edge + if ( it_edge==all_edges.begin() ) { - std::size_t eid1 = get(nm_features_map_1.e_nm_id, e_1); - if (eid1 != NM_NID) - { - CGAL_assertion( nm_features_map_1.non_manifold_edges[eid1][0]==e_1 ); - for (std::size_t k=1; - ksecond.erase(f_2); } } - } - 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()), - // all_edges.end(),std::get<1>(res),type); - // #endif - - typename std::vector::iterator it_edge=all_edges.begin(); - switch(type){ - case COPLANAR_TRIANGLES: - #ifndef DO_NOT_HANDLE_COPLANAR_FACES - CGAL_error_msg("COPLANAR_TRIANGLES : this point should never be reached!"); - #else - //nothing needs to be done, cf. comments at the beginning of the file - #endif - break; - case EMPTY: - fset.erase(fset.begin()); - break; + } // end case ON_FACE + break; - // Case when the edge pierces the face in its interior. - case ON_FACE: - { - CGAL_assertion(f_2==face(std::get<1>(res),tm2)); - - 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)); - 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); - //erase face from the list to test intersection with it_edge - if ( it_edge==all_edges.begin() ) - { - fset.erase(fset.begin()); - } + // Case when the edge intersect one edge of the face. + 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); + + std::size_t eid2 = nm_features_map_2.non_manifold_edges.empty() + ? NM_NID + : get(nm_features_map_2.e_nm_id, edge(h_2, tm2)); + + if (eid2!=NM_NID) + h_2 = halfedge(nm_features_map_2.non_manifold_edges[eid2].front(), tm2); + + 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() ){ + CGAL_SCOPED_LOCK(m); + 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; + if( eid2 == NM_NID ) + cip_handle_case_edge(node_id,fset_bis,*it_edge,h_2,tm1,tm2); 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); + for (edge_descriptor e2 : nm_features_map_2.non_manifold_edges[eid2]) + cip_handle_case_edge(node_id,fset_bis,*it_edge,halfedge(e2, tm2),tm1,tm2); } } - } // end case ON_FACE - break; - - // Case when the edge intersect one edge of the face. - 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); - - std::size_t eid2 = nm_features_map_2.non_manifold_edges.empty() - ? NM_NID - : get(nm_features_map_2.e_nm_id, edge(h_2, tm2)); - - if (eid2!=NM_NID) - h_2 = halfedge(nm_features_map_2.non_manifold_edges[eid2].front(), tm2); - - 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=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; - if( eid2 == NM_NID ) - cip_handle_case_edge(node_id,fset_bis,*it_edge,h_2,tm1,tm2); - else - { - for (edge_descriptor e2 : nm_features_map_2.non_manifold_edges[eid2]) - cip_handle_case_edge(node_id,fset_bis,*it_edge,halfedge(e2, tm2),tm1,tm2); - } - } + else + { + if( eid2 == NM_NID ) + cip_handle_case_edge(node_id,&fset,*it_edge,h_2,tm1,tm2); else - { - if( eid2 == NM_NID ) - cip_handle_case_edge(node_id,&fset,*it_edge,h_2,tm1,tm2); - else - for (edge_descriptor e2 : nm_features_map_2.non_manifold_edges[eid2]) - cip_handle_case_edge(node_id,&fset,*it_edge,halfedge(e2, tm2),tm1,tm2); - } + for (edge_descriptor e2 : nm_features_map_2.non_manifold_edges[eid2]) + cip_handle_case_edge(node_id,&fset,*it_edge,halfedge(e2, tm2),tm1,tm2); } - } // end case ON_EDGE - break; + } + } // end case ON_EDGE + break; - 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 - //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)); - - std::size_t vid2 = nm_features_map_2.non_manifold_vertices.empty() - ? NM_NID - : get(nm_features_map_2.v_nm_id, target(h_2, tm2)); - - 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; - if( vid2 == NM_NID ) - cip_handle_case_vertex(node_id,fset_bis,*it_edge,h_2,tm1,tm2); - else - for (vertex_descriptor vd2 : nm_features_map_2.non_manifold_vertices[vid2]) - cip_handle_case_vertex(node_id,fset_bis,*it_edge,halfedge(vd2, tm2),tm1,tm2); - } + 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 + //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)); + + std::size_t vid2 = nm_features_map_2.non_manifold_vertices.empty() + ? NM_NID + : get(nm_features_map_2.v_nm_id, target(h_2, tm2)); + + for (;it_edge!=all_edges.end();++it_edge){ + if ( it_edge!=all_edges.begin() ){ + CGAL_SCOPED_LOCK(m); + 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; + if( vid2 == NM_NID ) + cip_handle_case_vertex(node_id,fset_bis,*it_edge,h_2,tm1,tm2); else - if( vid2 == NM_NID ) - cip_handle_case_vertex(node_id,&fset,*it_edge,h_2,tm1,tm2); - else - for (vertex_descriptor vd2 : nm_features_map_2.non_manifold_vertices[vid2]) - cip_handle_case_vertex(node_id,&fset,*it_edge,halfedge(vd2, tm2),tm1,tm2); + for (vertex_descriptor vd2 : nm_features_map_2.non_manifold_vertices[vid2]) + cip_handle_case_vertex(node_id,fset_bis,*it_edge,halfedge(vd2, tm2),tm1,tm2); } } // end case ON_VERTEX break; } // end switch on the type of the intersection } // end loop on all faces that intersect the edge - } // end loop on all entries (edges) in 'edge_to_face' + }; // end step on all entries (edges) in 'edge_to_face' + +#ifdef CGAL_LINKED_WITH_TBB + // if constexpr(std::is_same_v){ + // oneapi::tbb::parallel_for_each(tm1_edge_to_tm2_faces.begin(), tm1_edge_to_tm2_faces.end(), loop_step); + // } + // else +#endif + { + for(typename Edge_to_faces::iterator it=tm1_edge_to_tm2_faces.begin(); + it!=tm1_edge_to_tm2_faces.end();++it){ + loop_step(*it); + } + } + + CGAL_assertion(nodes.size()==unsigned(current_node+1)); visitor.end_handling_edge_face_intersections(); } @@ -1712,7 +1599,7 @@ class Intersection_of_triangle_meshes CGAL_assertion_code( doing_autorefinement=true; ) } -// setting maps of non manifold features + // setting maps of non manifold features void set_non_manifold_feature_map_1(internal_np::Param_not_found){} void set_non_manifold_feature_map_2(internal_np::Param_not_found){} void set_non_manifold_feature_map_1(const Non_manifold_feature_map& m) @@ -1757,192 +1644,13 @@ class Intersection_of_triangle_meshes std::set tm1_faces; std::set tm2_faces; - // preprocessing to detect identical faces/edges and - // ignore them in the intersection tests - std::vector tm1_verts, tm2_verts; - tm1_verts.reserve(num_vertices(tm1)); - tm2_verts.reserve(num_vertices(tm2)); - for (vertex_descriptor v : vertices(tm1)) - { - const auto& p = get(vpm1, v); - if (p.x()<=bb12.xmax() && p.x()>=bb12.xmin() && - p.y()<=bb12.ymax() && p.y()>=bb12.ymin() && - p.z()<=bb12.zmax() && p.z()>=bb12.zmin()) - { - tm1_verts.push_back(v); - } - } - for (vertex_descriptor v : vertices(tm2)) - { - const auto& p = get(vpm2, v); - if (p.x()<=bb12.xmax() && p.x()>=bb12.xmin() && - p.y()<=bb12.ymax() && p.y()>=bb12.ymin() && - p.z()<=bb12.zmax() && p.z()>=bb12.zmin()) - { - tm2_verts.push_back(v); - } - } - - std::sort(tm1_verts.begin(), tm1_verts.end(), - [&vpm1](vertex_descriptor v1, vertex_descriptor v2) - { return get(vpm1, v1) < get(vpm1, v2);}); - std::sort(tm2_verts.begin(), tm2_verts.end(), - [&vpm2](vertex_descriptor v1, vertex_descriptor v2) - { return get(vpm2, v1) < get(vpm2, v2);}); - - static constexpr std::size_t NOT_SHARED = std::size_t(-1); - auto is_vshared_map1 = get(CGAL::dynamic_vertex_property_t(), tm1, NOT_SHARED); - auto is_vshared_map2 = get(CGAL::dynamic_vertex_property_t(), tm2, NOT_SHARED); - std::vector > common_vertices; - - auto itv1=tm1_verts.begin(), itv1_end=tm1_verts.end(), - itv2=tm2_verts.begin(), itv2_end=tm2_verts.end(); - while (itv1!=itv1_end && itv2!=itv2_end) - { - if (get(vpm1, *itv1) < get(vpm2, *itv2)) - ++itv1; - else - { - if (!(get(vpm2, *itv2) < get(vpm1, *itv1))) - { - put(is_vshared_map1, *itv1, common_vertices.size()); - put(is_vshared_map2, *itv2, common_vertices.size()); - common_vertices.emplace_back(*itv1, *itv2); - ++itv1; - } - ++itv2; - } - } - - std::vector> identical_patches; - if (common_vertices.size()>=4) // one triangle is useless - { - auto is_fshared_map1 = get(CGAL::dynamic_face_property_t(), tm1, false); - auto is_fshared_map2 = get(CGAL::dynamic_face_property_t(), tm2, false); - auto is_eshared_map1 = get(CGAL::dynamic_edge_property_t(), tm1, false); - auto is_eshared_map2 = get(CGAL::dynamic_edge_property_t(), tm2, false); - bool identical_meshes = tm1_bb==tm2_bb && tm1_verts.size()==tm2_verts.size() && tm1_verts.size()==common_vertices.size(); - - // detect and make identical faces in tm1 and tm2 - std::size_t k0=0; - for (const std::pair& vp : common_vertices) - { - for (halfedge_descriptor h1 : halfedges_around_target(vp.first, tm1)) - { - face_descriptor f1 = face(h1, tm1); - if (f1 == graph_traits::null_face()) continue; - std::size_t k1 = get(is_vshared_map1, source(h1,tm1)); - if (k1 == NOT_SHARED) continue; - h1 = next(h1, tm1); - std::size_t k2 = get(is_vshared_map1, target(h1,tm1)); - if (k2 == NOT_SHARED) continue; - if (k0>k1 || k0>k2) continue; // report faces only once - - //looking for (k0,k1,k2) in tm2 - vertex_descriptor v0=common_vertices[k0].second, - v1=common_vertices[k1].second, - v2=common_vertices[k2].second; - auto [h2, found] = halfedge(v0,v1,tm2); - if (found) - { - if (target(next(h2, tm2), tm2)!=v2) - { - h2=opposite(h2,tm2); - if (target(next(h2, tm2), tm2)!=v2) - { - identical_meshes=false; - continue; - } - } - face_descriptor f2 = face(h2,tm2); - if (f2 == graph_traits::null_face()) - { - identical_meshes=false; - continue; - } - - identical_patches.emplace_back(f1,f2); - visitor.input_have_coplanar_faces(); - put(is_fshared_map1, f1, true); - put(is_fshared_map2, f2, true); - - for (int i=0;i<3;++i) - { - put(is_eshared_map1, edge(h1, tm1), true); - h1=next(h1,tm1); - } - for (int i=0;i<3;++i) - { - put(is_eshared_map2, edge(h2, tm2), true); - h2=next(h2,tm2); - } - } - else - identical_meshes=false; - } - ++k0; - } - - // TODO: we should also be able to fill coplanar_patches_of_tm[1/2] as we might exclude entirely some components - if (!identical_meshes) - { - visitor.start_filtering_intersections(); - filter_intersections(tm1, tm2, vpm1, vpm2, non_manifold_feature_map_2, throw_on_self_intersection, tm1_faces, tm2_faces, bb12, - is_fshared_map1, is_eshared_map2, false); - filter_intersections(tm2, tm1, vpm2, vpm1, non_manifold_feature_map_1, throw_on_self_intersection, tm2_faces, tm1_faces, bb12, - is_fshared_map2, is_eshared_map1, true); - visitor.end_filtering_intersections(); - - // dumping shared edges in output - if constexpr (!std::is_same_v) - { - typedef typename boost::property_traits::value_type Point_3; - for (edge_descriptor e : edges(tm1)) - { - if (get(is_eshared_map1, e)) - { - std::vector polyline; - polyline.push_back(get(vpm1, source(e,tm1))); - polyline.push_back(get(vpm1, target(e,tm1))); - *output++=polyline; - } - } - } - } - else - { - if constexpr (!std::is_same_v) - { - typedef typename boost::property_traits::value_type Point_3; - for (edge_descriptor e : edges(tm1)) - { - - std::vector polyline; - polyline.push_back(get(vpm1, source(e,tm1))); - polyline.push_back(get(vpm1, target(e,tm1))); - *output++=polyline; - } - } - visitor.finalize(nodes,tm1,tm2,vpm1,vpm2, identical_patches); - - return output; - } - } - else - { - Static_boolean_property_map is_fshared_map; - Static_boolean_property_map is_eshared_map; - visitor.start_filtering_intersections(); - filter_intersections(tm1, tm2, vpm1, vpm2, non_manifold_feature_map_2, throw_on_self_intersection, tm1_faces, tm2_faces, bb12, - is_fshared_map, is_eshared_map, false); - filter_intersections(tm2, tm1, vpm2, vpm1, non_manifold_feature_map_1, throw_on_self_intersection, tm2_faces, tm1_faces, bb12, - is_fshared_map, is_eshared_map, true); - visitor.end_filtering_intersections(); - } + visitor.start_filtering_intersections(); + filter_intersections(tm1, tm2, vpm1, vpm2, non_manifold_feature_map_1, non_manifold_feature_map_2, throw_on_self_intersection, tm1_faces, tm2_faces, tm1_bb, tm2_bb); + visitor.end_filtering_intersections(); Node_id current_node((std::numeric_limits::max)()); CGAL_assertion(current_node+1==0); -// TODO: handle non-manifold edges in coplanar + // TODO: handle non-manifold edges in coplanar #ifndef DO_NOT_HANDLE_COPLANAR_FACES //first handle coplanar triangles if (&tm1<&tm2) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/connected_components.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/connected_components.h index a389e983957b..667d6139cfb9 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/connected_components.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/connected_components.h @@ -216,34 +216,46 @@ connected_components(const PolygonMesh& pmesh, EdgeConstraintMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained)); - typedef typename GetInitializedFaceIndexMap::const_type FaceIndexMap; - FaceIndexMap fimap = get_initialized_face_index_map(pmesh, np); - + // typedef typename GetInitializedFaceIndexMap::const_type FaceIndexMap; + // FaceIndexMap fimap = get_initialized_face_index_map(pmesh, np); + + // std::vector handled(num_faces(pmesh), false); + // auto is_handled = [&](face_descriptor f){ + // return handled[get(fimap,f)]; + // }; + // auto set_handled = [&](face_descriptor f){ + // handled[get(fimap,f)] = true; + // }; + + using F2B_tag = typename CGAL::dynamic_face_property_t; + using Face_to_bool_map = typename boost::property_map::const_type; + Face_to_bool_map handled = get(F2B_tag(), pmesh, false); + auto is_handled = [&](face_descriptor f){ + return get(handled, f); + }; + auto set_handled = [&](face_descriptor f){ + put(handled, f, true); + }; typename boost::property_traits::value_type i=0; - std::vector handled(num_faces(pmesh), false); for (face_descriptor f : faces(pmesh)) { - if (handled[get(fimap,f)]) continue; + if (is_handled(f)) continue; std::vector queue; queue.push_back(f); while(!queue.empty()) { face_descriptor fq = queue.back(); queue.pop_back(); - typename boost::property_traits::value_type fq_id = get(fimap,fq); - if ( handled[fq_id]) continue; - handled[fq_id]=true; + if (is_handled(fq)) continue; + set_handled(fq); put(fcm, fq, i); for (halfedge_descriptor h : halfedges_around_face(halfedge(fq, pmesh), pmesh)) { if ( get(ecmap, edge(h, pmesh)) ) continue; halfedge_descriptor opp = opposite(h, pmesh); face_descriptor fqo = face(opp, pmesh); - if ( fqo != GT::null_face() ) - { - if ( !handled[get(fimap,fqo)] ) + if ( fqo != GT::null_face() && !is_handled(fqo)) queue.push_back(fqo); - } } } ++i; diff --git a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h index 148785e4d90f..0612020b68f4 100644 --- a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h +++ b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -48,6 +49,10 @@ #include #include +#ifdef CGAL_OUTPUT_BUILDER_RUNNING_TIME +#include +#endif + namespace CGAL { #ifndef DOXYGEN_RUNNING @@ -1300,6 +1305,453 @@ class Surface_mesh return true; } + template < bool reverse_patch_orientation, + class PatchDescription, + class VertexPointMap, + class VertexPointMapOut, + // class EdgeMarkMapOut, + // class EdgeMarkMapIn , + class VertexToVertexMap, + class EdgeToEdgeMap , + class UserVisitor> + void append_patch( + PatchDescription& patch, + const Surface_mesh &tm, + const VertexPointMapOut& vpm_out, + const VertexPointMap& vpm_tm, + // EdgeMarkMapOut& edge_mark_map_out, + // const EdgeMarkMapIn& edge_mark_map_in, + VertexToVertexMap& tm_to_output_vertices, + EdgeToEdgeMap& tm_to_output_edges, + UserVisitor& user_visitor, + std::size_t vertices_idx_begin, + std::size_t edges_idx_begin, + std::size_t face_idx_begin) + { + Self &output = *this; + + using vertex_descriptor = SM_Vertex_index; + using edge_descriptor = SM_Edge_index; + using halfedge_descriptor = SM_Halfedge_index; + using face_descriptor = SM_Face_index; + + auto get_halfedge = [&](halfedge_descriptor h){ + edge_descriptor e = get(tm_to_output_edges, tm.edge(h)); + halfedge_descriptor h_out = output.halfedge(e); + if( output.target(h_out) == get(tm_to_output_vertices, tm.target(h))){ + CGAL_assertion( source(h_out, output) == get(tm_to_output_vertices, source(h, tm)) ); + if constexpr(reverse_patch_orientation) + return output.opposite(h_out); + else + return h_out; + } + CGAL_assertion( target(h_out, output) == get(tm_to_output_vertices, source(h, tm)) ); + CGAL_assertion( source(h_out, output) == get(tm_to_output_vertices, target(h, tm)) ); + if constexpr(reverse_patch_orientation) + return h_out; + else + return output.opposite(h_out); + }; + + #ifdef CGAL_COREFINEMENT_POLYHEDRA_DEBUG + #warning the size of tm_to_output_edges will increase at each step \ + when adding new patches by the size of internal edges. \ + Maybe the use of a copy would be better since we do not need \ + the internal edges added? + #endif + + // Fill Vertex to vertex map and vpm for interior vertices + auto fill_vertex = [&](std::size_t i){ + vertex_descriptor v = patch.interior_vertices[i]; + vertex_descriptor new_v(vertices_idx_begin + i); + + put(tm_to_output_vertices, v, new_v); + put(vpm_out, new_v, get(vpm_tm, v)); + output.set_halfedge(new_v, null_halfedge()); + }; +#ifdef CGAL_LINKED_WITH_TBB + if constexpr(true){ + if(patch.interior_vertices.size() > 10000){ + tbb::parallel_for(std::size_t(0), patch.interior_vertices.size(), fill_vertex); + } else { + for(std::size_t i=0; i(ed, new_edge, + // edge_mark_map_in, edge_mark_map_out); + + CGAL_assertion(is_border(new_h, output)); + CGAL_assertion(is_border(opposite(new_h, output), output)); + + // Fill the halfedge if interior vertices + if (tm.halfedge(tm.target(h)) == h && + output.halfedge(output.target(new_h)) == null_halfedge()){ + user_visitor.before_vertex_copy(tm.target(h), tm, output); // Call visitor when the vertex would be complete + output.set_halfedge(output.target(new_h), new_h); + user_visitor.after_vertex_copy(tm.target(h), tm, output.target(new_h), output); + } + + if (tm.halfedge(tm.source(h)) == tm.opposite(h) && + output.halfedge(output.source(new_h)) == null_halfedge()){ + user_visitor.before_vertex_copy(tm.source(h), tm, output); // Call visitor when the vertex would be complete + output.set_halfedge(output.source(new_h), new_h); + user_visitor.after_vertex_copy(tm.source(h), tm, output.target(new_h), output); + } + }; +#ifdef CGAL_LINKED_WITH_TBB + if constexpr(true){ + if(patch.interior_edges.size() > 10000){ + tbb::parallel_for(std::size_t(0), patch.interior_edges.size(), fill_edge); + } else { + for(std::size_t i=0; i hedges = { get_halfedge(h_in_1), get_halfedge(h_in_2), get_halfedge(h_in_3) }; + + user_visitor.before_face_copy(f, tm, output); + SM_Face_index new_f(face_idx_begin + i); + user_visitor.after_face_copy(f, tm, new_f, output); + output.set_halfedge(new_f, hedges[0]); + + for (int i=0;i<3;++i) + { + CGAL_assertion(hedges[i] != null_halfedge()); + if(reverse_patch_orientation) + output.set_next(hedges[i], hedges[(i+2)%3]); + else + output.set_next(hedges[i], hedges[(i+1)%3]); + output.set_face(hedges[i], new_f); + } + }; +#ifdef CGAL_LINKED_WITH_TBB + if constexpr(true){ + if(patch.faces.size() > 10000){ + tbb::parallel_for(std::size_t(0), patch.faces.size(), fill_face); + } else { + for(std::size_t i=0; i + typename boost::graph_traits::halfedge_descriptor + next_marked_halfedge_around_target_vertex( + typename boost::graph_traits::halfedge_descriptor h, + const PolygonMesh& pm, + const MarkedEdgeSet& marked_edges) + { + CGAL_assertion( marked_edges.count(pm.edge(h))!= 0 ); + typename boost::graph_traits::halfedge_descriptor nxt = + pm.next(h); + while( !marked_edges.count(pm.edge(nxt)) ) + { + nxt=pm.next(pm.opposite(nxt)); + } + CGAL_assertion(nxt!=h); + return nxt; + } + + + template + void import_polyline( + SM_Halfedge_index h1, + SM_Halfedge_index h2, + const Self& tm1, + const Self& tm2, + std::size_t nb_segments, + VertexToVertexMap& pm1_to_output_vertices, + VertexToVertexMap& pm2_to_output_vertices, + EdgeToEdgeMap& pm1_to_output_edges, + EdgeToEdgeMap& pm2_to_output_edges, + const IntersectionEdgeMap& intersection_edges1, + const IntersectionEdgeMap& intersection_edges2, + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2, + const VertexPointMapOut& vpm_out, + std::vector& output_shared_edges, + UserVisitor& user_visitor) + { + using TriangleMesh = Self; + using vertex_descriptor = SM_Vertex_index; + using halfedge_descriptor = SM_Halfedge_index; + using face_descriptor = SM_Face_index; + + auto set_output_vertex = [&](vertex_descriptor v1, vertex_descriptor v2, halfedge_descriptor h_out){ + user_visitor.before_vertex_copy(v1, tm1, *this); + vertex_descriptor new_v = this->add_vertex(); + this->set_halfedge(new_v, this->opposite(h_out)); + put(vpm_out, new_v, get(vpm1, v1)); + user_visitor.after_vertex_copy(v1, tm1, new_v, *this); + put(pm1_to_output_vertices, v1, new_v); + put(pm2_to_output_vertices, v2, new_v); + return new_v; + }; + + output_shared_edges.push_back(this->edge(this->add_edge())); + halfedge_descriptor h_out = this->halfedge(output_shared_edges.back()); + + // make sure the first vertex does not already exist + vertex_descriptor src; + if( get(pm1_to_output_vertices, tm1.source(h1)) == null_vertex()) + src = set_output_vertex(tm1.source(h1), tm2.source(h2), h_out); + else + src = get(pm1_to_output_vertices, tm1.source(h1)); + + //make sure the target vertex does not already exist if it is a polyline endpoint + vertex_descriptor tgt; + if ( nb_segments==1 ){ + if( get(pm1_to_output_vertices, tm1.target(h1)) == null_vertex()) + tgt = set_output_vertex(tm1.target(h1), tm2.target(h2), h_out); + else + tgt = get(pm1_to_output_vertices, tm1.target(h1)); + } else { + tgt = set_output_vertex(tm1.target(h1), tm2.target(h2), h_out); + } + + // update source and target vertex of the edge created + this->set_target(h_out, tgt); + this->set_target(this->opposite(h_out), src); + + //set the correspondence + put(pm1_to_output_edges, tm1.edge(h1), this->edge(h_out)); + put(pm2_to_output_edges, tm2.edge(h2), this->edge(h_out)); + + user_visitor.intersection_edge_copy(h1, tm1, h2, tm2, h_out, *this); + + src=tgt; + for (std::size_t i=1; iedge(this->add_edge())); + h_out = this->halfedge(output_shared_edges.back()); + //get the new edge + h1 = next_marked_halfedge_around_target_vertex(h1, tm1, intersection_edges1); + h2 = next_marked_halfedge_around_target_vertex(h2, tm2, intersection_edges2); + + user_visitor.intersection_edge_copy(h1, tm1, h2, tm2, h_out, *this); + + //if this is the final segment, only create a target vertex if it does not exist + if ( i+1 == nb_segments ){ + if( get(pm1_to_output_vertices, tm1.target(h1)) == null_vertex()) + tgt = set_output_vertex(tm1.target(h1), tm2.target(h2), h_out); + else + tgt = get(pm1_to_output_vertices, tm1.target(h1)); + } else { + tgt = set_output_vertex(tm1.target(h1), tm2.target(h2), h_out); + } + + this->set_target(h_out, tgt); + this->set_target(this->opposite(h_out), src); + + src = tgt; + + put(pm1_to_output_edges, tm1.edge(h1), this->edge(h_out)); + put(pm2_to_output_edges, tm2.edge(h2), this->edge(h_out)); + } + } + + // Version of filll_new_triangle_mesh of face_graph_utils spefically written for Surface_mesh + template < bool reverse_orientation_of_patches_from_tm1, + bool reverse_orientation_of_patches_from_tm2, + class IntersectionEdgeMap, + class VertexPointMap1, + class VertexPointMap2, + class VertexPointMapOut, + class EdgeMarkMap1, + class EdgeMarkMap2, + class EdgeMarkMapOut, + class IntersectionPolylines, + class PatchContainer1, + class PatchContainer2, + class UserVisitor> + void fill_new_triangle_mesh( + const boost::dynamic_bitset<>& patches_of_tm1_to_import, + const boost::dynamic_bitset<>& patches_of_tm2_to_import, + PatchContainer1& patches_of_tm1, + PatchContainer2& patches_of_tm2, + const IntersectionPolylines& polylines, + const IntersectionEdgeMap& intersection_edges1, + const IntersectionEdgeMap& intersection_edges2, + const VertexPointMap1& vpm1, + const VertexPointMap2& vpm2, + const VertexPointMapOut& vpm_out, + const EdgeMarkMap1& edge_mark_map1, + const EdgeMarkMap2& edge_mark_map2, + EdgeMarkMapOut& edge_mark_map_out, + std::vector< SM_Edge_index >& output_shared_edges, + UserVisitor& user_visitor) + { + using TriangleMesh = Self; + using Point = typename VertexPointMap1::value_type; + using vertex_descriptor = SM_Vertex_index; + using edge_descriptor = SM_Edge_index; + using face_descriptor = SM_Face_index; + + using V2V_tag = typename CGAL::dynamic_vertex_property_t; + using Vertex_to_vertex_map = typename boost::property_map::const_type; + + using E2E_tag = typename CGAL::dynamic_edge_property_t; + using Edge_to_edge_map = typename boost::property_map::const_type; + + const TriangleMesh &tm1 = patches_of_tm1.pm; + const TriangleMesh &tm2 = patches_of_tm2.pm; + + Vertex_to_vertex_map tm1_to_output_vertices = get(V2V_tag(), tm1, null_vertex()); + Vertex_to_vertex_map tm2_to_output_vertices = get(V2V_tag(), tm2, null_vertex()); + Edge_to_edge_map tm1_to_output_edges = get(E2E_tag(), tm1, null_edge()); + Edge_to_edge_map tm2_to_output_edges = get(E2E_tag(), tm2, null_edge()); + + #ifdef CGAL_OUTPUT_BUILDER_RUNNING_TIME + Real_timer t; + t.start(); + #endif + + // Import polylines + // TODO Fully sequential, a parallelization is possible but for most example not worth it + output_shared_edges.reserve( std::accumulate(polylines.lengths.begin(), polylines.lengths.end(), std::size_t(0)) ); + size_type nb_polylines = polylines.lengths.size(); + for (size_type i=0; i < nb_polylines; ++i) + if (!polylines.to_skip.test(i)) + import_polyline(polylines.tm1[i], polylines.tm2[i], + tm1, tm2, + polylines.lengths[i], + tm1_to_output_vertices, tm2_to_output_vertices, + tm1_to_output_edges, tm2_to_output_edges, + intersection_edges1, intersection_edges2, + vpm1, vpm2, vpm_out, + output_shared_edges, + user_visitor); + + #ifdef CGAL_OUTPUT_BUILDER_RUNNING_TIME + std::cout << "Import polyline " << t.time() << std::endl; + #endif + + // Get ids of patch to append + std::vector ids_of_patches_to_append_from_tm1; + std::vector ids_of_patches_to_append_from_tm2; + ids_of_patches_to_append_from_tm1.reserve(patches_of_tm1_to_import.count()); + for (std::size_t i= patches_of_tm1_to_import.find_first(); + i < patches_of_tm1_to_import.npos; + i = patches_of_tm1_to_import.find_next(i)){ + ids_of_patches_to_append_from_tm1.push_back(i); + } + ids_of_patches_to_append_from_tm2.reserve(patches_of_tm2_to_import.count()); + for (std::size_t i= patches_of_tm2_to_import.find_first(); + i < patches_of_tm2_to_import.npos; + i = patches_of_tm2_to_import.find_next(i)){ + ids_of_patches_to_append_from_tm2.push_back(i); + } + + // Compute final sizes + size_type nv = num_vertices(), ne = num_edges(), nf = num_faces(); + size_type tnv = num_vertices(), tne = num_edges(), tnf = num_faces(); + for (std::size_t i : ids_of_patches_to_append_from_tm1){ + tnv += patches_of_tm1[i].interior_vertices.size(); + tne += patches_of_tm1[i].interior_edges.size(); + tnf += patches_of_tm1[i].faces.size(); + } + for (std::size_t i : ids_of_patches_to_append_from_tm2){ + tnv += patches_of_tm2[i].interior_vertices.size(); + tne += patches_of_tm2[i].interior_edges.size(); + tnf += patches_of_tm2[i].faces.size(); + } + resize(tnv, tne, tnf); + + // Append patches + for (std::size_t i : ids_of_patches_to_append_from_tm1){ + append_patch( + patches_of_tm1[i], + tm1, + vpm_out, + vpm1, + tm1_to_output_vertices, + tm1_to_output_edges, + user_visitor, + nv, ne, nf); + nv += patches_of_tm1[i].interior_vertices.size(); + ne += patches_of_tm1[i].interior_edges.size(); + nf += patches_of_tm1[i].faces.size(); + } + + for (std::size_t i : ids_of_patches_to_append_from_tm2){ + append_patch( + patches_of_tm2[i], + tm2, + vpm_out, + vpm2, + tm2_to_output_vertices, + tm2_to_output_edges, + user_visitor, + nv, ne, nf); + nv += patches_of_tm2[i].interior_vertices.size(); + ne += patches_of_tm2[i].interior_edges.size(); + nf += patches_of_tm2[i].faces.size(); + } + + #ifdef CGAL_OUTPUT_BUILDER_RUNNING_TIME + std::cout << "Import patch: " << t.time() << std::endl; + #endif + } + ///@}