diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h index 6b6478f2a65..94e43149daf 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -74,43 +75,47 @@ class CollapseTriangulation typedef std::array Facet; typedef std::array Tet; - std::unordered_set vertices_to_insert; - for (Cell_handle ch : cells_to_insert) + /*vertex of main tr - vertex of collapse tr*/ + // the star holds around twenty distinct vertices, so a linear scan over + // handles compares cheaper than hashing them : hashing a compact-container + // iterator goes through Time_stamper, and this runs on every attempt + boost::container::small_vector, 32> v2i; + const auto index_of = [&v2i](const Vertex_handle vh) -> int { - for(int i = 0; i < 4; ++i) - vertices_to_insert.insert(ch->vertex(i)); - } - - CGAL_expensive_assertion(vertices_to_insert.end() - != std::find(vertices_to_insert.begin(), vertices_to_insert.end(), v1_init)); - CGAL_expensive_assertion(vertices_to_insert.end() - != std::find(vertices_to_insert.begin(), vertices_to_insert.end(), v0_init)); - - std::unordered_map v2i;/*vertex of main tr - vertex of collapse tr*/ + for (const std::pair& p : v2i) + if (p.first == vh) + return p.second; + return -1; + }; - //To add the vertices only once std::vector points; - int index = 0; - for (Vertex_handle vh : vertices_to_insert) - { - if (v2i.find(vh) == v2i.end()) - { - points.push_back(vh->point()); - v2i.insert(std::make_pair(vh, index++)); - } - } - std::vector finite_cells; std::vector subdomains; + finite_cells.reserve(cells_to_insert.size()); + subdomains.reserve(cells_to_insert.size()); + for (Cell_handle ch : cells_to_insert) { - finite_cells.push_back( { v2i.at(ch->vertex(0)), - v2i.at(ch->vertex(1)), - v2i.at(ch->vertex(2)), - v2i.at(ch->vertex(3)) } ); + Tet tet; + for (int i = 0; i < 4; ++i) + { + const Vertex_handle vh = ch->vertex(i); + int id = index_of(vh); + if (id == -1) + { + id = static_cast(points.size()); + v2i.emplace_back(vh, id); + points.push_back(vh->point()); + } + tet[i] = id; + } + finite_cells.push_back(tet); subdomains.push_back(ch->subdomain_index()); } + CGAL_expensive_assertion(index_of(v1_init) != -1); + CGAL_expensive_assertion(index_of(v0_init) != -1); + // finished std::vector new_vertices; std::map border_facets; @@ -124,8 +129,8 @@ class CollapseTriangulation CGAL_assertion(triangulation.infinite_vertex() == new_vertices[0]); // update() - vh0 = new_vertices[v2i.at(v0_init) + 1]; - vh1 = new_vertices[v2i.at(v1_init) + 1]; + vh0 = new_vertices[index_of(v0_init) + 1]; + vh1 = new_vertices[index_of(v1_init) + 1]; Cell_handle ch; int i0, i1; @@ -1039,36 +1044,33 @@ bool is_cells_set_manifold(const C3t3&, typedef std::array FV; typedef std::pair EV; - std::unordered_map> facets; + // A facet is shared by exactly two cells, so it bounds the set when its + // neighbour is outside : the triangulation already answers that, and asking + // it costs one lookup of a cell handle where counting the facets of the set + // meant hashing a triple of vertex handles for every facet of every cell. + std::unordered_map> edges; + edges.reserve(4 * cells.size()); + for (Cell_handle c : cells) { for (int i = 0; i < 4; ++i) { + if (cells.find(c->neighbor(i)) != cells.end()) + continue; // shared with another cell of the set + const FV fvi = make_vertex_array(c->vertex((i + 1) % 4), c->vertex((i + 2) % 4), c->vertex((i + 3) % 4)); - typename std::unordered_map>::iterator fit = facets.find(fvi); - if (fit == facets.end()) - facets.insert(std::make_pair(fvi, 1)); - else - fit->second++; - } - } - std::unordered_map> edges; - for (const auto& fvv : facets) - { - if (fvv.second != 1) - continue; - - for (int i = 0; i < 3; ++i) - { - const EV evi = make_vertex_pair(fvv.first[i], fvv.first[(i + 1) % 3]); - typename std::unordered_map>::iterator eit = edges.find(evi); - if (eit == edges.end()) - edges.insert(std::make_pair(evi, 1)); - else - eit->second++; + for (int k = 0; k < 3; ++k) + { + const EV evi = make_vertex_pair(fvi[k], fvi[(k + 1) % 3]); + typename std::unordered_map>::iterator eit = edges.find(evi); + if (eit == edges.end()) + edges.insert(std::make_pair(evi, 1)); + else + eit->second++; + } } } @@ -1252,8 +1254,12 @@ void collapse_short_edges(C3T3& c3t3, typedef typename T3::Vertex_handle Vertex_handle; typedef typename T3::Geom_traits::FT FT; + // the element side is only ever searched, never walked in order, so it is + // hashed : keeping it sorted meant comparing pairs of vertex handles + // O(log n) times for every edge the last collapse touched. The priority + // side keeps its order, which is what decides what runs next. typedef boost::bimap< - boost::bimaps::set_of >, + boost::bimaps::unordered_set_of, Equal_edges >, boost::bimaps::multiset_of > > Boost_bimap; typedef typename Boost_bimap::value_type short_edge; diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h index 15f8a868ca3..080bb511aae 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h @@ -29,8 +29,10 @@ #include #include +#include #include #include +#include #include #include @@ -615,6 +617,26 @@ struct Compare_edges } }; +// Same equivalence as `Compare_edges`, for the indices that only need to find +// an edge again rather than keep the edges in order. +template +struct Hash_edges +{ + std::size_t operator()(const Edge& e) const + { + return boost::hash()(make_vertex_pair(e)); + } +}; + +template +struct Equal_edges +{ + bool operator()(const Edge& e1, const Edge& e2) const + { + return make_vertex_pair(e1) == make_vertex_pair(e2); + } +}; + template CGAL::Triple make_vertex_triple(const Vh vh0, const Vh vh1, const Vh vh2) @@ -964,16 +986,33 @@ OutputIterator incident_surface_patches(const typename C3t3::Vertex_handle& v, return oit; } +// Counting how many distinct indices a simplex touches needs no container of +// its own : a vertex or an edge sees a handful of them, so the ones already +// seen are kept inline and scanned. +template +struct Distinct_index_counter +{ + boost::container::small_vector seen; + + void operator()(const Index& i) + { + for (const Index& s : seen) + if (s == i) + return; + seen.push_back(i); + } +}; + template std::size_t nb_incident_subdomains(const typename C3t3::Vertex_handle v, const C3t3& c3t3) { typedef typename C3t3::Subdomain_index Subdomain_index; - std::unordered_set indices; - incident_subdomains(v, c3t3, std::inserter(indices, indices.begin())); + Distinct_index_counter counter; + incident_subdomains(v, c3t3, boost::make_function_output_iterator(std::ref(counter))); - return indices.size(); + return counter.seen.size(); } template @@ -982,10 +1021,10 @@ std::size_t nb_incident_subdomains(const typename C3t3::Edge& e, { typedef typename C3t3::Subdomain_index Subdomain_index; - std::unordered_set indices; - incident_subdomains(e, c3t3, std::inserter(indices, indices.begin())); + Distinct_index_counter counter; + incident_subdomains(e, c3t3, boost::make_function_output_iterator(std::ref(counter))); - return indices.size(); + return counter.seen.size(); } template @@ -994,10 +1033,10 @@ std::size_t nb_incident_surface_patches(const typename C3t3::Edge& e, { typedef typename C3t3::Surface_patch_index Surface_patch_index; - std::unordered_set> indices; - incident_surface_patches(e, c3t3, std::inserter(indices, indices.begin())); + Distinct_index_counter counter; + incident_surface_patches(e, c3t3, boost::make_function_output_iterator(std::ref(counter))); - return indices.size(); + return counter.seen.size(); } template