From 60567eccbbe95f670a74b7d1841e26fa4ea32cd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 28 Sep 2021 16:25:43 +0200 Subject: [PATCH 1/2] Add Kd_tree-based vertex-vertex snapping + PMP::snap visitor --- .../internal/Snapping/helper.h | 127 ++- .../internal/Snapping/snap.h | 223 +++-- .../internal/Snapping/snap_vertices.h | 812 +++++++++++++----- 3 files changed, 885 insertions(+), 277 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/helper.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/helper.h index 5352ca1a8dd..0c3fc465392 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/helper.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/helper.h @@ -37,7 +37,8 @@ void vertices_as_halfedges(const VertexRange& vertex_range, *out++ = halfedge(v, pmesh); } -// Assigns at each vertex the 'tolerance' value as tolerance, but bounded by a percentage of the length of its shortest incident edge +// Assigns at each vertex the 'tolerance' value as tolerance, +// but bounded by a percentage of the length of its shortest incident edge template = sq_va_l * sq_vb_l * sq_cos); } +template +struct Snapping_default_visitor +{ + bool m_stop = false; + + bool stop() { return m_stop; } + + // ------------------------------- Preprocessing ------------------------------ + + // Called when starting preprocessing phase of the mesh. + void start_mesh_preprocessing() { } + + // Called at the end of the preprocessing phase of the mesh. + void end_mesh_preprocessing() { } + + // ------------------------------- Vertex - Vertex ------------------------------- + + // Called when starting snapping a range of vertices of `A` and a range of vertices of `B`. + void start_vertex_vertex_phase() { } + + // Called when finished snapping a range of vertices of `A` and a range of vertices of `B`. + void end_vertex_vertex_phase() { } + + // Called when trying to find `v`-`v_other` pair (with `v` in `A` and `v_other` in `B`) + // + // Available only if vertex-vertex pair detection is based on the kd-tree. + // Must be threadsafe if parallelism is used. + template + void on_vertex_vertex_inquiry(const Vertex /*v*/, const Mesh&) { } + + // Called when a match between two ranges of vertices is found. + // All vertices in `va` have the same position (and same for `vb`). + // + // Must be threadsafe if parallelism is used. + template + void on_vertex_vertex_match(const VertexContainer& /*va*/, const Mesh&, + const VertexContainer& /*vb*/, const Mesh&) { } + + // Called before proceeding with actual vertex-vertex snapping. + void start_vertex_vertex_snapping() { } + + // Called after proceeding with actual vertex-vertex snapping. + void end_vertex_vertex_snapping() { } + + // Called before two ranges of vertices are snapped together. + // All vertices in `va` have the same position (and same for `vb`). + template + void before_vertex_vertex_snap(const VertexContainer& /*va*/, const Mesh&, + const VertexContainer& /*vb*/, const Mesh&) { } + + // Called after two ranges of vertices have been snapped. + template + void after_vertex_vertex_snap(const VertexContainer& /*va*/, const Mesh&, + const VertexContainer& /*vb*/, const Mesh&) { } + + // ------------------------------- Vertex - Edge ------------------------------- + + // Called when starting snapping a range of vertices of A onto edges of B + void start_first_vertex_edge_phase() { } + + // Called when finished snapping a range of vertices of A onto edges of B + void end_first_vertex_edge_phase() { } + + // Called when starting snapping a range of vertices of B onto edges of A. + // Not called if A and B are the same mesh. + void start_second_vertex_edge_phase() { } + + // Called when starting snapping a range of vertices of B onto edges of A. + // Not called if A and B are the same mesh. + void end_second_vertex_edge_phase() { } + + // Called when trying to find `v`-`edge` pair + // + // Available only if vertex-vertex pair detection is based on the kd-tree. + // Must be threadsafe if parallelism is used. + template + void on_vertex_edge_inquiry(const Vertex /*v*/, const Mesh& /*tm*/) { } + + // Called when a match between a vertex and an edge is found. + // + // Must be threadsafe if parallelism is used. + template + void on_vertex_edge_match(const Vertex, const Mesh&, const Edge, const Mesh&, const Point&) { } + + // Called before proceeding with actual snapping. + void start_vertex_edge_snapping() { } + + // Called after proceeding with actual snapping. + void end_vertex_edge_snapping() { } + + // Called before a new vertex is created. + template + void before_vertex_edge_snap(const Edge, const Mesh&, const Point&) { } + + // Called after a new vertex has been created, splitting an edge. + template + void after_vertex_edge_snap(const Vertex /*new_vertex*/, const Mesh&) { } + + // ------------------------------- Two passes (segmentation or not) ------------------------------ + + // Called at the start of the snapping pass that is restricted to compatible patch (first pass). + void start_patch_snapping() { } + + // Called at the end of the snapping pass that is restricted to compatible patch (first pass). + void end_patch_snapping() { } + + // Called at the start of the second snapping pass, all elements are potentially compatible (second pass). + void start_global_snapping() { } + + // Called at the end of the second snapping pass, all elements are potentially compatible (second pass). + void end_global_snapping() { } +}; + } // namespace internal } // namespace Polygon_mesh_processing } // namespace CGAL diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h index 3571ae48b52..5f5e807d49e 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h @@ -453,7 +453,8 @@ template + typename VPMS, typename VPMT, typename GT, + typename Visitor> void find_splittable_edge(const VertexWithTolerance& vertex_with_tolerance, EdgeToSplitMap& edges_to_split, const AABBTree* aabb_tree_ptr, @@ -463,7 +464,8 @@ void find_splittable_edge(const VertexWithTolerance& vertex_with_tolerance, const TriangleMesh& tm_T, FacePatchMap_T face_patch_map_T, VPMT vpm_T, - const GT& gt) + const GT& gt, + Visitor& visitor) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; @@ -477,6 +479,9 @@ void find_splittable_edge(const VertexWithTolerance& vertex_with_tolerance, typedef internal::Projection_traits Projection_traits; + if(visitor.stop()) + throw CGAL::internal::Throw_at_output_exception(); + // by construction the whole range has the same position const halfedge_descriptor h = vertex_with_tolerance.first; const vertex_descriptor v = target(h, tm_S); @@ -488,6 +493,8 @@ void find_splittable_edge(const VertexWithTolerance& vertex_with_tolerance, std::cout << "--------------------------- Query: " << v << " (" << query << ")" << std::endl; #endif + visitor.on_vertex_edge_inquiry(v, tm_S); + Projection_traits traversal_traits(h, patch_id, sq_tolerance, tm_S, vpm_S, tm_T, face_patch_map_T, vpm_T, aabb_tree_ptr->traits()); aabb_tree_ptr->traversal(query, traversal_traits); @@ -522,7 +529,6 @@ void find_splittable_edge(const VertexWithTolerance& vertex_with_tolerance, if(!is_close_enough) return; - CGAL_assertion(get(vpm_T, source(closest_e, tm_T)) != query && get(vpm_T, target(closest_e, tm_T)) != query); @@ -534,43 +540,17 @@ void find_splittable_edge(const VertexWithTolerance& vertex_with_tolerance, // Using a map because we need to know if the same halfedge is split multiple times Edges_to_split_map_inserter()(edges_to_split, closest_h, vertex_with_tolerance.first, closest_p); + + visitor.on_vertex_edge_match(v, tm_S, closest_h, tm_T, closest_p); } #ifdef CGAL_LINKED_WITH_TBB template -struct Find_splittable_edge_for_parallel_for + typename VPMS, typename VPMT, typename GT, typename Visitor> +struct Parallel_find_splittable_edge { - Find_splittable_edge_for_parallel_for(const PointWithToleranceContainer& points_with_tolerance, - EdgeToSplitMap& edges_to_split, - const AABBTree* aabb_tree_ptr, - const TriangleMesh& tm_S, - const VertexPatchMap_S vertex_patch_map_S, - const VPMS vpm_S, - const TriangleMesh& tm_T, - const FacePatchMap_T face_patch_map_T, - const VPMT vpm_T, - const GT& gt) - : - m_points_with_tolerance(points_with_tolerance), - m_edges_to_split(edges_to_split), m_aabb_tree_ptr(aabb_tree_ptr), - m_tm_S(tm_S), m_vertex_patch_map_S(vertex_patch_map_S), m_vpm_S(vpm_S), - m_tm_T(tm_T), m_face_patch_map_T(face_patch_map_T), m_vpm_T(vpm_T), - m_gt(gt) - { } - - void operator()(const tbb::blocked_range& r) const - { - for(std::size_t i=r.begin(); i!=r.end(); ++i) - { - find_splittable_edge(m_points_with_tolerance[i], m_edges_to_split, m_aabb_tree_ptr, - m_tm_S, m_vertex_patch_map_S, m_vpm_S, - m_tm_T, m_face_patch_map_T, m_vpm_T, m_gt); - } - } - private: const PointWithToleranceContainer& m_points_with_tolerance; EdgeToSplitMap& m_edges_to_split; @@ -582,18 +562,51 @@ private: const FacePatchMap_T m_face_patch_map_T; const VPMT m_vpm_T; const GT& m_gt; + Visitor& m_visitor; + +public: + Parallel_find_splittable_edge(const PointWithToleranceContainer& points_with_tolerance, + EdgeToSplitMap& edges_to_split, + const AABBTree* aabb_tree_ptr, + const TriangleMesh& tm_S, + const VertexPatchMap_S vertex_patch_map_S, + const VPMS vpm_S, + const TriangleMesh& tm_T, + const FacePatchMap_T face_patch_map_T, + const VPMT vpm_T, + const GT& gt, + Visitor& visitor) + : + m_points_with_tolerance(points_with_tolerance), + m_edges_to_split(edges_to_split), m_aabb_tree_ptr(aabb_tree_ptr), + m_tm_S(tm_S), m_vertex_patch_map_S(vertex_patch_map_S), m_vpm_S(vpm_S), + m_tm_T(tm_T), m_face_patch_map_T(face_patch_map_T), m_vpm_T(vpm_T), + m_gt(gt), m_visitor(visitor) + { } + + void operator()(const tbb::blocked_range& r) const + { + for(std::size_t i=r.begin(); i!=r.end(); ++i) + { + find_splittable_edge(m_points_with_tolerance[i], m_edges_to_split, m_aabb_tree_ptr, + m_tm_S, m_vertex_patch_map_S, m_vpm_S, + m_tm_T, m_face_patch_map_T, m_vpm_T, m_gt, m_visitor); + } + } }; #endif template + typename TriangleMesh, + typename VPMS, typename VPMT, typename GeomTraits, + typename Visitor> std::size_t split_edges(EdgesToSplitContainer& edges_to_split, TriangleMesh& tm_S, VPMS vpm_S, TriangleMesh& tm_T, VPMT vpm_T, const GeomTraits& gt, + Visitor& visitor, const bool is_source_mesh_fixed) // when snapping is B --> A and the mesh B is fixed { #ifdef CGAL_PMP_SNAP_DEBUG @@ -611,10 +624,14 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split, typedef std::vector Vertices_with_new_position; typedef std::pair Edge_and_splitters; - std::size_t snapped_n = 0; - +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME CGAL::Real_timer timer; timer.start(); +#endif + + visitor.start_vertex_edge_snapping(); + + std::size_t snapped_n = 0; for(Edge_and_splitters& es : edges_to_split) { @@ -642,6 +659,9 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split, Point previous_split_position = get(vpm_S, *(vertices(tm_S).begin())); // dummy value to avoid "used uninitialized" warning for(const Vertex_with_new_position& vnp : splitters) { + if(visitor.stop()) + return snapped_n; + const halfedge_descriptor splitter_h = vnp.first; const vertex_descriptor splitter_v = target(splitter_h, tm_S); const Point new_position = is_source_mesh_fixed ? get(vpm_S, splitter_v) : vnp.second; @@ -668,9 +688,13 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split, vertex_descriptor new_v = boost::graph_traits::null_vertex(); if(do_split) { + visitor.before_vertex_edge_snap(h_to_split, tm_T, vnp.second); + CGAL::Euler::split_edge(h_to_split, tm_T); new_v = source(h_to_split, tm_T); put(vpm_T, new_v, new_position); // position of the new point on the target mesh + + visitor.after_vertex_edge_snap(new_v, tm_T); } if(!is_source_mesh_fixed) @@ -740,6 +764,13 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split, } } + visitor.end_vertex_edge_snapping(); + +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + timer.stop(); + std::cout << "time for actual snapping (vertex-edge) " << timer.time() << " s." << std::endl; +#endif + return snapped_n; } @@ -757,6 +788,7 @@ template std::size_t snap_non_conformal_one_way(const HalfedgeRange& halfedge_range_S, TriangleMesh& tm_S, @@ -768,6 +800,7 @@ std::size_t snap_non_conformal_one_way(const HalfedgeRange& halfedge_range_S, FacePatchMap_T face_patch_map_T, LockedHalfedgeMap locked_halfedges_T, const bool is_source_mesh_fixed, + Visitor& visitor, const SourceNamedParameters& snp, const TargetNamedParameters& tnp) { @@ -839,7 +872,7 @@ std::size_t snap_non_conformal_one_way(const HalfedgeRange& halfedge_range_S, for(auto& p : vertices_to_snap) p.second = point_tolerance_map[get(vpm_S, target(p.first, tm_S))]; - // Since we're inserting primitives one by one, we can't pass this shared data in the constructor of the tree + // Since primitives are inserted one by one, the shared data cannot be in the constructor of the tree AABB_Traits aabb_traits; aabb_traits.set_shared_data(tm_T, vpm_T); AABB_tree aabb_tree(aabb_traits); @@ -863,10 +896,15 @@ std::size_t snap_non_conformal_one_way(const HalfedgeRange& halfedge_range_S, std::cout << "Collect edges to split with " << vertices_to_snap.size() << " vertices" << std::endl; #endif +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + CGAL::Real_timer timer; + timer.start(); +#endif + #ifndef CGAL_LINKED_WITH_TBB CGAL_static_assertion_msg (!(std::is_convertible::value), "Parallel_tag is enabled but TBB is unavailable."); -#else +#else // CGAL_LINKED_WITH_TBB if(std::is_convertible::value) { #ifdef CGAL_PMP_SNAP_DEBUG @@ -875,22 +913,45 @@ std::size_t snap_non_conformal_one_way(const HalfedgeRange& halfedge_range_S, typedef tbb::concurrent_hash_map Concurrent_edge_to_split_container; - typedef internal::Find_splittable_edge_for_parallel_for< + typedef internal::Parallel_find_splittable_edge< Vertices_with_tolerance, TriangleMesh, Concurrent_edge_to_split_container, AABB_tree, - VertexPatchMap_S, FacePatchMap_T, VPMS, VPMT, GT> Functor; - - CGAL::Real_timer timer; - timer.start(); + VertexPatchMap_S, FacePatchMap_T, VPMS, VPMT, GT, Visitor> Functor; Concurrent_edge_to_split_container edges_to_split; Functor f(vertices_to_snap, edges_to_split, &aabb_tree, - tm_S, vertex_patch_map_S, vpm_S, tm_T, face_patch_map_T, vpm_T, gt); - tbb::parallel_for(tbb::blocked_range(0, vertices_to_snap.size()), f); + tm_S, vertex_patch_map_S, vpm_S, tm_T, face_patch_map_T, vpm_T, + gt, visitor); - std::cout << "Time to gather edges: " << timer.time() << std::endl; + try + { + tbb::parallel_for(tbb::blocked_range(0, vertices_to_snap.size()), f); + } + catch(const CGAL::internal::Throw_at_output_exception&) + { + return 0; + } +#if TBB_USE_CAPTURED_EXCEPTION + catch(const tbb::captured_exception& e) + { + const std::string tn1(e.name()); + const std::string tn2(typeid(const CGAL::internal::Throw_at_output_exception&).name()); + if(tn1 != tn2) + { + std::cerr << "Unexpected throw: " << tn1 << std::endl; + throw; + } - return split_edges(edges_to_split, tm_S, vpm_S, tm_T, vpm_T, gt, is_source_mesh_fixed); + return 0; + } +#endif + +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + timer.stop(); + std::cout << "time for find split edges (parallel): " << timer.time() << std::endl; +#endif + + return split_edges(edges_to_split, tm_S, vpm_S, tm_T, vpm_T, gt, visitor, is_source_mesh_fixed); } else #endif // CGAL_LINKED_WITH_TBB @@ -900,14 +961,28 @@ std::size_t snap_non_conformal_one_way(const HalfedgeRange& halfedge_range_S, #endif std::map edges_to_split; - for(const Vertex_with_tolerance& vt : vertices_to_snap) + + try { - internal::find_splittable_edge(vt, edges_to_split, &aabb_tree, - tm_S, vertex_patch_map_S, vpm_S, - tm_T, face_patch_map_T, vpm_T, gt); + for(const Vertex_with_tolerance& vt : vertices_to_snap) + { + internal::find_splittable_edge(vt, edges_to_split, &aabb_tree, + tm_S, vertex_patch_map_S, vpm_S, + tm_T, face_patch_map_T, vpm_T, + gt, visitor); + } + } + catch(const CGAL::internal::Throw_at_output_exception&) + { + return 0; } - return split_edges(edges_to_split, tm_S, vpm_S, tm_T, vpm_T, gt, is_source_mesh_fixed); +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + timer.stop(); + std::cout << "time for find split edges (sequential): " << timer.time() << std::endl; +#endif + + return split_edges(edges_to_split, tm_S, vpm_S, tm_T, vpm_T, gt, visitor, is_source_mesh_fixed); } } @@ -999,15 +1074,28 @@ std::size_t snap_non_conformal(HalfedgeRange& halfedge_range_A, internal_np::face_patch_t, NamedParameters_B, Constant_property_map /*default*/ >::type Face_patch_map_B; + typedef typename internal_np::Lookup_named_param_def < + internal_np::visitor_t, + NamedParameters_A, + internal::Snapping_default_visitor // default + >::reference Visitor; + using CGAL::parameters::choose_parameter; using CGAL::parameters::is_default_parameter; using CGAL::parameters::get_parameter; + using CGAL::parameters::get_parameter_reference; const bool is_same_mesh = (&tm_A == &tm_B); const bool simplify_first_mesh = choose_parameter(get_parameter(np_A, internal_np::do_simplify_border), false); const bool simplify_second_mesh = choose_parameter(get_parameter(np_B, internal_np::do_simplify_border), false); const bool is_second_mesh_fixed = choose_parameter(get_parameter(np_B, internal_np::do_lock_mesh), false); + internal::Snapping_default_visitor default_visitor; + Visitor& visitor = choose_parameter(get_parameter_reference(np_A, internal_np::visitor), default_visitor); + + if(visitor.stop()) + return 0; + // vertex-vertex and vertex-edge snapping is only considered within compatible patches Face_patch_map_A face_patch_map_A = choose_parameter(get_parameter(np_A, internal_np::face_patch), Constant_property_map(-1)); @@ -1053,6 +1141,9 @@ std::size_t snap_non_conformal(HalfedgeRange& halfedge_range_A, /// #1 and #1bis (Simplification of borders) ////////////////////////////////////////////////////////////////////////////////////////////////// + if(visitor.stop()) + return 0; + #ifdef CGAL_PMP_SNAP_DEBUG std::cout << "Simplify ranges (" << simplify_first_mesh << " " << simplify_second_mesh << ")..." << std::endl; #endif @@ -1079,6 +1170,8 @@ std::size_t snap_non_conformal(HalfedgeRange& halfedge_range_A, /// #2 (Two-way vertex-vertex snapping) ////////////////////////////////////////////////////////////////////////////////////////////////// + std::size_t snapped_n = 0; + // We keep in memory pairs of source/target edges that are stitchable after vertex-vertex snapping // --> these halfedges should not be considered as targets in non-conformal snapping // Similarly, matching vertices whose incident edges have matching directions are also locked @@ -1087,8 +1180,6 @@ std::size_t snap_non_conformal(HalfedgeRange& halfedge_range_A, Locked_halfedges locked_halfedges_A = get(Halfedge_bool_tag(), tm_A); Locked_halfedges locked_halfedges_B = get(Halfedge_bool_tag(), tm_B); - std::size_t snapped_n = 0; - std::vector > locked_vertices; std::vector locked_halfedges_A_vector, locked_halfedges_B_vector; @@ -1130,17 +1221,24 @@ std::size_t snap_non_conformal(HalfedgeRange& halfedge_range_A, #endif ////////////////////////////////////////////////////////////////////////////////////////////////// - /// #3 (Two one-way vertex-edge snapping) + /// #3 (Two one-way vertex-edge snappings) ////////////////////////////////////////////////////////////////////////////////////////////////// + if(visitor.stop()) + return snapped_n; + #ifdef CGAL_PMP_SNAP_DEBUG std::cout << " ///////////// Two one-way vertex-edge snapping (A --> B) " << std::endl; #endif + visitor.start_first_vertex_edge_phase(); + snapped_n += internal::snap_non_conformal_one_way( halfedge_range_A, tm_A, tolerance_map_A, vertex_patch_map_A, locked_vertices_A, halfedge_range_B, tm_B, face_patch_map_B, locked_halfedges_B, - false /*source is never fixed*/, np_A, np_B); + false /*source is never fixed*/, visitor, np_A, np_B); + + visitor.end_first_vertex_edge_phase(); #ifdef CGAL_PMP_SNAP_DEBUG_OUTPUT std::ofstream("results/vertex_edge_A.off") << std::setprecision(17) << tm_A; @@ -1149,14 +1247,21 @@ std::size_t snap_non_conformal(HalfedgeRange& halfedge_range_A, if(!is_self_snapping) { + if(visitor.stop()) + return snapped_n; + #ifdef CGAL_PMP_SNAP_DEBUG - std::cout << " ///////////// Two one-way vertex-edge snapping (B --> A) " << std::endl; + std::cout << " ///////////// Two one-way vertex-edge snapping (B --> A) " << std::endl; #endif + visitor.start_second_vertex_edge_phase(); + snapped_n += internal::snap_non_conformal_one_way( halfedge_range_B, tm_B, tolerance_map_B, vertex_patch_map_B, locked_vertices_B, halfedge_range_A, tm_A, face_patch_map_A, locked_halfedges_A, - is_second_mesh_fixed, np_B, np_A); + is_second_mesh_fixed, visitor, np_B, np_A); + + visitor.end_second_vertex_edge_phase(); } return snapped_n; diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap_vertices.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap_vertices.h index 90df8a56406..802dc64aa48 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap_vertices.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap_vertices.h @@ -34,12 +34,18 @@ #include #include #include +#include +#include +#include #include +#include #include #include -#include #include +#include +#include #include +#include #include #include @@ -47,7 +53,9 @@ #include #ifdef CGAL_LINKED_WITH_TBB +#define TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS 1 # include +# include # include # include #endif @@ -88,26 +96,41 @@ struct Snapping_pair FT sq_dist; }; -// Functor that just forwards the pair of the two intersecting boxes -template -struct Intersecting_boxes_pairs_report +#ifdef CGAL_LINKED_WITH_TBB +// Functor that forwards the pair of the two intersecting boxes +// Note that Box_intersection_d does a lot of copies of the callback functor, but we are passing it +// with std::ref, so is always refering to the same reporter object (and importantly, the same counter) +template +struct Intersecting_boxes_pairs_parallel_report { - Intersecting_boxes_pairs_report(OutputIterator it) : m_iterator(it) { } + mutable OutputIterator m_iterator; + Visitor& m_visitor; + + Intersecting_boxes_pairs_parallel_report(OutputIterator it, Visitor& visitor) + : m_iterator(it), + m_visitor(visitor) + { } + + Intersecting_boxes_pairs_parallel_report(const Intersecting_boxes_pairs_parallel_report& other) = delete; + Intersecting_boxes_pairs_parallel_report& operator=(const Intersecting_boxes_pairs_parallel_report&) = delete; template - void operator()(const Box* b, const Box* c) const { - *m_iterator++ = std::make_pair(b->info(), c->info()); - } + void operator()(const Box* b, const Box* c) const + { + *m_iterator++ = std::make_pair(b->info(), c->info()); // writing into a concurrent vector - mutable OutputIterator m_iterator; + if(m_visitor.stop()) + throw CGAL::internal::Throw_at_output_exception(); + } }; +#endif template Unique_vertex_pair; +private: + const PolygonMesh& m_tm_A; + ToleranceMap_A m_tolerance_map_A; + VertexPatchMap_A m_vertex_patch_map_A; + VPM_A m_vpm_A; + const PolygonMesh& m_tm_B; + ToleranceMap_B m_tolerance_map_B; + VertexPatchMap_B m_vertex_patch_map_B; + VPM_B m_vpm_B; + const GeomTraits& m_gt; + + Visitor& m_visitor; + SnappingPairContainer& m_snapping_pairs; // will only be filled here in the sequential setting + +#ifdef CGAL_LINKED_WITH_TBB + const UniqueVertexPairContainer* m_uv_pairs; + ToKeepContainer* m_to_keep; +#endif + +public: Vertex_proximity_report(SnappingPairContainer& snapping_pairs, const PolygonMesh& tm_A, ToleranceMap_A tolerance_map_A, @@ -135,7 +178,8 @@ struct Vertex_proximity_report ToleranceMap_B tolerance_map_B, VertexPatchMap_B vertex_patch_map_B, const VPM_B& vpm_B, - const GeomTraits& gt + const GeomTraits& gt, + Visitor& visitor #ifdef CGAL_LINKED_WITH_TBB , const UniqueVertexPairContainer* uv_pairs = nullptr , ToKeepContainer* to_keep = nullptr @@ -151,6 +195,7 @@ struct Vertex_proximity_report m_vertex_patch_map_B(vertex_patch_map_B), m_vpm_B(vpm_B), m_gt(gt), + m_visitor(visitor), m_snapping_pairs(snapping_pairs) #ifdef CGAL_LINKED_WITH_TBB , m_uv_pairs(uv_pairs) @@ -158,7 +203,7 @@ struct Vertex_proximity_report #endif { } - bool are_equal_vertices(Unique_vertex_ptr uv_a, Unique_vertex_ptr uv_b) const + bool are_equal_vertices(const Unique_vertex_ptr uv_a, const Unique_vertex_ptr uv_b) const { if(&m_tm_A != &m_tm_B) // must be the same mesh return false; @@ -171,10 +216,13 @@ struct Vertex_proximity_report return (std::find(uv_b->first.begin(), uv_b->first.end(), ha) != uv_b->first.end()); } - std::pair do_keep(Unique_vertex_ptr uv_a, Unique_vertex_ptr uv_b) const + std::pair do_keep(const Unique_vertex_ptr uv_a, const Unique_vertex_ptr uv_b) const { + CGAL_precondition(uv_a != nullptr && uv_b != nullptr); + const Vertex_container& vs_a = uv_a->first; const Vertex_container& vs_b = uv_b->first; + const vertex_descriptor va = target(vs_a.front(), m_tm_A); const vertex_descriptor vb = target(vs_b.front(), m_tm_B); @@ -204,7 +252,9 @@ struct Vertex_proximity_report #endif if(res == CGAL::LARGER) + { return std::make_pair(-1, false); // ignore + } else { const FT sq_dist = (min)(m_gt.compute_squared_distance_3_object()(get(m_vpm_A, va), get(m_vpm_B, vb)), @@ -213,17 +263,29 @@ struct Vertex_proximity_report } } - // that's the sequential version - void operator()(const Box* a, const Box* b) + // sequential version (kd_tree) + void operator()(const Unique_vertex_ptr uv_a, const Unique_vertex_ptr uv_b) { - Unique_vertex_ptr uv_a = a->info(), uv_b = b->info(); + if(m_visitor.stop()) + throw CGAL::internal::Throw_at_output_exception(); + const std::pair res = do_keep(uv_a, uv_b); if(res.second) - m_snapping_pairs.insert(Snapping_pair(uv_a, uv_b, res.first)); + { + m_snapping_pairs.emplace(uv_a, uv_b, res.first); + m_visitor.on_vertex_vertex_match(uv_a->first, m_tm_A, uv_b->first, m_tm_B); + } } + // sequential version (box_d) + template + void operator()(const BoxPtr a, const BoxPtr b) + { + return this->operator()(a->info(), b->info()); + } + + // parallel version #ifdef CGAL_LINKED_WITH_TBB - // that's the parallel version void operator()(const tbb::blocked_range& r) const { CGAL_assertion(m_uv_pairs != nullptr); @@ -236,26 +298,384 @@ struct Vertex_proximity_report } } #endif - -private: - const PolygonMesh& m_tm_A; - ToleranceMap_A m_tolerance_map_A; - VertexPatchMap_A m_vertex_patch_map_A; - VPM_A m_vpm_A; - const PolygonMesh& m_tm_B; - ToleranceMap_B m_tolerance_map_B; - VertexPatchMap_B m_vertex_patch_map_B; - VPM_B m_vpm_B; - const GeomTraits& m_gt; - - SnappingPairContainer& m_snapping_pairs; // will only be filled here in the sequential setting - -#ifdef CGAL_LINKED_WITH_TBB - const UniqueVertexPairContainer* m_uv_pairs; - ToKeepContainer* m_to_keep; -#endif }; +template +struct Unique_point_position_accessor +{ + typedef Key key_type; + typedef ValueType value_type; + typedef const value_type& reference; + typedef boost::lvalue_property_map_tag category; + + friend reference get(const Unique_point_position_accessor&, const key_type& k) + { + return k->first.first; + } +}; + +template +void find_vertex_vertex_matches_with_kd_tree(Unique_positions& unique_positions_A, // not const because it complicates parallel loops + const Unique_positions& unique_positions_B, + const PolygonMesh& tm_A, + VPM_A vpm_A, + ToleranceMap_A tolerance_map_A, + VertexPatchMap_A vertex_patch_map_A, + const PolygonMesh& tm_B, + VPM_B vpm_B, + ToleranceMap_B tolerance_map_B, + VertexPatchMap_B vertex_patch_map_B, + const GT& gt, + Visitor& visitor, + Snapping_pair_container& snapping_pairs) +{ + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + typedef typename GT::FT FT; + typedef typename boost::property_traits::value_type Point; + + typedef std::vector Vertex_container; + typedef std::pair Unique_vertex; + typedef const Unique_vertex* Unique_vertex_ptr; + typedef std::pair Patch_point; + + // Put iterators into the tree: + // - to be lighter + // - because kd_tree.search() returns copies + typedef typename Unique_positions::const_iterator Key; + + typedef CGAL::Search_traits_3 Base_search_traits; + typedef CGAL::Search_traits_adapter, Base_search_traits> Search_traits; + typedef CGAL::Kd_tree Kd_tree; + typedef CGAL::Fuzzy_sphere Fuzzy_sphere; + + typedef std::vector Unique_positions_vector; + +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + CGAL::Real_timer timer; + timer.start(); +#endif + + Unique_positions_vector unique_positions_Bv; + unique_positions_Bv.reserve(unique_positions_B.size()); + for(auto it=unique_positions_B.cbegin(); it!=unique_positions_B.cend(); ++it) + unique_positions_Bv.push_back(it); + + Kd_tree tree; + tree.insert(unique_positions_Bv.begin(), unique_positions_Bv.end()); + +#if !defined(CGAL_LINKED_WITH_TBB) + CGAL_static_assertion_msg (!(std::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + // parallel + if(std::is_convertible::value) + { + typedef std::pair Unique_vertex_pair; + typedef tbb::concurrent_vector Unique_vertex_pairs; + + Unique_vertex_pairs uv_pairs; + + // Functor for the parallel_for over the conccurrent map + auto vertex_inquiry = [&](const typename Unique_positions::range_type &r) + { + for(const auto& e : r) + { + if(visitor.stop()) + return; + + const Patch_point& pp = e.first; + const Point& q = pp.first; + const Unique_vertex& uv = e.second; + const FT tolerance = uv.second; + + visitor.on_vertex_vertex_inquiry(uv, tm_A); + + // @fixme this assumes constant and equal tolerances, because querying the kd tree is one way + // and we don't know the tolerance at the potential match since we don't know the potential match... + Fuzzy_sphere fq(q, 2 * tolerance); + + std::vector res; + tree.search(std::back_inserter(res), fq); + + for(const auto rit : res) + { + uv_pairs.emplace_back(&uv, &(rit->second)); + CGAL_assertion(!rit->second.first.empty()); + } + } + }; + + tbb::parallel_for(unique_positions_A.range(), vertex_inquiry); + + // Filter the range of possible matches (also in parallel) + typedef std::vector > Filters; + typedef Vertex_proximity_report Reporter; + + Filters to_keep(uv_pairs.size()); + Reporter proximity_filterer(snapping_pairs, + tm_A, tolerance_map_A, vertex_patch_map_A, vpm_A, + tm_B, tolerance_map_B, vertex_patch_map_B, vpm_B, + gt, visitor, &uv_pairs, &to_keep); + + tbb::parallel_for(tbb::blocked_range(0, uv_pairs.size()), proximity_filterer); + + // Now fill the multi-index, sequentially + for(std::size_t i=0, uvps = uv_pairs.size(); ifirst, tm_A, uv_pairs[i].second->first, tm_B); + snapping_pairs.emplace(uv_pairs[i].first, uv_pairs[i].second, to_keep[i].first); + } + } + +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + timer.stop(); + std::cout << "time for KD Tree (parallel): " << timer.time() << std::endl; +#endif + } + else +#endif // CGAL_LINKED_WITH_TBB + // sequential + { + typedef std::pair Unique_vertex_pair; + typedef std::vector Unique_vertex_pairs; + + Unique_vertex_pairs uv_pairs; + + for(const auto& e : unique_positions_A) + { + if(visitor.stop()) + return; + + const Patch_point& pp = e.first; + const Point& q = pp.first; + const Unique_vertex& uv = e.second; + const FT tolerance = uv.second; + + visitor.on_vertex_vertex_inquiry(uv, tm_A); + + // @fixme this assumes constant and equal tolerances, because querying the kd tree is one way + // and we don't know the tolerance at the potential match since we don't know the potential match... + Fuzzy_sphere fq(q, 2 * tolerance); + + std::vector res; + tree.search(std::back_inserter(res), fq); + + for(const auto rit : res) + { + uv_pairs.emplace_back(&uv, &(rit->second)); + CGAL_assertion(!rit->second.first.empty()); + } + } + + typedef Vertex_proximity_report Reporter; + + Reporter proximity_filterer(snapping_pairs, + tm_A, tolerance_map_A, vertex_patch_map_A, vpm_A, + tm_B, tolerance_map_B, vertex_patch_map_B, vpm_B, + gt, visitor); + + for(const auto& p : uv_pairs) + proximity_filterer(p.first, p.second); + +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + timer.stop(); + std::cout << "time for KD Tree (sequential): " << timer.time() << std::endl; +#endif + } +} + +template +void find_vertex_vertex_matches_with_box_d(const Unique_positions& unique_positions_A, + const Unique_positions& unique_positions_B, + const PolygonMesh& tm_A, + VPM_A vpm_A, + ToleranceMap_A tolerance_map_A, + VertexPatchMap_A vertex_patch_map_A, + const PolygonMesh& tm_B, + VPM_B vpm_B, + ToleranceMap_B tolerance_map_B, + VertexPatchMap_B vertex_patch_map_B, + const GT& gt, + Visitor& visitor, + Snapping_pair_container& snapping_pairs) +{ + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + typedef typename GT::FT FT; + + typedef std::vector Vertex_container; + typedef std::pair Unique_vertex; + typedef const Unique_vertex* Unique_vertex_ptr; + + typedef Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy; + typedef CGAL::Box_intersection_d::Box_with_info_d Box; + +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + CGAL::Real_timer timer; + timer.start(); +#endif + + std::vector boxes_A; + boxes_A.reserve(unique_positions_A.size()); + std::vector boxes_B; + boxes_B.reserve(unique_positions_B.size()); + + // Actually build the boxes now + for(const auto& p : unique_positions_A) + { +#ifdef CGAL_PMP_SNAP_DEBUG_PP + std::cout << "Unique_vertex (A), pos: " << p.first.first << " vertices:"; + for(const halfedge_descriptor h : p.second.first) + std::cout << " " << target(h, tm_A); + std::cout << std::endl; +#endif + + const Unique_vertex& ev = p.second; + CGAL_assertion(!ev.first.empty()); + + // this only makes the box a little larger to ease intersection computations, + // the final tolerance is not changed + const double eps = 1.01 * CGAL::to_double(ev.second); + const Bbox_3 pb = gt.construct_bbox_3_object()(p.first.first); + const Bbox_3 b(pb.xmin() - eps, pb.ymin() - eps, pb.zmin() - eps, + pb.xmax() + eps, pb.ymax() + eps, pb.zmax() + eps); + boxes_A.emplace_back(b, &ev); + } + + for(const auto& p : unique_positions_B) + { +#ifdef CGAL_PMP_SNAP_DEBUG_PP + std::cout << "Unique_vertex (B), pos: " << p.first.first << " vertices:"; + for(const halfedge_descriptor h : p.second.first) + std::cout << " " << target(h, tm_B); + std::cout << std::endl; +#endif + + const Unique_vertex& ev = p.second; + CGAL_assertion(!ev.first.empty()); + + const double eps = 1.01 * CGAL::to_double(ev.second); + const Bbox_3 pb = gt.construct_bbox_3_object()(p.first.first); + const Bbox_3 b(pb.xmin() - eps, pb.ymin() - eps, pb.zmin() - eps, + pb.xmax() + eps, pb.ymax() + eps, pb.zmax() + eps); + boxes_B.emplace_back(b, &ev); + } + + // @fixme bench and don't use ptrs if not useful + std::vector boxes_A_ptr; + boxes_A_ptr.reserve(boxes_A.size()); + for(const Box& b : boxes_A) + boxes_A_ptr.push_back(&b); + + std::vector boxes_B_ptr; + boxes_B_ptr.reserve(boxes_B.size()); + for(const Box& b : boxes_B) + boxes_B_ptr.push_back(&b); + +#if !defined(CGAL_LINKED_WITH_TBB) + CGAL_static_assertion_msg (!(std::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else // CGAL_LINKED_WITH_TBB + if(std::is_convertible::value) + { + typedef std::pair Unique_vertex_pair; + typedef tbb::concurrent_vector Unique_vertex_pairs; + typedef std::back_insert_iterator UVP_output_iterator; + + Unique_vertex_pairs uv_pairs; + Intersecting_boxes_pairs_parallel_report box_callback(std::back_inserter(uv_pairs), visitor); + + // Shenanigans to pass a reference as callback (which is copied by value by 'box_intersection_d') + std::function callback(std::ref(box_callback)); + + // Grab the boxes that are interesecting but don't do any extra filtering (in parallel) + CGAL::box_intersection_d(boxes_A_ptr.begin(), boxes_A_ptr.end(), + boxes_B_ptr.begin(), boxes_B_ptr.end(), + callback); + + // Actually filter the range of intersecting boxes now (also in parallel) + typedef std::vector > Filters; + typedef Vertex_proximity_report Reporter; + + Filters to_keep(uv_pairs.size()); + Reporter proximity_filterer(snapping_pairs, + tm_A, tolerance_map_A, vertex_patch_map_A, vpm_A, + tm_B, tolerance_map_B, vertex_patch_map_B, vpm_B, + gt, visitor, &uv_pairs, &to_keep); + + tbb::parallel_for(tbb::blocked_range(0, uv_pairs.size()), proximity_filterer); + + // Now fill the multi-index, sequentially + for(std::size_t i=0, uvps = uv_pairs.size(); ifirst, uv_pairs[i].second->first); + snapping_pairs.emplace(uv_pairs[i].first, uv_pairs[i].second, to_keep[i].first); + } + } + +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + timer.stop(); + std::cout << "time for box_d (parallel): " << timer.time() << std::endl; +#endif + } + else +#endif // CGAL_LINKED_WITH_TBB + // Sequential code + { + typedef Vertex_proximity_report Reporter; + + Reporter proximity_filterer(snapping_pairs, + tm_A, tolerance_map_A, vertex_patch_map_A, vpm_A, + tm_B, tolerance_map_B, vertex_patch_map_B, vpm_B, gt, + visitor); + + // Shenanigans to pass a reference as callback (which is copied by value by 'box_intersection_d') + std::function callback(std::ref(proximity_filterer)); + + CGAL::box_intersection_d(boxes_A_ptr.begin(), boxes_A_ptr.end(), + boxes_B_ptr.begin(), boxes_B_ptr.end(), + callback); + +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + timer.stop(); + std::cout << "time for box_d (sequential): " << timer.time() << std::endl; +#endif + } +} + template ::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename GetVertexPointMap::type VPM_A; - typedef typename GetVertexPointMap::type VPM_B; + typedef typename GetVertexPointMap::type VPM_A; + typedef typename GetVertexPointMap::type VPM_B; - typedef typename GetGeomTraits::type GT; - typedef typename GT::FT FT; - typedef typename boost::property_traits::value_type Point; + typedef typename GetGeomTraits::type GT; + typedef typename GT::FT FT; + typedef typename boost::property_traits::value_type Point; - typedef std::vector Vertex_container; - typedef std::pair Unique_vertex; - typedef const Unique_vertex* Unique_vertex_ptr; - typedef std::map, Unique_vertex> Unique_positions; + typedef std::vector Vertex_container; + typedef std::pair Unique_vertex; + typedef const Unique_vertex* Unique_vertex_ptr; + typedef std::pair Patch_point; - typedef Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy; - typedef CGAL::Box_intersection_d::Box_with_info_d Box; + // @todo in theory could be unordered maps, but this requires defining the hashes & stuff +#ifdef CGAL_LINKED_WITH_TBB + typedef tbb::concurrent_map Unique_positions; +#else + typedef std::map Unique_positions; +#endif + + typedef typename internal_np::Lookup_named_param_def < + internal_np::visitor_t, + NamedParameters_A, + internal::Snapping_default_visitor // default + >::reference Visitor; - using parameters::get_parameter; using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::get_parameter_reference; CGAL_static_assertion((std::is_same::value)); @@ -306,6 +737,13 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, VPM_B vpm_B = choose_parameter(get_parameter(np_B, internal_np::vertex_point), get_property_map(vertex_point, tm_B)); + internal::Snapping_default_visitor default_visitor; + Visitor& visitor = choose_parameter(get_parameter_reference(np_A, internal_np::visitor), default_visitor); + + visitor.start_vertex_vertex_phase(); + if(visitor.stop()) + return 0; + #ifdef CGAL_PMP_SNAP_DEBUG std::cout << "Finding snappables vertices. Range sizes: " << std::distance(halfedge_range_A.begin(), halfedge_range_A.end()) << " and " @@ -316,15 +754,10 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, is_empty_range(halfedge_range_B.begin(), halfedge_range_B.end())) return 0; - // Vertex-Vertex snapping is performed as follows: - // - Identify points which are already equal and group them together so that they are moved together - // - Create a single box for these points - - std::vector boxes_A; - boxes_A.reserve(halfedge_range_A.size()); - std::vector boxes_B; - boxes_B.reserve(halfedge_range_B.size()); + // Vertex-Vertex snapping is performed by identify points which are already equal + // and grouping them together so that they are moved together + // @todo all range building could be parallel (worth it?) Unique_positions unique_positions_A; for(halfedge_descriptor h : halfedge_range_A) { @@ -336,10 +769,9 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, #endif Vertex_container nvc {{ h }}; - std::pair is_insert_successful = - unique_positions_A.insert(std::make_pair( - std::make_pair(get(vpm_A, v), get(vertex_patch_map_A, v)), // point and patch id - std::make_pair(nvc, tolerance))); + auto is_insert_successful = unique_positions_A.emplace(std::make_pair(get(vpm_A, v), + get(vertex_patch_map_A, v)), + std::make_pair(nvc, tolerance)); if(!is_insert_successful.second) // point was already met { @@ -351,7 +783,7 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, } } - // same for tm_B (@todo avoid all that for self snapping + use self_intersection_d) + // same for tm_B (@todo when doing boxes, avoid all that for self snapping + use self_intersection_d) Unique_positions unique_positions_B; for(halfedge_descriptor h : halfedge_range_B) { @@ -364,9 +796,8 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, Vertex_container nvc {{ h }}; std::pair is_insert_successful = - unique_positions_B.insert(std::make_pair( - std::make_pair(get(vpm_B, v), get(vertex_patch_map_B, v)), // point and patch id - std::make_pair(nvc, tolerance))); + unique_positions_B.emplace(std::make_pair(get(vpm_B, v), get(vertex_patch_map_B, v)), // point and patch id + std::make_pair(nvc, tolerance)); if(!is_insert_successful.second) // point was already met { @@ -378,62 +809,10 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, } } - // Actually build the boxes now - for(const auto& p : unique_positions_A) - { -#ifdef CGAL_PMP_SNAP_DEBUG_PP - std::cout << "Unique_vertex (A), pos: " << p.first.first << " vertices:"; - for(const halfedge_descriptor h : p.second.first) - std::cout << " " << target(h, tm_A); - std::cout << std::endl; -#endif - - const Unique_vertex& ev = p.second; - CGAL_assertion(!ev.first.empty()); - - // this only makes the box a little larger to ease intersection computations, - // the final tolerance is not changed - const double eps = 1.01 * CGAL::to_double(ev.second); - const Bbox_3 pb = gt.construct_bbox_3_object()(p.first.first); - const Bbox_3 b(pb.xmin() - eps, pb.ymin() - eps, pb.zmin() - eps, - pb.xmax() + eps, pb.ymax() + eps, pb.zmax() + eps); - boxes_A.push_back(Box(b, &ev)); - } - - for(const auto& p : unique_positions_B) - { -#ifdef CGAL_PMP_SNAP_DEBUG_PP - std::cout << "Unique_vertex (B), pos: " << p.first.first << " vertices:"; - for(const halfedge_descriptor h : p.second.first) - std::cout << " " << target(h, tm_B); - std::cout << std::endl; -#endif - - const Unique_vertex& ev = p.second; - CGAL_assertion(!ev.first.empty()); - - const double eps = 1.01 * CGAL::to_double(ev.second); - const Bbox_3 pb = gt.construct_bbox_3_object()(p.first.first); - const Bbox_3 b(pb.xmin() - eps, pb.ymin() - eps, pb.zmin() - eps, - pb.xmax() + eps, pb.ymax() + eps, pb.zmax() + eps); - boxes_B.push_back(Box(b, &ev)); - } - - // @fixme bench and don't use ptrs if not useful - std::vector boxes_A_ptr; - boxes_A_ptr.reserve(boxes_A.size()); - for(const Box& b : boxes_A) - boxes_A_ptr.push_back(&b); - - std::vector boxes_B_ptr; - boxes_B_ptr.reserve(boxes_B.size()); - for(const Box& b : boxes_B) - boxes_B_ptr.push_back(&b); - // Use a multi index to sort easily by sources, targets, AND distance. // Then, look up the distances in increasing order, and snap whenever the source and the target // have both not been snapped yet. - typedef internal::Snapping_pair Snapping_pair; + typedef internal::Snapping_pair Snapping_pair; typedef boost::multi_index::multi_index_container< Snapping_pair, boost::multi_index::indexed_by< @@ -444,77 +823,55 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, boost::multi_index::ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(Snapping_pair, FT, sq_dist)> > - > Snapping_pair_container; + > Snapping_pair_container; Snapping_pair_container snapping_pairs; - -#if !defined(CGAL_LINKED_WITH_TBB) - CGAL_static_assertion_msg (!(std::is_convertible::value), - "Parallel_tag is enabled but TBB is unavailable."); -#else - if(std::is_convertible::value) + try { - typedef std::pair Unique_vertex_pair; - typedef tbb::concurrent_vector Unique_vertex_pairs; - typedef std::back_insert_iterator UVP_output_iterator; - - Unique_vertex_pairs uv_pairs; - Intersecting_boxes_pairs_report callback(std::back_inserter(uv_pairs)); - - CGAL::Real_timer timer; - timer.start(); - - // Grab the boxes that are interesecting but don't do any extra filtering (in parallel) - CGAL::box_intersection_d(boxes_A_ptr.begin(), boxes_A_ptr.end(), - boxes_B_ptr.begin(), boxes_B_ptr.end(), - callback); - - std::cout << "time for box_d: " << timer.time() << std::endl; - - // Actually filter the range of intersecting boxes now (in parallel) - typedef std::vector > Filters; - typedef Vertex_proximity_report Reporter; - - Filters to_keep(uv_pairs.size()); - Reporter proximity_filterer(snapping_pairs, tm_A, tolerance_map_A, vertex_patch_map_A, vpm_A, - tm_B, tolerance_map_B, vertex_patch_map_B, vpm_B, - gt, &uv_pairs, &to_keep); - tbb::parallel_for(tbb::blocked_range(0, uv_pairs.size()), proximity_filterer); - - // Now fill the multi index, sequentially - for(std::size_t i=0, uvps = uv_pairs.size(); i(unique_positions_A, unique_positions_B, + tm_A, vpm_A, tolerance_map_A, vertex_patch_map_A, + tm_B, vpm_B, tolerance_map_B, vertex_patch_map_B, gt, + visitor, snapping_pairs); +#else // CGAL_PMP_SNAP_VERTICES_USE_KD_TREE + find_vertex_vertex_matches_with_box_d(unique_positions_A, unique_positions_B, + tm_A, vpm_A, tolerance_map_A, vertex_patch_map_A, + tm_B, vpm_B, tolerance_map_B, vertex_patch_map_B, gt, + visitor, snapping_pairs); +#endif // CGAL_PMP_SNAP_VERTICES_USE_KD_TREE } - else -#endif + catch(const CGAL::internal::Throw_at_output_exception&) { - typedef Vertex_proximity_report Reporter; - - Reporter vpr(snapping_pairs, tm_A, tolerance_map_A, vertex_patch_map_A, vpm_A, - tm_B, tolerance_map_B, vertex_patch_map_B, vpm_B, gt); - - // Shenanigans to pass a reference as callback (which is copied by value by 'box_intersection_d') - std::function callback(std::ref(vpr)); - - CGAL::box_intersection_d(boxes_A_ptr.begin(), boxes_A_ptr.end(), - boxes_B_ptr.begin(), boxes_B_ptr.end(), - callback); + return 0; } +#if TBB_USE_CAPTURED_EXCEPTION + catch(const tbb::captured_exception& e) + { + const std::string tn1(e.name()); + const std::string tn2(typeid(const CGAL::internal::Throw_at_output_exception&).name()); + if(tn1 != tn2) + { + std::cerr << "Unexpected throw: " << tn1 << std::endl; + throw; + } + + return 0; + } +#endif // TBB_USE_CAPTURED_EXCEPTION ////////////////////////////////////////////////////////////////////////////////////////////////// /// Done collecting; start matching ////////////////////////////////////////////////////////////////////////////////////////////////// - if(snapping_pairs.empty()) + if(snapping_pairs.empty() || visitor.stop()) return 0; +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + CGAL::Real_timer timer; + timer.start(); +#endif + typedef std::pair Unique_vertex_pair; std::vector snappable_vertices_pairs; @@ -566,15 +923,21 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, /// Done matching; start snapping ////////////////////////////////////////////////////////////////////////////////////////////////// - std::size_t counter = 0; + visitor.start_vertex_vertex_snapping(); + if(snappable_vertices_pairs.empty() || visitor.stop()) + return 0; #ifdef CGAL_PMP_SNAP_DEBUG_OUTPUT std::ofstream out_edges("results/snappable.polylines.txt"); out_edges.precision(17); #endif + std::size_t counter = 0; for(const Unique_vertex_pair& uvp : snappable_vertices_pairs) { + if(visitor.stop()) + return counter; + Unique_vertex_ptr uv_a = uvp.first; Unique_vertex_ptr uv_b = uvp.second; const Vertex_container& vs_a = uv_a->first; @@ -582,6 +945,8 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, const vertex_descriptor va = target(vs_a.front(), tm_A); const vertex_descriptor vb = target(vs_b.front(), tm_B); + visitor.before_vertex_vertex_snap(vs_a, tm_A, vs_b, tm_B); + #ifdef CGAL_PMP_SNAP_DEBUG_OUTPUT out_edges << "2 " << tm_A.point(va) << " " << tm_B.point(vb) << std::endl; #endif @@ -615,12 +980,16 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, ++counter; } + + visitor.after_vertex_vertex_snap(vs_a, tm_A, vs_b, tm_B); } #ifdef CGAL_PMP_SNAP_DEBUG_OUTPUT out_edges.close(); #endif + visitor.end_vertex_vertex_snapping(); + #ifdef CGAL_PMP_SNAP_DEBUG std::cout << "Snapped " << counter << " pair(s)!" << std::endl; #endif @@ -629,7 +998,7 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, /// Done snapping; start analyzing ////////////////////////////////////////////////////////////////////////////////////////////////// - // Below is used in non-conformal snapping + // Below is performed to improve non-conformal snapping // @todo could avoid doing it if not required // // Now that vertex-vertex snapping has been performed, look around to see if we can already @@ -640,8 +1009,7 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, // #2 : If pairs on either side of the two matching vertices are compatible (not necessary fully matching), // then the two vertices should be locked as no better match can be obtained // #3 : If a pair of incident edges are not fully matching, but still have compatible directions - // (i.e. collinear and opposite directions), then we don't want to project onto the shorter - // of the two + // (i.e. collinear and opposite directions), then we don't want to project onto the shorter of the two for(const Unique_vertex_pair& uvp : snappable_vertices_pairs) { Unique_vertex_ptr uv_a = uvp.first; @@ -724,6 +1092,13 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, } } + visitor.end_vertex_vertex_phase(); + +#ifdef CGAL_PMP_SNAPPING_PRINT_RUNTIME + timer.stop(); + std::cout << "time for the actual snapping (vertex-vertex): " << timer.time() << std::endl; +#endif + return counter; } @@ -750,15 +1125,15 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A, using CGAL::parameters::choose_parameter; using CGAL::parameters::get_parameter; - return snap_vertices_two_way(halfedge_range_A, tm_A, tolerance_map_A, - choose_parameter(get_parameter(np_A, internal_np::vertex_incident_patches), - Constant_property_map(-1)), - halfedge_range_B, tm_B, tolerance_map_B, - choose_parameter(get_parameter(np_B, internal_np::vertex_incident_patches), - Constant_property_map(-1)), - lockable_vps_out, lockable_ha_out, lockable_hb_out, - choose_parameter(get_parameter(np_B, internal_np::do_lock_mesh), false), - np_A, np_B); + return snap_vertices_two_way(halfedge_range_A, tm_A, tolerance_map_A, + choose_parameter(get_parameter(np_A, internal_np::vertex_incident_patches), + Constant_property_map(-1)), + halfedge_range_B, tm_B, tolerance_map_B, + choose_parameter(get_parameter(np_B, internal_np::vertex_incident_patches), + Constant_property_map(-1)), + lockable_vps_out, lockable_ha_out, lockable_hb_out, + choose_parameter(get_parameter(np_B, internal_np::do_lock_mesh), false), + np_A, np_B); } } // namespace internal @@ -818,7 +1193,8 @@ namespace experimental { // // \return the number of snapped vertex pairs // -template std::size_t snap_vertices(const HalfedgeRange_A& halfedge_range_A, @@ -832,13 +1208,14 @@ std::size_t snap_vertices(const HalfedgeRange_A& halfedge_range_A, { CGAL::Emptyset_iterator unused_output_iterator; - return internal::snap_vertices_two_way(halfedge_range_A, tm_A, tolerance_map_A, - halfedge_range_B, tm_B, tolerance_map_B, - unused_output_iterator, unused_output_iterator, - unused_output_iterator, np_A, np_B); + return internal::snap_vertices_two_way(halfedge_range_A, tm_A, tolerance_map_A, + halfedge_range_B, tm_B, tolerance_map_B, + unused_output_iterator, unused_output_iterator, + unused_output_iterator, np_A, np_B); } -template std::size_t snap_vertices(const HalfedgeRange_A& halfedge_range_A, PolygonMesh& tm_A, @@ -847,11 +1224,13 @@ std::size_t snap_vertices(const HalfedgeRange_A& halfedge_range_A, PolygonMesh& tm_B, ToleranceMap_B tolerance_map_B) { - return snap_vertices(halfedge_range_A, tm_A, tolerance_map_A, halfedge_range_B, tm_B, tolerance_map_B, - CGAL::parameters::all_default(), CGAL::parameters::all_default()); + return snap_vertices(halfedge_range_A, tm_A, tolerance_map_A, + halfedge_range_B, tm_B, tolerance_map_B, + CGAL::parameters::all_default(), CGAL::parameters::all_default()); } -template std::size_t snap_vertices(const HalfedgeRange_A& halfedge_range_A, @@ -877,21 +1256,23 @@ std::size_t snap_vertices(const HalfedgeRange_A& halfedge_range_A, return snap_vertices(halfedge_range_A, tm_A, tolerance_map_A, halfedge_range_B, tm_B, tolerance_map_B, np_A, np_B); } -template +template std::size_t snap_vertices(const HalfedgeRange_A& halfedge_range_A, PolygonMesh& tm_A, const HalfedgeRange_B& halfedge_range_B, PolygonMesh& tm_B) { - return snap_vertices(halfedge_range_A, tm_A, halfedge_range_B, tm_B, - parameters::all_default(), parameters::all_default()); + return snap_vertices(halfedge_range_A, tm_A, halfedge_range_B, tm_B, + parameters::all_default(), parameters::all_default()); } /////////////////////////////////////////////////////////////////////////////////////////////////// /// Border convenience overloads /////////////////////////////////////////////////////////////////////////////////////////////////// -template +template std::size_t snap_border_vertices(PolygonMesh& tm_A, ToleranceMap_A tolerance_map_A, PolygonMesh& tm_B, @@ -904,11 +1285,12 @@ std::size_t snap_border_vertices(PolygonMesh& tm_A, std::vector border_B; border_halfedges(tm_B, std::back_inserter(border_B)); - return snap_vertices(border_A, tm_A, tolerance_map_A, - border_B, tm_B, tolerance_map_B); + return snap_vertices(border_A, tm_A, tolerance_map_A, + border_B, tm_B, tolerance_map_B); } -template +template std::size_t snap_border_vertices(PolygonMesh& tm_A, PolygonMesh& tm_B) { @@ -920,25 +1302,28 @@ std::size_t snap_border_vertices(PolygonMesh& tm_A, std::vector border_vertices_B; border_halfedges(tm_B, std::back_inserter(border_vertices_B)); - return snap_vertices(border_vertices_A, tm_A, border_vertices_B, tm_B); + return snap_vertices(border_vertices_A, tm_A, border_vertices_B, tm_B); } -template +template std::size_t snap_border_vertices(PolygonMesh& tm, ToleranceMap tolerance_map) { - return snap_border_vertices(tm, tolerance_map, tm, tolerance_map); + return snap_border_vertices(tm, tolerance_map, tm, tolerance_map); } -template +template std::size_t snap_border_vertices(PolygonMesh& tm) { - return snap_border_vertices(tm, tm); + return snap_border_vertices(tm, tm); } /////////////////////////////////////////////////////////////////////////////////////////////////// /// Other convenience overloads /////////////////////////////////////////////////////////////////////////////////////////////////// -template +template std::size_t snap_all_vertices(PolygonMesh& tm, ToleranceMap tolerance_map) { typedef boost::graph_traits GT; @@ -952,10 +1337,9 @@ std::size_t snap_all_vertices(PolygonMesh& tm, ToleranceMap tolerance_map) boost::make_transform_iterator(vertices(tm).begin(), get_halfedge), boost::make_transform_iterator(vertices(tm).end(), get_halfedge) ); - return snap_vertices(hedges, tm, tolerance_map, hedges, tm, tolerance_map); + return snap_vertices(hedges, tm, tolerance_map, hedges, tm, tolerance_map); } - } // namespace experimental } // namespace Polygon_mesh_processing } // namespace CGAL From b92c22d0a9acb8903e90c6e3a12c00a612b1a086 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 28 Sep 2021 16:49:19 +0200 Subject: [PATCH 2/2] Improve preprocessing of snapping (don't move fixed vertices during collapsing) --- .../internal/Snapping/snap.h | 27 ++++++++++++++----- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h index 5f5e807d49e..490cc457b13 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h @@ -92,6 +92,7 @@ void simplify_range(HalfedgeRange& halfedge_range, typedef typename GT::FT FT; typedef typename GetVertexPointMap::type VPM; + typedef typename boost::property_traits::value_type Point; typedef typename boost::property_traits::reference Point_ref; using parameters::get_parameter; @@ -125,10 +126,10 @@ void simplify_range(HalfedgeRange& halfedge_range, // @fixme what if the source vertex is not to be snapped? Tolerance cannot be obtained... // and where should the post-collapse vertex be since we can't move the source vertex... // --> simply don't collapse? - const FT min_tol = (std::min)(get(tolerance_map, vs), get(tolerance_map, vt)); - const FT max_tol = (std::max)(get(tolerance_map, vs), get(tolerance_map, vt)); + const FT tol_s = get(tolerance_map, vs), tol_t = get(tolerance_map, vt); + const FT max_tol = (std::max)(tol_s, tol_t); - if(gt.compare_squared_distance_3_object()(ps,pt,CGAL::square(max_tol))==SMALLER) + if(gt.compare_squared_distance_3_object()(ps,pt,CGAL::square(max_tol)) == SMALLER) { const halfedge_descriptor prev_h = prev(h, tm); const halfedge_descriptor next_h = next(h, tm); @@ -136,11 +137,23 @@ void simplify_range(HalfedgeRange& halfedge_range, // check that the border has at least 4 edges not to create degenerate volumes if(border_size(h, tm) >= 4) { - const FT h_sq_length = gt.compute_squared_distance_3_object()(ps, pt); - vertex_descriptor v = Euler::collapse_edge(edge(h, tm), tm); + const FT lambda = tol_s / (tol_s + tol_t); + const Point new_p = ps + lambda * (pt - ps); - put(vpm, v, gt.construct_midpoint_3_object()(ps, pt)); - put(tolerance_map, v, min_tol + 0.5 * CGAL::approximate_sqrt(h_sq_length)); + // If we're collapsing onto a vertex that is not allowed to move, keep it fixed + const FT min_tol = (std::min)(tol_s, tol_t); + FT new_tolerance = min_tol; + if(!is_zero(min_tol)) + { + if(tol_t > tol_s) // the new point is closer to ps than to pt + new_tolerance += CGAL::approximate_sqrt(CGAL::squared_distance(new_p, ps)); + else + new_tolerance += CGAL::approximate_sqrt(CGAL::squared_distance(new_p, pt)); + } + + vertex_descriptor v = Euler::collapse_edge(edge(h, tm), tm); + put(vpm, v, new_p); + put(tolerance_map, v, new_tolerance); if(get(range_halfedges, prev_h)) edges_to_test.insert(prev_h);