diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/Elementary_operation.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/Elementary_operation.h new file mode 100644 index 00000000000..14aa854b178 --- /dev/null +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/Elementary_operation.h @@ -0,0 +1,98 @@ +// Copyright (c) 2025 GeometryFactory (France) and Telecom Paris (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Iasonas Manolas, Jane Tournois + +#ifndef CGAL_TETRAHEDRAL_REMESHING_ELEMENTARY_OPERATIONS_H +#define CGAL_TETRAHEDRAL_REMESHING_ELEMENTARY_OPERATIONS_H + +#include + +#include + +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE +#include +#include +#include +#endif + +namespace CGAL { +namespace Tetrahedral_remeshing { +namespace internal { + +template +class Elementary_operation +{ +public: + using C3t3 = C3t3_; + using Triangulation = typename C3t3::Triangulation; + using Element_type = ElementType; + using Element_range = ElementRange; + + Elementary_operation() = default; + virtual ~Elementary_operation() = default; + + virtual Element_range get_elements(const C3t3& c3t3) const = 0; + virtual bool execute_operation(const Element_type& e, C3t3& c3t3) = 0; + virtual std::string operation_name() const = 0; +}; + +template +class Elementary_operation_execution_sequential +{ +public: + using C3t3 = typename Operation::C3t3; + using Element_range = typename Operation::Element_range; + + bool execute(Operation& op, C3t3& c3t3) const + { + Element_range candidates = op.get_elements(c3t3); + if (candidates.empty()) + return false; + +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + std::size_t nb_done = 0; + const std::size_t nb_candidates = candidates.size(); + CGAL::Real_timer timer; + timer.start(); +#endif +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS + std::size_t nb_processed = 0; +#endif + for (const auto& element : candidates) + { + if (op.execute_operation(element, c3t3)) + { +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + ++nb_done; +#endif + } +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS + std::cout << "\r" << op.operation_name() << "... (" + << ++nb_processed << "/" << nb_candidates << ")"; + std::cout.flush(); +#endif + } + +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + timer.stop(); + std::cout << op.operation_name() << ": " << nb_done << "/" + << nb_candidates << " done (" + << timer.time() << " sec)." << std::endl; +#endif + return true; + } +}; + +} // namespace internal +} // namespace Tetrahedral_remeshing +} // namespace CGAL + +#endif // CGAL_TETRAHEDRAL_REMESHING_ELEMENTARY_OPERATIONS_H 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..c7cf9f1edd7 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 @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -31,6 +32,7 @@ #include #include +#include #include #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE @@ -977,7 +979,7 @@ collapse(const typename C3t3::Cell_handle ch, template -typename C3t3::Vertex_handle collapse(typename C3t3::Edge& edge, +typename C3t3::Vertex_handle collapse(const typename C3t3::Edge& edge, const Collapse_type& collapse_type, CellSelector& cell_selector, C3t3& c3t3, @@ -1084,7 +1086,7 @@ template -typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge, +typename C3t3::Vertex_handle collapse_edge(const typename C3t3::Edge& edge, C3t3& c3t3, const Sizing& sizing, const bool /* protect_boundaries */, @@ -1237,125 +1239,163 @@ auto can_be_collapsed(const typename C3T3::Edge& e, return Collapsible {true, boundary}; } -template +using Short_edges_bimap = boost::bimap< + boost::bimaps::set_of >, + boost::bimaps::multiset_of > >; + +template -void collapse_short_edges(C3T3& c3t3, - const Sizing& sizing, - const bool protect_boundaries, - CellSelector cell_selector, - Visitor& visitor) +class Edge_collapse_operation + : public Elementary_operation > { - typedef typename C3T3::Triangulation T3; - typedef typename T3::Edge Edge; - typedef typename T3::Vertex_handle Vertex_handle; - - typedef typename T3::Geom_traits::FT FT; - typedef boost::bimap< - boost::bimaps::set_of >, - boost::bimaps::multiset_of > > Boost_bimap; - typedef typename Boost_bimap::value_type short_edge; - - T3& tr = c3t3.triangulation(); - -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - std::cout << "Collapse short edges..."; - std::cout.flush(); - std::size_t nb_collapses = 0; - CGAL::Real_timer timer; - timer.start(); -#endif +public: + using Tr = typename C3t3::Triangulation; + using Vertex_handle = typename Tr::Vertex_handle; + using Edge = typename Tr::Edge; + using FT = typename Tr::Geom_traits::FT; + + using Short_edges = Short_edges_bimap; + using Base_operation = Elementary_operation; + using Element_type = typename Base_operation::Element_type; + using Element_range = typename Base_operation::Element_range; + +private: + const SizingFunction& m_sizing; + const CellSelector& m_cell_selector; + bool m_protect_boundaries; + Visitor& m_visitor; - //collect long edges - Boost_bimap short_edges; - for (const Edge& e : tr.finite_edges()) +public: + Edge_collapse_operation(const SizingFunction& sizing, + const CellSelector& cell_selector, + const bool protect_boundaries, + Visitor& visitor) + : m_sizing(sizing) + , m_cell_selector(cell_selector) + , m_protect_boundaries(protect_boundaries) + , m_visitor(visitor) {} + + Element_range get_elements(const C3t3& c3t3) const override { - auto [collapsible, boundary] = can_be_collapsed(e, c3t3, protect_boundaries, cell_selector); - if (!collapsible) - continue; - - const auto sqlen = is_too_short(e, boundary, sizing, c3t3, cell_selector); - if(sqlen != std::nullopt) - short_edges.insert(short_edge(e, sqlen.value())); + Short_edges short_edges; + for (const Edge& e : c3t3.triangulation().finite_edges()) + { + auto [collapsible, boundary] + = can_be_collapsed(e, c3t3, m_protect_boundaries, m_cell_selector); + if (!collapsible) + continue; + + const auto sqlen = is_too_short(e, boundary, m_sizing, c3t3, m_cell_selector); + if (sqlen != std::nullopt) + short_edges.insert(typename Short_edges::value_type(e, sqlen.value())); + } + return short_edges; } -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - debug::dump_edges(short_edges, "short_edges.polylines.txt"); - - std::ofstream short_success("short_collapse_success.polylines.txt"); - std::ofstream short_fail("short_collapse_fail.polylines.txt"); - std::ofstream short_cancel("short_collapse_canceled.polylines.txt"); -#endif + bool execute_operation(const Element_type& edge, C3t3& c3t3) override + { + Short_edges no_short_edges; // no work list to keep up to date + return execute_operation(edge, c3t3, no_short_edges); + } - while(!short_edges.empty()) + /** + * Collapses `edge`, and keeps `short_edges` up to date : `collapse_edge()` + * removes from it the edges it destroys, and the edges incident to the + * vertex it keeps are re-evaluated here, since their length has changed. + */ + bool execute_operation(const Element_type& edge, C3t3& c3t3, + Short_edges& short_edges) { - //the edge with shortest length - typename Boost_bimap::right_map::iterator eit = short_edges.right.begin(); - Edge e = eit->second; - -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS - FT sqlen = eit->first; - std::cout << "\rCollapse... (" << short_edges.left.size() << " short edges, "; - std::cout << std::sqrt(sqlen) << ", "; - std::cout << nb_collapses << " collapses)"; - std::cout.flush(); -#endif + const Vertex_handle vh = collapse_edge(edge, c3t3, m_sizing, m_protect_boundaries, + m_cell_selector, short_edges, m_visitor); + if (vh == Vertex_handle()) + return false; - short_edges.right.erase(eit); + std::vector incident_short; + c3t3.triangulation().finite_incident_edges(vh, std::back_inserter(incident_short)); + for (const Edge& eshort : incident_short) + { + const auto [collapsible, boundary] + = can_be_collapsed(eshort, c3t3, m_protect_boundaries, m_cell_selector); - CGAL_expensive_assertion_code(const bool bd = is_boundary_edge(e)); - CGAL_expensive_assertion(!!is_too_short(e, bd, sizing, c3t3, cell_selector)); - CGAL_expensive_assertion(can_be_collapsed(e, c3t3, protect_boundaries, cell_selector)); + // an edge that can no longer be collapsed leaves the work list, rather + // than being taken out of it later and refused by collapse_edge() + std::optional sqlen; + if (collapsible) + sqlen = is_too_short(eshort, boundary, m_sizing, c3t3, m_cell_selector); -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - const auto p1 = e.first->vertex(e.second)->point(); - const auto p2 = e.first->vertex(e.third)->point(); -#endif + update_bimap(eshort, short_edges, sqlen); + } + return true; + } - Vertex_handle vh = collapse_edge(e, c3t3, sizing, - protect_boundaries, cell_selector, - short_edges, - visitor); - if (vh != Vertex_handle()) - { - std::vector incident_short; - c3t3.triangulation().finite_incident_edges(vh, - std::back_inserter(incident_short)); - for (const Edge& eshort : incident_short) - { - const auto [collapsible, boundary] - = can_be_collapsed(eshort, c3t3, protect_boundaries, cell_selector); - if (!collapsible) - continue; + std::string operation_name() const override { return "Collapse short edges"; } +}; - const auto sqlen = is_too_short(eshort, boundary, sizing, c3t3, cell_selector); - update_bimap(eshort, short_edges, sqlen); - } +/** +* Collapse is the only operation whose elements change as it runs : collapsing +* an edge shortens the ones around the vertex it keeps, and destroys others. +* Its elements are therefore taken from a work list that `execute_operation()` +* keeps up to date, shortest first, rather than from a list collected once. +*/ +template +class Elementary_operation_execution_sequential< + Edge_collapse_operation > +{ + using Operation = Edge_collapse_operation; + using Short_edges = typename Operation::Short_edges; + using Edge = typename Operation::Edge; + +public: + bool execute(Operation& op, C3t3& c3t3) const + { + Short_edges short_edges = op.get_elements(c3t3); + if (short_edges.empty()) + return false; #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - ++nb_collapses; + std::size_t nb_done = 0; + CGAL::Real_timer timer; + timer.start(); #endif + while (!short_edges.empty()) + { + // the edge with shortest length + typename Short_edges::right_map::iterator eit = short_edges.right.begin(); + const Edge e = eit->second; + short_edges.right.erase(eit); -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - if (vh != Vertex_handle()) - short_success << "2 " << point(p1) << " " << point(p2) << std::endl; - else - short_fail << "2 " << point(p1) << " " << point(p2) << std::endl; + if (op.execute_operation(e, c3t3, short_edges)) + { +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + ++nb_done; #endif + } } - }//end loop on short_edges -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - short_success.close(); - short_fail.close(); -#endif #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - timer.stop(); - std::cout << " done (" << nb_collapses << " collapses, in " - << timer.time() << " seconds)." << std::endl; + timer.stop(); + std::cout << op.operation_name() << ": " << nb_done << " done (" + << timer.time() << " sec)." << std::endl; #endif -} + return true; + } +}; + } } } diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h index 3e9f1902764..dc2e02df630 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h @@ -18,6 +18,7 @@ #include #include +#include #include #include @@ -28,10 +29,6 @@ #include #include -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE -#include -#endif - namespace CGAL { namespace Tetrahedral_remeshing @@ -1396,8 +1393,8 @@ template Sliver_removal_result flip_on_surface(C3T3& c3t3, typename C3T3::Edge& edge, - typename C3T3::Vertex_handle v0i,//v0 of new edge that will replace edge - typename C3T3::Vertex_handle v1i,//v1 of new edge that will replace edge + const typename C3T3::Vertex_handle v0i,//v0 of new edge that will replace edge + const typename C3T3::Vertex_handle v1i,//v1 of new edge that will replace edge IncCellsVectorMap& inc_cells, Flip_Criterion flip_criterion, Visitor& visitor) @@ -1755,20 +1752,18 @@ Sliver_removal_result flip_on_surface(C3T3& c3t3, template -std::size_t flipBoundaryEdges( - C3T3& c3t3, - const std::vector& boundary_edges, - SurfaceIndexMapMap& boundary_vertices_valences, - IncidentCellsVectorMap& inc_cells, - const Flip_Criterion& flip_criterion, - CellSelector& cell_selector, - Visitor& visitor) +std::size_t flipBoundaryEdges(C3T3& c3t3, + const std::vector& boundary_edges, + SurfaceIndexMapMap& boundary_vertices_valences, + IncidentCellsVectorMap& inc_cells, + const Flip_Criterion& flip_criterion, + CellSelector& cell_selector, + Visitor& visitor) { typedef typename C3T3::Vertex_handle Vertex_handle; - typedef typename C3T3::Cell_handle Cell_handle; - typedef typename C3T3::Facet Facet; - typedef typename C3T3::Edge Edge; - typedef typename C3T3::Surface_patch_index Surface_patch_index; + typedef typename C3T3::Cell_handle Cell_handle; + typedef typename C3T3::Facet Facet; + typedef typename C3T3::Edge Edge; typedef typename C3T3::Triangulation Tr; typedef std::pair Edge_vv; @@ -1777,56 +1772,88 @@ std::size_t flipBoundaryEdges( Tr& tr = c3t3.triangulation(); std::vector candidate_edges_for_flip; - for (const Edge& e : boundary_edges) - { - if (!c3t3.is_in_complex(e)) + for(const Edge& e : boundary_edges) { + if(!c3t3.is_in_complex(e)) candidate_edges_for_flip.push_back(make_vertex_pair(e)); } - for (const auto& [vh0, vh1] : candidate_edges_for_flip) - { + for(const auto& [vh0, vh1] : candidate_edges_for_flip) { boost::container::small_vector& inc_vh0 = inc_cells[vh0]; - if (inc_vh0.empty()) + if(inc_vh0.empty()) tr.incident_cells(vh0, std::back_inserter(inc_vh0)); Cell_handle c; int i, j; - if (!is_edge_uv(vh0, vh1, inc_vh0, c, i, j)) + if(!is_edge_uv(vh0, vh1, inc_vh0, c, i, j)) continue; Edge edge(c, i, j); std::vector boundary_facets; const bool on_boundary = is_boundary_edge(edge, c3t3, cell_selector, boundary_facets); -// if (on_boundary && boundary_facets.empty()) -// { -// std::cerr << vh0->point().point() << "\t " << vh1->point().point() << std::endl; -// bool b = is_boundary_edge(vh0, vh1, c3t3, cell_selector); -// CGAL::Tetrahedral_remeshing::debug::dump_c3t3(c3t3, "dump_c3t3_about_boundary_"); -// CGAL::Tetrahedral_remeshing::debug::dump_facets_in_complex(c3t3, "dump_facets_about_boundary_.off"); -// CGAL::Tetrahedral_remeshing::debug::dump_facets_from_selection( -// c3t3, cell_selector, "dump_facets_from_selection_.off"); -// std::cerr << "valid = " << tr.tds().is_valid(true) << std::endl; -// std::cerr << "boundary = " << b << std::endl; -// CGAL_assertion(on_boundary); -// } -// else if (on_boundary && boundary_facets.size() != 2) -// { -// std::cerr << vh0->point().point() << "\t " << vh1->point().point() << std::endl; -// CGAL::Tetrahedral_remeshing::debug::dump_c3t3(c3t3, "dump_c3t3_about_boundary_"); -// CGAL::Tetrahedral_remeshing::debug::dump_facets(boundary_facets, "dump_boundary_facets.polylines.txt"); -// std::vector dummy_facets; -// bool b = is_boundary_edge(edge, c3t3, cell_selector, dummy_facets, true/**/); -// std::cerr << "boundary = " << b << std::endl; -// } - - if (!on_boundary || boundary_facets.size() != 2) + // if (on_boundary && boundary_facets.empty()) + // { + // std::cerr << vh0->point().point() << "\t " << vh1->point().point() << std::endl; + // bool b = is_boundary_edge(vh0, vh1, c3t3, cell_selector); + // CGAL::Tetrahedral_remeshing::debug::dump_c3t3(c3t3, "dump_c3t3_about_boundary_"); + // CGAL::Tetrahedral_remeshing::debug::dump_facets_in_complex(c3t3, "dump_facets_about_boundary_.off"); + // CGAL::Tetrahedral_remeshing::debug::dump_facets_from_selection( + // c3t3, cell_selector, "dump_facets_from_selection_.off"); + // std::cerr << "valid = " << tr.tds().is_valid(true) << std::endl; + // std::cerr << "boundary = " << b << std::endl; + // CGAL_assertion(on_boundary); + // } + // else if (on_boundary && boundary_facets.size() != 2) + // { + // std::cerr << vh0->point().point() << "\t " << vh1->point().point() << std::endl; + // CGAL::Tetrahedral_remeshing::debug::dump_c3t3(c3t3, "dump_c3t3_about_boundary_"); + // CGAL::Tetrahedral_remeshing::debug::dump_facets(boundary_facets, "dump_boundary_facets.polylines.txt"); + // std::vector dummy_facets; + // bool b = is_boundary_edge(edge, c3t3, cell_selector, dummy_facets, true/**/); + // std::cerr << "boundary = " << b << std::endl; + // } + + if(!on_boundary || boundary_facets.size() != 2) continue; CGAL_assertion(boundary_facets.size() == 2); + if(flip_surface_edge(c3t3, edge, boundary_facets, boundary_vertices_valences, + inc_cells, flip_criterion, visitor)) + { + ++nb_success; + CGAL_expensive_assertion(tr.tds().is_valid()); + } + } + return nb_success; +} + +//todo : write this function +template +bool flip_surface_edge(C3t3& c3t3, + typename C3t3::Edge& edge, + const std::vector& boundary_facets, + BV_valences& boundary_vertices_valences, + IncidentCellsVectorMap& inc_cells, + const Flip_criterion& flip_criterion, + Visitor& visitor) +{ + using Vertex_handle = typename C3t3::Vertex_handle; + using Facet = typename C3t3::Facet; + using Cell_handle = typename C3t3::Cell_handle; + using Surface_patch_index = typename C3t3::Surface_patch_index; + + auto& tr = c3t3.triangulation(); + const Facet& f0 = boundary_facets[0]; const Facet& f1 = boundary_facets[1]; + const Vertex_handle vh0 = edge.first->vertex(edge.second); + const Vertex_handle vh1 = edge.first->vertex(edge.third); + // find 3rd and 4th vertices to flip on surface const Vertex_handle vh2 = third_vertex(f0, vh0, vh1, tr); const Vertex_handle vh3 = third_vertex(f1, vh0, vh1, tr); @@ -1843,7 +1870,7 @@ std::size_t flipBoundaryEdges( int v3 = boundary_vertices_valences.at(vh3)[surfi]; if(v0 < 2 || v1 < 2 || v2 < 2 || v3 < 2) - continue; + return false; int m0 = (boundary_vertices_valences.at(vh0).size() > 1 ? 4 : 6); int m1 = (boundary_vertices_valences.at(vh1).size() > 1 ? 4 : 6); @@ -1905,109 +1932,202 @@ std::size_t flipBoundaryEdges( boundary_vertices_valences[vh2][surfi]++; boundary_vertices_valences[vh3][surfi]++; - nb_success++; + return true; } - else - continue; } } - } - CGAL_expensive_assertion(tr.tds().is_valid()); - - return nb_success; + return false; } -template -void flip_edges(C3T3& c3t3, - const bool protect_boundaries, - CellSelector& cell_selector, - Visitor& visitor) +// Shared state for the internal and boundary edge-flip operations: the cell +// selector, the visitor, and the incident-cells cache used by find_best_flip +// and flip_on_surface. +template +class Edge_flip_operation_base { - CGAL_USE(protect_boundaries); - typedef typename C3T3::Triangulation T3; - typedef typename T3::Vertex_handle Vertex_handle; - typedef typename T3::Cell_handle Cell_handle; - typedef typename T3::Edge Edge; - typedef typename std::pair Edge_vv; - typedef typename C3T3::Subdomain_index Subdomain_index; - -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - std::cout << "Flip edges..."; - std::cout.flush(); - std::size_t nb_flips_in_volume = 0; - std::size_t nb_flips_on_surface = 0; - CGAL::Real_timer timer; - timer.start(); -#endif - - for (auto c : c3t3.cells_in_complex()) - c->reset_cache_validity();//we will use sliver_value - //to store the cos_dihedral_angle +protected: + using Tr = typename C3t3::Triangulation; + using Cell_handle = typename Tr::Cell_handle; + using Vertex_handle = typename Tr::Vertex_handle; + using Edge = typename Tr::Edge; + using Facet = typename Tr::Facet; + using Cells_vector = boost::container::small_vector; + using Incident_cells_map = std::unordered_map; + + CellSelector& m_cell_selector; + Visitor& m_visitor; + // Shared across the internal and boundary flip passes, exactly as the former + // flip_edges() shared a single inc_cells map between flip_all_edges() and + // flipBoundaryEdges(). + Incident_cells_map& inc_cells; + + Edge_flip_operation_base(CellSelector& cell_selector, + Visitor& visitor, + Incident_cells_map& incident_cells) + : m_cell_selector(cell_selector) + , m_visitor(visitor) + , inc_cells(incident_cells) {} +}; + +// Flip of internal (non-boundary) edges. Mirrors the former flip_all_edges(): +// reset the cell caches, collect the internal edges, then run find_best_flip +// on each in turn, sharing the incident-cells cache. +template +class Internal_edge_flip_operation + : public Edge_flip_operation_base, + public Elementary_operation, + std::vector>> +{ + using BaseClass = Edge_flip_operation_base; + using typename BaseClass::Cell_handle; + using typename BaseClass::Vertex_handle; + using typename BaseClass::Edge; + using typename BaseClass::Cells_vector; + using BaseClass::m_cell_selector; + using BaseClass::m_visitor; + using BaseClass::inc_cells; + +public: + using Incident_cells_map = typename BaseClass::Incident_cells_map; + using Edge_vv = std::pair; + using Base_operation = Elementary_operation>; + using Element_type = typename Base_operation::Element_type; + static_assert(std::is_same_v, "Element_type must be Edge_vv"); + using Element_range = typename Base_operation::Element_range; + + Internal_edge_flip_operation(CellSelector& cell_selector, + Visitor& visitor, + Incident_cells_map& incident_cells) + : BaseClass(cell_selector, visitor, incident_cells) {} + + Element_range get_elements(const C3t3& c3t3) const override + { + for (auto c : c3t3.cells_in_complex()) + c->reset_cache_validity();//we will use sliver_value + //to store the cos_dihedral_angle - //const Flip_Criterion criterion = VALENCE_MIN_DH_BASED; + std::vector inside_edges; + get_internal_edges(c3t3, m_cell_selector, std::back_inserter(inside_edges)); + return inside_edges; + } - std::vector inside_edges; - get_internal_edges(c3t3, - cell_selector, - std::back_inserter(inside_edges)); + bool execute_operation(const Element_type& vp, C3t3& c3t3) override + { + Cells_vector& o_inc_vh = inc_cells[vp.first]; + if (o_inc_vh.empty()) + c3t3.triangulation().incident_cells(vp.first, std::back_inserter(o_inc_vh)); -// //if (criterion == VALENCE_BASED) -// // flip_inside_edges(inside_edges); -// //else -// //{ -//#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE -// nb_flips += -//#endif -// flip_all_edges(inside_edges, c3t3, MIN_ANGLE_BASED, visitor); -// //} + Cell_handle ch; + int i0, i1; + if (!is_edge_uv(vp.first, vp.second, o_inc_vh, ch, i0, i1)) + return false; - std::unordered_map > inc_cells; + Edge edge(ch, i0, i1); + const Sliver_removal_result res + = find_best_flip(edge, c3t3, MIN_ANGLE_BASED, inc_cells, m_cell_selector, m_visitor); + return (res == VALID_FLIP); + } -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - nb_flips_in_volume += -#endif - flip_all_edges(inside_edges, c3t3, inc_cells, MIN_ANGLE_BASED, cell_selector, visitor); - if (!protect_boundaries) + std::string operation_name() const override { return "Flip edges (internal)"; } +}; + +// Flip of boundary edges. Mirrors the former flipBoundaryEdges(): compute the +// per-vertex boundary valences once, then, for each boundary edge, flip on the +// surface when it lowers the valence cost, updating the valences accordingly. +template +class Boundary_edge_flip_operation + : public Edge_flip_operation_base, + public Elementary_operation, + std::vector>> +{ + using BaseClass = Edge_flip_operation_base; + using typename BaseClass::Cell_handle; + using typename BaseClass::Vertex_handle; + using typename BaseClass::Edge; + using typename BaseClass::Facet; + using typename BaseClass::Cells_vector; + using BaseClass::m_cell_selector; + using BaseClass::m_visitor; + using BaseClass::inc_cells; + + using Subdomain_index = typename C3t3::Subdomain_index; + using Surface_patch_index = typename C3t3::Surface_patch_index; + using Spi_map = boost::unordered_map; + + mutable boost::unordered_map m_boundary_vertices_valences; + +public: + using Incident_cells_map = typename BaseClass::Incident_cells_map; + using Edge_vv = std::pair; + using Base_operation = Elementary_operation>; + using Element_type = typename Base_operation::Element_type; + static_assert(std::is_same_v, "Element_type must be Edge_vv"); + using Element_range = typename Base_operation::Element_range; + + Boundary_edge_flip_operation(CellSelector& cell_selector, Visitor& visitor, Incident_cells_map& incident_cells) + : BaseClass(cell_selector, visitor, incident_cells) {} + + Element_range get_elements(const C3t3& c3t3) const override { - typedef typename C3T3::Surface_patch_index Surface_patch_index; - typedef boost::unordered_map Spi_map; - - //Boundary flip std::vector boundary_edges; - boost::unordered_map boundary_vertices_valences; - boost::unordered_map > vertices_subdomain_indices; + boost::unordered_map> vertices_subdomain_indices; + m_boundary_vertices_valences.clear(); collectBoundaryEdgesAndComputeVerticesValences(c3t3, - cell_selector, + m_cell_selector, boundary_edges, - boundary_vertices_valences, + m_boundary_vertices_valences, vertices_subdomain_indices); #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - if(!debug::are_cell_orientations_valid(c3t3.triangulation())) + if (!debug::are_cell_orientations_valid(c3t3.triangulation())) std::cerr << "ERROR in ORIENTATION" << std::endl; #endif - // if (criterion == VALENCE_BASED) - // flipBoundaryEdges(boundary_edges, boundary_vertices_valences, VALENCE_BASED); - // else -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - nb_flips_on_surface += -#endif - flipBoundaryEdges(c3t3, boundary_edges, boundary_vertices_valences, - inc_cells, - MIN_ANGLE_BASED, - cell_selector, visitor); + std::vector candidate_edges_for_flip; + for (const Edge& e : boundary_edges) + { + if (!c3t3.is_in_complex(e)) + candidate_edges_for_flip.push_back(make_vertex_pair(e)); + } + return candidate_edges_for_flip; } -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - timer.stop(); - std::cout << "\rFlip edges... done (" - << nb_flips_on_surface << "/" - << nb_flips_in_volume << " surface/volume flips done, in " - << timer.time() << " seconds)." << std::endl; -#endif -} + bool execute_operation(const Element_type& vp, C3t3& c3t3) override + { + const Vertex_handle vh0 = vp.first; + const Vertex_handle vh1 = vp.second; + typename C3t3::Triangulation& tr = c3t3.triangulation(); + + Cells_vector& inc_vh0 = inc_cells[vh0]; + if (inc_vh0.empty()) + tr.incident_cells(vh0, std::back_inserter(inc_vh0)); + + Cell_handle c; + int i, j; + if (!is_edge_uv(vh0, vh1, inc_vh0, c, i, j)) + return false; + + Edge edge(c, i, j); + std::vector boundary_facets; + const bool on_boundary = is_boundary_edge(edge, c3t3, m_cell_selector, boundary_facets); + + if (!on_boundary || boundary_facets.size() != 2) + return false; + + return flip_surface_edge(c3t3, edge, + boundary_facets, + m_boundary_vertices_valences, + inc_cells, + MIN_ANGLE_BASED, + m_visitor); + } + + std::string operation_name() const override { return "Flip edges (boundary)"; } +}; }//namespace internal }//namespace Tetrahedral_remeshing diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h index cf73ab88891..9923a9eb39e 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h @@ -19,6 +19,7 @@ #include #include +#include #include #include @@ -41,20 +42,21 @@ namespace Tetrahedral_remeshing { namespace internal { -template -class Tetrahedral_remeshing_smoother +template +class Vertex_smoothing_context { - typedef typename C3t3::Triangulation Tr; - typedef typename C3t3::Surface_patch_index Surface_patch_index; - typedef typename Tr::Cell_handle Cell_handle; - typedef typename Tr::Vertex_handle Vertex_handle; - typedef typename Tr::Edge Edge; - typedef typename Tr::Facet Facet; - - typedef typename Tr::Geom_traits Gt; - typedef typename Gt::Vector_3 Vector_3; - typedef typename Gt::Point_3 Point_3; - typedef typename Gt::FT FT; +public: + using Tr = typename C3t3::Triangulation; + using Surface_patch_index = typename C3t3::Surface_patch_index; + using Cell_handle = typename Tr::Cell_handle; + using Vertex_handle = typename Tr::Vertex_handle; + using Edge = typename Tr::Edge; + using Facet = typename Tr::Facet; + + using Gt = typename Tr::Geom_traits; + using Vector_3 = typename Gt::Vector_3; + using Point_3 = typename Gt::Point_3; + using FT = typename Gt::FT; using Triangle_vec = std::vector; using Triangle_iter = typename Triangle_vec::iterator; @@ -69,19 +71,35 @@ class Tetrahedral_remeshing_smoother using AABB_segment_tree = CGAL::AABB_tree; private: - typedef CGAL::Tetrahedral_remeshing::internal::FMLS FMLS; + Triangle_vec m_aabb_triangles; + Segment_vec m_aabb_segments; + +public: + using FMLS = CGAL::Tetrahedral_remeshing::internal::FMLS; std::vector subdomain_FMLS; std::unordered_map> subdomain_FMLS_indices; - Triangle_vec m_aabb_triangles; +public: AABB_triangle_tree m_triangles_aabb_tree; - Segment_vec m_aabb_segments; AABB_segment_tree m_segments_aabb_tree; FT m_aabb_epsilon; const SizingFunction& m_sizing; + + using Incident_cells_vector = boost::container::small_vector; + std::vector m_inc_cells; + + using Vertices_surface_indices_map = std::unordered_map>; + using Vertices_normals_map = + std::unordered_map>>; + + Vertices_surface_indices_map m_vertices_surface_indices; + Vertices_normals_map m_vertices_normals; + const CellSelector& m_cell_selector; const bool m_protect_boundaries; + const bool m_smooth_constrained_edges; // the 2 following variables become useful and valid @@ -91,48 +109,55 @@ class Tetrahedral_remeshing_smoother std::vector m_free_vertices{}; bool m_flip_smooth_steps{false}; +public: struct Move { Vector_3 move; int neighbors; FT mass; }; - -public: - Tetrahedral_remeshing_smoother(const SizingFunction& sizing, - const CellSelector& cell_selector, - const bool protect_boundaries, - const bool smooth_constrained_edges) - : m_sizing(sizing) - , m_cell_selector(cell_selector) - , m_protect_boundaries(protect_boundaries) - , m_smooth_constrained_edges(smooth_constrained_edges) - {} - - void init(const C3t3& c3t3) + std::vector m_moves{}; + FT m_total_move{0}; + + Vertex_smoothing_context(C3t3& c3t3, + const SizingFunction& sizing, + const CellSelector& cell_selector, + const bool protect_boundaries, + const bool smooth_constrained_edges) + : m_sizing(sizing) + , m_cell_selector(cell_selector) + , m_protect_boundaries(protect_boundaries) + , m_smooth_constrained_edges(smooth_constrained_edges) { + refresh(c3t3); #ifdef CGAL_TET_REMESHING_SMOOTHING_WITH_MLS - //collect a map of vertices surface indices - std::unordered_map > vertices_surface_indices; - collect_vertices_surface_indices(c3t3, vertices_surface_indices); - - //collect a map of normals at surface vertices - std::unordered_map>> vertices_normals; - compute_vertices_normals(c3t3, vertices_normals); - - // Build MLS Surfaces + if (m_protect_boundaries) + { + collect_vertices_surface_indices(c3t3); + compute_vertices_normals(c3t3); + } createMLSSurfaces(subdomain_FMLS, subdomain_FMLS_indices, - vertices_normals, - vertices_surface_indices, + m_vertices_normals, + m_vertices_surface_indices, c3t3); #else - // Build AABB tree build_aabb_trees(c3t3); #endif } + void refresh(C3t3& c3t3) + { + if (!m_protect_boundaries) + { + collect_vertices_surface_indices(c3t3); + compute_vertices_normals(c3t3); + } + reset_vertex_id_map(c3t3.triangulation()); + reset_free_vertices(c3t3.triangulation()); + collect_incident_cells(c3t3.triangulation()); + } + void start_flip_smooth_steps(const C3t3& c3t3) { CGAL_assertion(!m_flip_smooth_steps); @@ -145,26 +170,28 @@ class Tetrahedral_remeshing_smoother m_flip_smooth_steps = true; } -private: - - bool is_selected(const Cell_handle c) const - { - return get(m_cell_selector, c); - } - - std::size_t vertex_id(const Vertex_handle v) const + std::size_t vertex_id(const Vertex_handle v) const { CGAL_expensive_assertion(m_vertex_id.find(v) != m_vertex_id.end()); return m_vertex_id.at(v); } - bool is_free(const Vertex_handle v) const + bool is_free(const Vertex_handle v) const { return m_free_vertices[vertex_id(v)]; } + bool is_free(const std::size_t& vid) const { return m_free_vertices[vid]; } + + const Incident_cells_vector& incident_cells(const Vertex_handle v) const { - return m_free_vertices[vertex_id(v)]; + return m_inc_cells[vertex_id(v)]; } - bool is_free(const std::size_t & vid) const + const Incident_cells_vector& incident_cells(const std::size_t& vid) const + { + return m_inc_cells[vid]; + } + +private: + bool is_selected(const Cell_handle c) const { - return m_free_vertices[vid]; + return get(m_cell_selector, c); } // this function can be used iff m_vertex_id @@ -240,29 +267,22 @@ class Tetrahedral_remeshing_smoother bb.zmax() - bb.zmin())); } - template - void collect_incident_cells(const Tr& tr, - IncCellsVector& inc_cells) + void collect_incident_cells(const Tr& tr) { + m_inc_cells.clear(); + const std::size_t nbv = tr.number_of_vertices(); + m_inc_cells.resize(nbv, Incident_cells_vector{}); for (const Cell_handle c : tr.finite_cell_handles()) { for (auto vi : tr.vertices(c)) { const std::size_t idi = vertex_id(vi); if(is_free(idi)) - inc_cells[idi].push_back(c); + m_inc_cells[idi].push_back(c); } } } - Point_3 project_on_tangent_plane(const Point_3& gi, - const Point_3& pi, - const Vector_3& normal) - { - Vector_3 diff(gi, pi); - return gi + (normal * diff) * normal; - } - std::optional find_adjacent_facet_on_surface(const Facet& f, const Edge& edge, @@ -296,14 +316,14 @@ class Tetrahedral_remeshing_smoother return {}; } - template + template Vector_3 compute_normal(const Facet& f, const Vector_3& reference_normal, - const Gt& gt) + const Gt_& gt) { - typename Gt::Construct_opposite_vector_3 + typename Gt_::Construct_opposite_vector_3 opp = gt.construct_opposite_vector_3_object(); - typename Gt::Compute_scalar_product_3 + typename Gt_::Compute_scalar_product_3 scalar_product = gt.compute_scalar_product_3_object(); Vector_3 n = CGAL::Tetrahedral_remeshing::normal(f, gt); @@ -327,10 +347,9 @@ class Tetrahedral_remeshing_smoother return str; } - template - void compute_vertices_normals(const C3t3& c3t3, - VertexNormalsMap& normals_map) + void compute_vertices_normals(const C3t3& c3t3) { + m_vertices_normals.clear(); typename Tr::Geom_traits gt = c3t3.triangulation().geom_traits(); typename Tr::Geom_traits::Construct_opposite_vector_3 opp = gt.construct_opposite_vector_3_object(); @@ -400,16 +419,16 @@ class Tetrahedral_remeshing_smoother for (const Vertex_handle vi : tr.vertices(f)) { - typename VertexNormalsMap::iterator patch_vector_it = normals_map.find(vi); + typename Vertices_normals_map::iterator patch_vector_it = m_vertices_normals.find(vi); - if (patch_vector_it == normals_map.end() + if (patch_vector_it == m_vertices_normals.end() || patch_vector_it->second.find(surf_i) == patch_vector_it->second.end()) { - normals_map[vi][surf_i] = n; + m_vertices_normals[vi][surf_i] = n; } else { - normals_map[vi][surf_i] += n; + m_vertices_normals[vi][surf_i] += n; } } } @@ -422,24 +441,20 @@ class Tetrahedral_remeshing_smoother #endif //normalize the computed normals - for (typename VertexNormalsMap::iterator vnm_it = normals_map.begin(); - vnm_it != normals_map.end(); ++vnm_it) + for (auto& [v, patch_normals] : m_vertices_normals) { //value type is map - for (typename VertexNormalsMap::mapped_type::iterator it = vnm_it->second.begin(); - it != vnm_it->second.end(); ++it) + for (auto& [surf_i, n] : patch_normals) { - Vector_3& n = it->second; - #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - auto p = point(vnm_it->first->point()); + auto p = point(v->point()); os << "2 " << p << " " << (p + n) << std::endl; #endif CGAL::Tetrahedral_remeshing::normalize(n, gt); #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - const Surface_patch_index si = it->first; + const Surface_patch_index si = surf_i; if (ons_map.find(si) == ons_map.end()) ons_map[si] = std::vector(); ons_map[si].push_back(typename Tr::Geom_traits::Segment_3(p, p + n)); @@ -462,42 +477,91 @@ class Tetrahedral_remeshing_smoother #endif } - std::optional project(const Surface_patch_index& si, - const Point_3& gi) + void collect_vertices_surface_indices(const C3t3& c3t3) { - CGAL_expensive_assertion(subdomain_FMLS_indices.find(si) != subdomain_FMLS_indices.end()); - CGAL_assertion(!std::isnan(gi.x()) && !std::isnan(gi.y()) && !std::isnan(gi.z())); + m_vertices_surface_indices.clear(); + for (Facet fit : c3t3.facets_in_complex()) + { + const Surface_patch_index& surface_index = c3t3.surface_patch_index(fit); - Vector_3 point(gi.x(), gi.y(), gi.z()); - Vector_3 res_normal = CGAL::NULL_VECTOR; - Vector_3 result(CGAL::ORIGIN, gi); + for (const Vertex_handle vi : c3t3.triangulation().vertices(fit)) + { + std::vector& v_surface_indices = m_vertices_surface_indices[vi]; + if (std::find(v_surface_indices.begin(), v_surface_indices.end(), surface_index) == v_surface_indices.end()) + v_surface_indices.push_back(surface_index); + } + } + } - const FMLS& fmls = subdomain_FMLS[subdomain_FMLS_indices.at(si)]; + void reset_vertex_id_map(const Tr& tr) + { + // when flip-smooth steps start, + // m_vertex_id should already be initialized, + // done by the last smoothing step. + // Then, it does not need to be recomputed + // because no vertices are inserted nor removed anymore + if(m_flip_smooth_steps) + return; + m_vertex_id.clear(); + std::size_t id = 0; + for (const Vertex_handle v : tr.finite_vertex_handles()) + { + m_vertex_id[v] = id++; + } + } +}; - int it_nb = 0; - const int max_it_nb = 5; - const double epsilon = fmls.getPNScale() / 1000.; - const double sq_eps = CGAL::square(epsilon); - do - { - point = result; +template +class Vertex_smooth_operation_base +{ +protected: + typedef typename C3t3::Triangulation Tr; + typedef typename C3t3::Surface_patch_index Surface_patch_index; + typedef typename Tr::Cell_handle Cell_handle; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Edge Edge; + typedef typename Tr::Facet Facet; - fmls.fastProjectionCPU(point, result, res_normal); + typedef typename Tr::Geom_traits Gt; + typedef typename Gt::Vector_3 Vector_3; + typedef typename Gt::Point_3 Point_3; + typedef typename Gt::FT FT; - if (std::isnan(result[0]) || std::isnan(result[1]) || std::isnan(result[2])) { - std::cout << "MLS error detected si " //<< si - << "\t(size : " << fmls.getPNSize() << ")" - << "\t(point = " << point << " )" << std::endl; - return {}; - } - } while ((result - point).squared_length() > sq_eps && ++it_nb < max_it_nb); + using Triangle_vec = std::vector; + using Triangle_iter = typename Triangle_vec::iterator; + using Triangle_primitive = CGAL::AABB_triangle_primitive_3; + using AABB_triangle_traits = CGAL::AABB_traits_3; + using AABB_triangle_tree = CGAL::AABB_tree; - return Point_3(result[0], result[1], result[2]); + using Context = Vertex_smoothing_context; + +public: + std::shared_ptr m_context{nullptr}; + + Vertex_smooth_operation_base(std::shared_ptr context) + : m_context(context) {} + + void set_context(std::shared_ptr p_context) { m_context = p_context; } + +protected: + Point_3 project_on_tangent_plane(const Point_3& gi, const Point_3& pi, const Vector_3& normal) + { + Vector_3 diff(gi, pi); + return gi + (normal * diff) * normal; } - Dihedral_angle_cosine max_cosine(const Tr& tr, - const boost::container::small_vector& cells) + FT density_along_segment(const Edge& e, const C3t3& c3t3, bool boundary_edge = false) const + { + const auto [pt, dim, index] = midpoint_with_info(e, boundary_edge, c3t3); + const FT s = sizing_at_midpoint(e, pt, dim, index, m_context->m_sizing, c3t3, m_context->m_cell_selector); + return 1. / s; + } + + bool is_selected(const Cell_handle c) const { return get(m_context->m_cell_selector, c); } + + template + Dihedral_angle_cosine max_cosine(const Tr& tr, const CellRange& cells) const { Dihedral_angle_cosine max_cos_dh = cosine_of_90_degrees();// = 0. for (Cell_handle c : cells) @@ -519,12 +583,12 @@ class Tetrahedral_remeshing_smoother const CellRange& inc_cells, const Tr& tr, #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - FT& total_move) + FT& total_move) const #else - FT&) + FT&) const #endif { - const typename Tr::Point backup = v->point(); //backup v's position + const typename Tr::Point backup = v->point();//backup v's position const typename Tr::Geom_traits::Point_3 pv = point(backup); bool valid_orientation = false; @@ -532,7 +596,7 @@ class Tetrahedral_remeshing_smoother double frac = 1.0; typename Tr::Geom_traits::Vector_3 move(pv, final_pos); - const Dihedral_angle_cosine curr_max_cos = m_flip_smooth_steps + const Dihedral_angle_cosine curr_max_cos = m_context->m_flip_smooth_steps ? max_cosine(tr, inc_cells) : Dihedral_angle_cosine(CGAL::ZERO, 0., 1.);//Dummy unused value @@ -557,7 +621,7 @@ class Tetrahedral_remeshing_smoother valid_orientation = false; break; } - else if (m_flip_smooth_steps) //check that dihedral angles get improved + else if (m_context->m_flip_smooth_steps) //check that dihedral angles get improved { if(is_selected(ci)) { @@ -578,7 +642,7 @@ class Tetrahedral_remeshing_smoother while(!valid_try && frac > 0.1); // if move failed, cancel move - bool valid_move = valid_orientation && angles_improved; + bool valid_move = valid_orientation && angles_improved; if(!valid_move) v->set_point(backup); @@ -590,80 +654,101 @@ class Tetrahedral_remeshing_smoother return valid_move; } +}; + +template +class Complex_edge_vertex_smooth_operation + : public Vertex_smooth_operation_base, + public Elementary_operation +{ +public: + using BaseClass = Vertex_smooth_operation_base; + using Vertex_handle = typename C3t3::Triangulation::Vertex_handle; + using Surface_patch_index = typename C3t3::Surface_patch_index; + + using Base_operation = Elementary_operation; + using Element_type = typename Base_operation::Element_type; + static_assert(std::is_same_v, "Element_type should be Vertex_handle"); + using Element_range = typename Base_operation::Element_range; + + using BaseClass::m_context; + + using typename BaseClass::Cell_handle; + using typename BaseClass::Edge; + using typename BaseClass::FT; + using typename BaseClass::Point_3; + using typename BaseClass::Tr; + using typename BaseClass::Vector_3; - void collect_vertices_surface_indices( - const C3t3& c3t3, - std::unordered_map >& vertices_surface_indices) - { - for (Facet fit : c3t3.facets_in_complex()) - { - const Surface_patch_index& surface_index = c3t3.surface_patch_index(fit); +public: + Complex_edge_vertex_smooth_operation(std::shared_ptr context) + : BaseClass(context) {} - for (const Vertex_handle vi : c3t3.triangulation().vertices(fit)) - { - std::vector& v_surface_indices = vertices_surface_indices[vi]; - if (std::find(v_surface_indices.begin(), v_surface_indices.end(), surface_index) == v_surface_indices.end()) - v_surface_indices.push_back(surface_index); - } - } + Element_range get_elements(const C3t3& c3t3) const override + { + perform_global_preprocessing(c3t3); + return c3t3.triangulation().finite_vertex_handles(); } - void reset_vertex_id_map(const Tr& tr) + bool execute_operation(const Element_type& v, C3t3& c3t3) override { - // when flip-smooth steps start, - // m_vertex_id should already be initialized, - // done by the last smoothing step. - // Then, it does not need to be recomputed - // because no vertices are inserted nor removed anymore - if(m_flip_smooth_steps) - return; - m_vertex_id.clear(); - std::size_t id = 0; - for (const Vertex_handle v : tr.finite_vertex_handles()) - { - m_vertex_id[v] = id++; - } - } + auto& tr = c3t3.triangulation(); + const std::size_t vid = m_context->vertex_id(v); + if (!m_context->is_free(vid) || !is_on_feature(v)) + return false; + + const Point_3 current_pos = point(v->point()); + const auto& moves = m_context->m_moves; + const std::size_t nb_neighbors = moves[vid].neighbors; + if (nb_neighbors == 0) + return false; - // boundary_edge is set to false by default because - // we may not care about this information, for example while collecting - // weights in the smoothing inside volume step - FT density_along_segment(const Edge& e, - const C3t3& c3t3, - const bool boundary_edge = false) const - { - const auto [pt, dim, index] = midpoint_with_info(e, boundary_edge, c3t3); - const FT s = sizing_at_midpoint(e, pt, dim, index, m_sizing, c3t3, m_cell_selector); - const FT density = 1. / s; //density = 1 / size^(dimension) - //edge dimension is 1, so density = 1 / size - //to have mass = length * density with no dimension - return density; - } + CGAL_assertion(moves[vid].mass > 0); + const Vector_3 move = (nb_neighbors > 0) + ? moves[vid].move / moves[vid].mass + : CGAL::NULL_VECTOR; + const Point_3 smoothed_position = current_pos + move; - template - std::size_t smooth_edges_in_complex(C3t3& c3t3, #ifdef CGAL_TET_REMESHING_SMOOTHING_WITH_MLS - const SurfaceIndices& vertices_surface_indices, -#else - const SurfaceIndices&, + Vector_3 sum_projections = CGAL::NULL_VECTOR; + Point_3 tmp_pos = current_pos; + +#ifndef CGAL_TET_REMESHING_EDGE_SMOOTHING_DISABLE_PROJECTION + const std::vector& v_surface_indices = m_context->m_vertices_surface_indices.at(v); + for (const Surface_patch_index& si : v_surface_indices) + { + Point_3 normal_projection = BaseClass::project_on_tangent_plane(smoothed_position, current_pos, + m_context->m_vertices_normals.at(v).at(si)); + sum_projections += Vector_3(tmp_pos, normal_projection); + tmp_pos = normal_projection; + } #endif - const IncidentCells& inc_cells, - const NormalsMap& vertices_normals, - FT& total_move -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - , std::ofstream& os_surf + + const Point_3 new_pos = current_pos + sum_projections; +#else + const Point_3 new_pos = m_context->m_segments_aabb_tree.closest_point(smoothed_position); #endif - ) + + const auto& inc_cells = m_context->m_inc_cells[vid]; + return BaseClass::check_inversion_and_move(v, new_pos, inc_cells, tr, m_context->m_total_move); + } + + std::string operation_name() const override { return "Vertex Smooth (Complex Edge Vertices)"; } + + void perform_global_preprocessing(const C3t3& c3t3) const { - std::size_t nb_done_1d = 0; auto& tr = c3t3.triangulation(); + auto& moves = m_context->m_moves; const std::size_t nbv = tr.number_of_vertices(); + using Move = typename BaseClass::Context::Move; const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/}; - std::vector moves(nbv, default_move); + moves.assign(nbv, default_move); //collect neighbors for (const Edge& e : c3t3.edges_in_complex()) @@ -674,18 +759,18 @@ class Tetrahedral_remeshing_smoother CGAL_expensive_assertion(is_on_feature(vh0)); CGAL_expensive_assertion(is_on_feature(vh1)); - const std::size_t& i0 = vertex_id(vh0); - const std::size_t& i1 = vertex_id(vh1); + const std::size_t& i0 = m_context->vertex_id(vh0); + const std::size_t& i1 = m_context->vertex_id(vh1); - const bool vh0_moving = is_free(i0); - const bool vh1_moving = is_free(i1); + const bool vh0_moving = m_context->is_free(i0); + const bool vh1_moving = m_context->is_free(i1); if (!vh0_moving && !vh1_moving) continue; const Point_3& p0 = point(vh0->point()); const Point_3& p1 = point(vh1->point()); - const FT density = density_along_segment(e, c3t3, true); + const FT density = BaseClass::density_along_segment(e, c3t3, true); if (vh0_moving) { @@ -700,217 +785,208 @@ class Tetrahedral_remeshing_smoother ++moves[i1].neighbors; } } + } +}; + +template +class Surface_vertex_smooth_operation + : public Vertex_smooth_operation_base, + public Elementary_operation +{ +public: + using BaseClass = Vertex_smooth_operation_base; + using Base_operation = Elementary_operation; + using Element_type = typename Base_operation::Element_type; + static_assert(std::is_same_v, + "Element_type should be Vertex_handle"); + using Element_range = typename Base_operation::Element_range; + + using BaseClass::m_context; + + using typename BaseClass::AABB_triangle_tree; + using typename BaseClass::Cell_handle; + using typename BaseClass::Edge; + using typename BaseClass::FT; + using typename BaseClass::Gt; + using typename BaseClass::Point_3; + using typename BaseClass::Surface_patch_index; + using typename BaseClass::Tr; + using typename BaseClass::Vector_3; + using typename BaseClass::Vertex_handle; - // iterate over vertices and move - for(Vertex_handle v : tr.finite_vertex_handles()) - { - const std::size_t vid = vertex_id(v); - - if (!is_free(vid) || !is_on_feature(v)) - continue; - - const Point_3 current_pos = point(v->point()); - - const std::size_t nb_neighbors = moves[vid].neighbors; - if(nb_neighbors == 0) - continue; - - CGAL_assertion(moves[vid].mass > 0); - const Vector_3 move = (nb_neighbors > 0) - ? moves[vid].move / moves[vid].mass - : CGAL::NULL_VECTOR; - - const Point_3 smoothed_position = current_pos + move; - - CGAL_USE(vertices_normals); -#ifdef CGAL_TET_REMESHING_SMOOTHING_WITH_MLS - - Vector_3 sum_projections = CGAL::NULL_VECTOR; - Point_3 tmp_pos = current_pos; +private: + void perform_global_preprocessing(const C3t3& c3t3) const + { + auto& tr = c3t3.triangulation(); + auto& moves = m_context->m_moves; + using Move = typename BaseClass::Context::Move; + const std::size_t nbv = tr.number_of_vertices(); + const Move default_move{CGAL::NULL_VECTOR, 0/*neighbors*/, 0./*mass*/}; + moves.assign(nbv, default_move); -#ifndef CGAL_TET_REMESHING_EDGE_SMOOTHING_DISABLE_PROJECTION - const std::vector& v_surface_indices = vertices_surface_indices.at(v); - for (const Surface_patch_index& si : v_surface_indices) + for (const Edge& e : tr.finite_edges()) + { + if (!c3t3.is_in_complex(e) && is_boundary(c3t3, e, m_context->m_cell_selector)) { - Point_3 normal_projection = project_on_tangent_plane(smoothed_position, - current_pos, - vertices_normals.at(v).at(si)); + const Vertex_handle vh0 = e.first->vertex(e.second); + const Vertex_handle vh1 = e.first->vertex(e.third); - sum_projections += Vector_3(tmp_pos, normal_projection); - tmp_pos = normal_projection; - } -#endif //CGAL_TET_REMESHING_EDGE_SMOOTHING_DISABLE_PROJECTION + const std::size_t& i0 = m_context->vertex_id(vh0); + const std::size_t& i1 = m_context->vertex_id(vh1); - const Point_3 new_pos = current_pos + sum_projections; + const bool vh0_moving = !is_on_feature(vh0) && m_context->is_free(i0); + const bool vh1_moving = !is_on_feature(vh1) && m_context->is_free(i1); -#else // AABB_tree projection - - const Point_3 new_pos = m_segments_aabb_tree.closest_point(smoothed_position); + if (!vh0_moving && !vh1_moving) + continue; + const Point_3& p0 = point(vh0->point()); + const Point_3& p1 = point(vh1->point()); + const FT density = BaseClass::density_along_segment(e, c3t3, true); -#endif //CGAL_TET_REMESHING_SMOOTHING_WITH_MLS - -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - os_surf << "2 " << current_pos << " " << new_pos << std::endl; -#endif - // move vertex - if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){ - nb_done_1d++; + if (vh0_moving) + { + moves[i0].move += density * Vector_3(p0, p1); + moves[i0].mass += density; + ++moves[i0].neighbors; + } + if (vh1_moving) + { + moves[i1].move += density * Vector_3(p1, p0); + moves[i1].mass += density; + ++moves[i1].neighbors; + } } } - return nb_done_1d; } - -template -std::size_t smooth_vertices_on_surfaces(C3t3& c3t3, - const SurfaceIndices& vertices_surface_indices, - const IncidentCells& inc_cells, - const NormalsMap& vertices_normals, - FT& total_move -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - , std::ofstream& os_surf - , std::ofstream& os_surf0 -#endif - ) -{ - std::size_t nb_done_2d = 0; - auto& tr = c3t3.triangulation(); - - const std::size_t nbv = tr.number_of_vertices(); - const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/}; - std::vector moves(nbv, default_move); - - for (const Edge& e : tr.finite_edges()) + std::optional project(const Surface_patch_index& si, const Point_3& gi) { - if (!c3t3.is_in_complex(e) && is_boundary(c3t3, e, m_cell_selector)) - { - const Vertex_handle vh0 = e.first->vertex(e.second); - const Vertex_handle vh1 = e.first->vertex(e.third); - - const std::size_t& i0 = vertex_id(vh0); - const std::size_t& i1 = vertex_id(vh1); + CGAL_expensive_assertion(m_context->subdomain_FMLS_indices.find(si) != m_context->subdomain_FMLS_indices.end()); + CGAL_assertion(!std::isnan(gi.x()) && !std::isnan(gi.y()) && !std::isnan(gi.z())); - const bool vh0_moving = !is_on_feature(vh0) && is_free(i0); - const bool vh1_moving = !is_on_feature(vh1) && is_free(i1); + Vector_3 point_vec(gi.x(), gi.y(), gi.z()); + Vector_3 res_normal = CGAL::NULL_VECTOR; + Vector_3 result(CGAL::ORIGIN, gi); - if (!vh0_moving && !vh1_moving) - continue; + const typename BaseClass::Context::FMLS& fmls = m_context->subdomain_FMLS[m_context->subdomain_FMLS_indices.at(si)]; - const Point_3& p0 = point(vh0->point()); - const Point_3& p1 = point(vh1->point()); - const FT density = density_along_segment(e, c3t3, true); + int it_nb = 0; + const int max_it_nb = 5; + const double epsilon = fmls.getPNScale() / 1000.; + const double sq_eps = CGAL::square(epsilon); - if (vh0_moving) - { - moves[i0].move += density * Vector_3(p0, p1); - moves[i0].mass += density; - ++moves[i0].neighbors; - } - if (vh1_moving) + do + { + point_vec = result; + fmls.fastProjectionCPU(point_vec, result, res_normal); + if(std::isnan(result[0]) || std::isnan(result[1]) || std::isnan(result[2])) { - moves[i1].move += density * Vector_3(p1, p0); - moves[i1].mass += density; - ++moves[i1].neighbors; +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + std::cout << "MLS error detected si " + << "\t(size : " << fmls.getPNSize() << ")" + << "\t(point = " << point_vec << " )" << std::endl; +#endif + return {}; } - } + } while((result - point_vec).squared_length() > sq_eps && ++it_nb < max_it_nb); + + return Point_3(result[0], result[1], result[2]); } - // iterate over vertices and move - for(Vertex_handle v : tr.finite_vertex_handles()) +public: + Surface_vertex_smooth_operation(std::shared_ptr context) + : BaseClass(context) {} + + Element_range get_elements(const C3t3& c3t3) const override { - const std::size_t vid = vertex_id(v); + perform_global_preprocessing(c3t3); + return c3t3.triangulation().finite_vertex_handles(); + } - if (!is_free(vid) || v->in_dimension() != 2) - continue; + bool execute_operation(const Element_type& v, C3t3& c3t3) override + { + auto& tr = c3t3.triangulation(); + auto& moves = m_context->m_moves; + const std::size_t vid = m_context->vertex_id(v); + if (!m_context->is_free(vid) || v->in_dimension() != 2) + return false; const std::size_t nb_neighbors = moves[vid].neighbors; const Point_3 current_pos = point(v->point()); - const auto& incident_surface_patches = vertices_surface_indices.at(v); + CGAL_assertion(m_context->m_vertices_surface_indices.find(v) != m_context->m_vertices_surface_indices.end()); + const auto& incident_surface_patches = m_context->m_vertices_surface_indices.at(v); + if (incident_surface_patches.size() > 1) - continue; - const Surface_patch_index si = incident_surface_patches[0]; + return false; + const Surface_patch_index si = incident_surface_patches[0]; CGAL_assertion(si != Surface_patch_index()); CGAL_expensive_assertion_code(auto siv = surface_patch_index(v, c3t3)); CGAL_expensive_assertion(si == siv); + Point_3 new_pos; + bool result = false; + if (nb_neighbors > 1) { const Vector_3 move = moves[vid].move / moves[vid].mass; const Point_3 smoothed_position = point(v->point()) + move; #ifdef CGAL_TET_REMESHING_SMOOTHING_WITH_MLS - Point_3 normal_projection = project_on_tangent_plane(smoothed_position, - current_pos, - vertices_normals.at(v).at(si)); + Point_3 normal_projection = BaseClass::project_on_tangent_plane(smoothed_position, current_pos, + m_context->m_vertices_normals.at(v).at(si)); std::optional mls_projection = project(si, normal_projection); - - const Point_3 new_pos = (mls_projection != std::nullopt) - ? *mls_projection - : smoothed_position; - -#else // AABB_tree projection - Point_3 new_pos; - if (m_triangles_aabb_tree.squared_distance(smoothed_position) < m_aabb_epsilon) + new_pos = (mls_projection != std::nullopt) ? *mls_projection : smoothed_position; +#else + if(m_context->m_triangles_aabb_tree.squared_distance(smoothed_position) < m_context->m_aabb_epsilon) { - new_pos = m_triangles_aabb_tree.closest_point(smoothed_position); + new_pos = m_context->m_triangles_aabb_tree.closest_point(smoothed_position); } else { using Ray = typename Tr::Geom_traits::Ray_3; - using Projection = std::optional< - typename AABB_triangle_tree::template Intersection_and_primitive_id::Type>; - - auto get_intersection_point = - [](const Projection& proj) -> std::optional - { - const auto intersection = proj.value().first; - if (const Point_3* p = std::get_if(&intersection)) - return *p; - else - return std::nullopt; - }; - - // this lambda is called only when we are sure that proj is a Segment - auto get_intersection_midpoint = - [](const Projection& proj) -> std::optional - { - const auto intersection = proj.value().first; - using Segment = typename Tr::Geom_traits::Segment_3; - if (const Segment* s = std::get_if(&intersection)) - return CGAL::midpoint(s->source(), s->target()); - else - { - CGAL_assertion(false); - return std::nullopt; - } - }; - - const auto n = vertices_normals.at(v).at(si); + using Projection = + std::optional::Type>; + + auto get_intersection_point = [](const Projection& proj) -> std::optional { + const auto intersection = proj.value().first; + if(const Point_3* p = std::get_if(&intersection)) + return *p; + return std::nullopt; + }; + + auto get_intersection_midpoint = [](const Projection& proj) -> std::optional { + const auto intersection = proj.value().first; + using Segment = typename Tr::Geom_traits::Segment_3; + if(const Segment* s = std::get_if(&intersection)) + return CGAL::midpoint(s->source(), s->target()); + CGAL_assertion(false); + return std::nullopt; + }; + + const auto n = m_context->m_vertices_normals.at(v).at(si); const Ray ray = tr.geom_traits().construct_ray_3_object()(current_pos, n); - - const Projection proj = m_triangles_aabb_tree.first_intersection(ray); - const Projection proj_opp = m_triangles_aabb_tree.first_intersection( - tr.geom_traits().construct_opposite_ray_3_object()(ray)); + const Projection proj = m_context->m_triangles_aabb_tree.first_intersection(ray); + const Projection proj_opp = m_context->m_triangles_aabb_tree.first_intersection( + tr.geom_traits().construct_opposite_ray_3_object()(ray)); if(proj != std::nullopt && proj_opp == std::nullopt) { const auto p = get_intersection_point(proj); - if (p != std::nullopt) - new_pos = p.value(); - else - new_pos = get_intersection_midpoint(proj).value(); + new_pos = (p != std::nullopt) ? p.value() : get_intersection_midpoint(proj).value(); } else if(proj == std::nullopt && proj_opp != std::nullopt) { const auto p = get_intersection_point(proj_opp); - if (p != std::nullopt) - new_pos = p.value(); - else - new_pos = get_intersection_midpoint(proj_opp).value(); + new_pos = (p != std::nullopt) + ? p.value() + : get_intersection_midpoint(proj_opp).value(); } else if(proj != std::nullopt && proj_opp != std::nullopt) { @@ -930,229 +1006,144 @@ std::size_t smooth_vertices_on_surfaces(C3t3& c3t3, new_pos = smoothed_position; } else //no valid projection + { new_pos = smoothed_position; + } } #endif //CGAL_TET_REMESHING_SMOOTHING_WITH_MLS - if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){ - nb_done_2d++; - } -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - os_surf << "2 " << current_pos << " " << new_pos << std::endl; -#endif + const auto& inc_cells = m_context->m_inc_cells[vid]; + result = BaseClass::check_inversion_and_move(v, new_pos, inc_cells, tr, m_context->m_total_move); } else if (nb_neighbors > 0) { #ifdef CGAL_TET_REMESHING_SMOOTHING_WITH_MLS std::optional mls_proj = project(si, current_pos); - if (mls_proj == std::nullopt) - continue; + if(mls_proj == std::nullopt) + return false; - const Point_3 new_pos = *mls_proj; + new_pos = *mls_proj; #else // AABB_tree projection - const Point_3 new_pos = m_segments_aabb_tree.closest_point(current_pos); -#endif // CGAL_TET_REMESHING_SMOOTHING_WITH_MLS + new_pos = m_context->m_segments_aabb_tree.closest_point(current_pos); +#endif //CGAL_TET_REMESHING_SMOOTHING_WITH_MLS - if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){ - nb_done_2d++; - } -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - os_surf0 << "2 " << current_pos << " " << new_pos << std::endl; -#endif + const auto& inc_cells = m_context->m_inc_cells[vid]; + result = BaseClass::check_inversion_and_move(v, new_pos, inc_cells, tr, m_context->m_total_move); } + + return result; } - return nb_done_2d; -} + std::string operation_name() const override { return "Vertex Smooth (Surface Vertices)"; } +}; -template -std::size_t smooth_internal_vertices(C3t3& c3t3, - const IncidentCells& inc_cells, - FT& total_move -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - , std::ofstream& os_vol -#endif - ) +template +class Internal_vertex_smooth_operation + : public Vertex_smooth_operation_base, + public Elementary_operation { - std::size_t nb_done_3d = 0; - auto& tr = c3t3.triangulation(); +public: + using BaseClass = Vertex_smooth_operation_base; + using Base_operation = Elementary_operation; + using Element_type = typename Base_operation::Element_type; + static_assert(std::is_same_v, + "Element_type should be Vertex_handle"); + using Element_range = typename Base_operation::Element_range; + + using BaseClass::m_context; + + using typename BaseClass::Cell_handle; + using typename BaseClass::Edge; + using typename BaseClass::FT; + using typename BaseClass::Point_3; + using typename BaseClass::Tr; + using typename BaseClass::Vector_3; + using typename BaseClass::Vertex_handle; - const std::size_t nbv = tr.number_of_vertices(); - const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/}; - std::vector moves(nbv, default_move); - /*for dim 3 vertices, start counting neighbors directly from 0*/ +public: + Internal_vertex_smooth_operation(std::shared_ptr context) + : BaseClass(context) {} - for (const Edge& e : tr.finite_edges()) + void perform_global_preprocessing(const C3t3& c3t3) const { - if (is_outside(e, c3t3, m_cell_selector)) - continue; - else + auto& tr = c3t3.triangulation(); + auto& moves = m_context->m_moves; + + using Move = typename BaseClass::Context::Move; + const std::size_t nbv = tr.number_of_vertices(); + const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/}; + moves.assign(nbv, default_move); + /*for dim 3 vertices, start counting neighbors directly from 0*/ + + for (const Edge& e : tr.finite_edges()) { - const auto [vh0, vh1] = make_vertex_pair(e); + if (is_outside(e, c3t3, m_context->m_cell_selector)) + continue; + else + { + const auto [vh0, vh1] = make_vertex_pair(e); - const std::size_t& i0 = vertex_id(vh0); - const std::size_t& i1 = vertex_id(vh1); + const std::size_t& i0 = m_context->vertex_id(vh0); + const std::size_t& i1 = m_context->vertex_id(vh1); - const bool vh0_moving = (c3t3.in_dimension(vh0) == 3 && is_free(i0)); - const bool vh1_moving = (c3t3.in_dimension(vh1) == 3 && is_free(i1)); + const bool vh0_moving = (c3t3.in_dimension(vh0) == 3 && m_context->is_free(i0)); + const bool vh1_moving = (c3t3.in_dimension(vh1) == 3 && m_context->is_free(i1)); - if (!vh0_moving && !vh1_moving) - continue; + if (!vh0_moving && !vh1_moving) + continue; - const Point_3& p0 = point(vh0->point()); - const Point_3& p1 = point(vh1->point()); - const FT density = density_along_segment(e, c3t3); + const Point_3& p0 = point(vh0->point()); + const Point_3& p1 = point(vh1->point()); + const FT density = BaseClass::density_along_segment(e, c3t3); - if (vh0_moving) - { - moves[i0].move += density * Vector_3(p0, p1); - moves[i0].mass += density; - ++moves[i0].neighbors; - } - if (vh1_moving) - { - moves[i1].move += density * Vector_3(p1, p0); - moves[i1].mass += density; - ++moves[i1].neighbors; + if (vh0_moving) + { + moves[i0].move += density * Vector_3(p0, p1); + moves[i0].mass += density; + ++moves[i0].neighbors; + } + if (vh1_moving) + { + moves[i1].move += density * Vector_3(p1, p0); + moves[i1].mass += density; + ++moves[i1].neighbors; + } } } } - // iterate over vertices and move - for(Vertex_handle v : tr.finite_vertex_handles()) + Element_range get_elements(const C3t3& c3t3) const override { - const std::size_t vid = vertex_id(v); - if (!is_free(vid)) - continue; - - if (c3t3.in_dimension(v) == 3 && moves[vid].neighbors > 1) - { -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - os_vol << "2 " << point(v->point()); -#endif - const Vector_3 move = moves[vid].move / moves[vid].mass;// static_cast(neighbors[vid]); - Point_3 new_pos = point(v->point()) + move; - if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){ - nb_done_3d++; - } - -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - os_vol << " " << point(v->point()) << std::endl; -#endif - } + perform_global_preprocessing(c3t3); + return c3t3.triangulation().finite_vertex_handles(); } - return nb_done_3d; -} -public: - void smooth_vertices(C3t3& c3t3) + bool execute_operation(const Element_type& v, C3t3& c3t3) override { - typedef typename C3t3::Cell_handle Cell_handle; - -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - std::ofstream os_surf("smooth_surfaces.polylines.txt"); - std::ofstream os_surf0("smooth_surfaces0.polylines.txt"); - std::ofstream os_vol("smooth_volume.polylines.txt"); -#endif - -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - std::cout << "Smooth vertices..."; - std::cout.flush(); - - std::size_t nb_done_3d = 0; - std::size_t nb_done_2d = 0; - std::size_t nb_done_1d = 0; - CGAL::Real_timer timer; - timer.start(); -#endif - - FT total_move = 0.; - - Tr& tr = c3t3.triangulation(); - - //collect a map of vertices surface indices - std::unordered_map > vertices_surface_indices; - if(!m_protect_boundaries) - collect_vertices_surface_indices(c3t3, vertices_surface_indices); - - //collect a map of normals at surface vertices - std::unordered_map>> vertices_normals; - if(!m_protect_boundaries) - compute_vertices_normals(c3t3, vertices_normals); - - //collect ids - reset_vertex_id_map(tr); - - //are vertices free to move? indices are in `vertex_id` - reset_free_vertices(tr); - - //collect incident cells - using Incident_cells_vector = boost::container::small_vector; - const std::size_t nbv = tr.number_of_vertices(); - std::vector inc_cells(nbv, Incident_cells_vector{}); - collect_incident_cells(tr, inc_cells); + auto& tr = c3t3.triangulation(); + auto& moves = m_context->m_moves; - if (!m_protect_boundaries && m_smooth_constrained_edges) - { -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - nb_done_1d = -#endif - smooth_edges_in_complex(c3t3, - vertices_surface_indices, inc_cells, vertices_normals, total_move -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - , os_surf -#endif - ); - } + const std::size_t vid = m_context->vertex_id(v); + if (!m_context->is_free(vid)) + return false; - /////////////// EDGES ON SURFACE, BUT NOT IN COMPLEX ////////////////// - if (!m_protect_boundaries) + if (c3t3.in_dimension(v) == 3 && moves[vid].neighbors > 1) { -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - nb_done_2d = -#endif - smooth_vertices_on_surfaces(c3t3, - vertices_surface_indices, inc_cells, vertices_normals, - total_move -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - , os_surf, os_surf0 -#endif - ); + const Vector_3 move = moves[vid].move / moves[vid].mass; + const Point_3 new_pos = point(v->point()) + move; + return BaseClass::check_inversion_and_move(v, new_pos, m_context->incident_cells(vid), tr, + m_context->m_total_move); } - CGAL_expensive_assertion(CGAL::Tetrahedral_remeshing::debug::are_cell_orientations_valid(tr)); - //// end if(!protect_boundaries) - - ////////////// INTERNAL VERTICES /////////////////////// -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - nb_done_3d = -#endif - smooth_internal_vertices(c3t3, inc_cells, - total_move -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - , os_vol -#endif - ); - - CGAL_expensive_assertion(CGAL::Tetrahedral_remeshing::debug::are_cell_orientations_valid(tr)); - -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - timer.stop(); - std::size_t nb_done = nb_done_3d + nb_done_2d + nb_done_1d; - std::cout << " done (" - << nb_done_1d << "/" << nb_done_2d << "/" << nb_done_3d << " vertices smoothed," - << " average move = " << (total_move / nb_done) - << ", in "<< timer.time() << " seconds)." << std::endl; -#endif -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension( - c3t3.triangulation(), "c3t3_vertices_after_smoothing"); - os_surf.close(); - os_vol.close(); -#endif + return false; } -};//end class Tetrahedral_remeshing_smoother + std::string operation_name() const override { return "Vertex Smooth (Internal Vertices)"; } +}; + }//namespace internal }//namespace Tetrahedral_adaptive_remeshing }//namespace CGAL diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h index 8a53ed2ea9e..af56cc11df3 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h @@ -15,12 +15,10 @@ #include -#include -#include -#include #include #include +#include #include #include @@ -28,10 +26,6 @@ #include #include -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE -#include -#endif - namespace CGAL { namespace Tetrahedral_remeshing @@ -340,127 +334,147 @@ auto can_be_split(const typename C3T3::Edge& e, -template -void split_long_edges(C3T3& c3t3, - const Sizing& sizing, - const bool protect_boundaries, - CellSelector cell_selector, - Visitor& visitor) +class Edge_split_operation + : public Elementary_operation, + std::vector>> { - typedef typename C3T3::Triangulation T3; - typedef typename T3::Cell_handle Cell_handle; - typedef typename T3::Edge Edge; - typedef typename T3::Vertex_handle Vertex_handle; - typedef typename std::pair Edge_vv; - - typedef typename T3::Geom_traits::FT FT; - typedef boost::bimap< - boost::bimaps::set_of, - boost::bimaps::multiset_of > > Boost_bimap; - typedef typename Boost_bimap::value_type long_edge; - -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - std::cout << "Split long edges..."; - std::cout.flush(); - std::size_t nb_splits = 0; - CGAL::Real_timer timer; - timer.start(); +public: + using Tr = typename C3t3::Triangulation; + using Vertex_handle = typename Tr::Vertex_handle; + using Cell_handle = typename Tr::Cell_handle; + using Edge = typename Tr::Edge; + using Edge_vv = std::pair; + using FT = typename Tr::Geom_traits::FT; + + // Candidates are stored as vertex pairs, captured at collection time. A raw + // Edge (Cell_handle, i, j) would go stale: each split destroys and recycles + // cells, so by the time the executor reaches a later candidate its Cell_handle + // may point at a different cell. Vertices are never removed by a split, so the + // vertex pair stays valid and is re-resolved to the current edge via is_edge(). + using Long_edges = std::vector; + using Base_operation = Elementary_operation; + using Element_type = typename Base_operation::Element_type; + static_assert(std::is_same_v, "Element_type must be Edge_vv"); + using ElementSource = typename Base_operation::Element_range; + +private: + const SizingFunction& m_sizing; + const CellSelector& m_cell_selector; + bool m_protect_boundaries; + Visitor& m_visitor; + +#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG + mutable std::ofstream m_can_be_split_ofs; + mutable std::ofstream m_split_failed_ofs; + mutable std::ofstream m_midpoints_ofs; #endif - //collect long edges - T3& tr = c3t3.triangulation(); - Boost_bimap long_edges; - for (Edge e : tr.finite_edges()) +public: + Edge_split_operation(const SizingFunction& sizing, + const CellSelector& cell_selector, + const bool protect_boundaries, + Visitor& visitor) + : m_sizing(sizing) + , m_cell_selector(cell_selector) + , m_protect_boundaries(protect_boundaries) + , m_visitor(visitor) {} + + ElementSource get_elements(const C3t3& c3t3) const override { - auto [splittable, boundary] = can_be_split(e, c3t3, protect_boundaries, cell_selector); - if (!splittable) - continue; + struct Long_edge_with_length + { + Edge edge; + FT sqlength; + }; + std::vector long_edges_with_lengths; + const Tr& tr = c3t3.triangulation(); - const std::optional sqlen = is_too_long(e, boundary, sizing, c3t3, cell_selector); - if(sqlen != std::nullopt) - long_edges.insert(long_edge(make_vertex_pair(e), sqlen.value())); - } + for (Edge e : tr.finite_edges()) + { + auto [splittable, boundary] = can_be_split(e, c3t3, m_protect_boundaries, m_cell_selector); + if (!splittable) + continue; -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - debug::dump_edges(long_edges, "long_edges.polylines.txt"); + const std::optional sqlen = is_too_long(e, boundary, m_sizing, c3t3, m_cell_selector); + if (sqlen != std::nullopt) + long_edges_with_lengths.push_back(Long_edge_with_length{e, sqlen.value()}); + } - std::ofstream can_be_split_ofs("can_be_split_edges.polylines.txt"); - std::ofstream split_failed_ofs("split_failed.polylines.txt"); + // longest first; stable to match the original bimap's ordering + std::stable_sort(long_edges_with_lengths.begin(), long_edges_with_lengths.end(), + [](const Long_edge_with_length& a, const Long_edge_with_length& b) { + return a.sqlength > b.sqlength; + }); - std::ofstream ofs("midpoints.off"); - ofs << "OFF" << std::endl; - ofs << long_edges.size() << " 0 0" << std::endl; +#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG + { + std::ofstream ofs("long_edges.polylines.txt"); + for (const auto& le : long_edges_with_lengths) + ofs << "2 " << point(le.edge.first->point()) + << " " << point(le.edge.second->point()) << std::endl; + } + m_can_be_split_ofs.open("can_be_split_edges.polylines.txt"); + m_split_failed_ofs.open("split_failed.polylines.txt"); + m_midpoints_ofs.open("midpoints.off"); + m_midpoints_ofs << "OFF" << std::endl; + m_midpoints_ofs << long_edges_with_lengths.size() << " 0 0" << std::endl; #endif - while(!long_edges.empty()) + + Long_edges long_edges; + long_edges.reserve(long_edges_with_lengths.size()); + for(const auto& ef : long_edges_with_lengths) + long_edges.push_back(make_vertex_pair(ef.edge)); + return long_edges; + } + + bool execute_operation(const Element_type& element, C3t3& c3t3) override { - //the edge with longest length - typename Boost_bimap::right_map::iterator eit = long_edges.right.begin(); - Edge_vv e = eit->second; -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS - const double sqlen = eit->first; -#endif - long_edges.right.erase(eit); + Tr& tr = c3t3.triangulation(); + const Edge_vv& e = element; Cell_handle cell; int i1, i2; - if ( tr.tds().is_edge(e.first, e.second, cell, i1, i2)) - { - Edge edge(cell, i1, i2); + if (!tr.tds().is_edge(e.first, e.second, cell, i1, i2)) + return false; - //check that splittability has not changed - auto [splittable, _] = can_be_split(edge, c3t3, protect_boundaries, cell_selector); - if (!splittable) - continue; + Edge edge(cell, i1, i2); + + // check that splittability has not changed + if (!can_be_split(edge, c3t3, m_protect_boundaries, m_cell_selector).can_be_split) + return false; #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - else - can_be_split_ofs << "2 " << edge.first->vertex(edge.second)->point() - << " " << edge.first->vertex(edge.third)->point() << std::endl; + m_can_be_split_ofs << "2 " << edge.first->vertex(edge.second)->point() + << " " << edge.first->vertex(edge.third)->point() << std::endl; #endif - visitor.before_split(tr, edge); - Vertex_handle vh = split_edge(edge, cell_selector, c3t3); - if(vh != Vertex_handle()) - visitor.after_split(tr, vh); + m_visitor.before_split(tr, edge); + Vertex_handle vh = split_edge(edge, m_cell_selector, c3t3); + if (vh == Vertex_handle()) + { #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - else - split_failed_ofs << "2 " << edge.first->vertex(edge.second)->point() << " " + m_split_failed_ofs << "2 " << edge.first->vertex(edge.second)->point() << " " << edge.first->vertex(edge.third)->point() << std::endl; - if (vh != Vertex_handle()) - ofs << vh->point() << std::endl; -#endif - -#if defined(CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS) \ -|| defined(CGAL_TETRAHEDRAL_REMESHING_VERBOSE) - if (vh != Vertex_handle()) - ++nb_splits; -#endif - -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS - std::cout << "\rSplit... (" - << long_edges.left.size() << " long edges, " - << "length = " << std::sqrt(sqlen) << ", " - << nb_splits << " splits)"; - std::cout.flush(); #endif + return false; } - }//end loop on long_edges #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - if(can_be_split_ofs.is_open()) can_be_split_ofs.close(); - if(split_failed_ofs.is_open()) split_failed_ofs.close(); - if(ofs.is_open()) ofs.close(); + m_midpoints_ofs << vh->point() << std::endl; #endif + m_visitor.after_split(tr, vh); + return true; + } -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - timer.stop(); - std::cout << " done (" << nb_splits << " splits, in " - << timer.time() << " sec)." << std::endl; -#endif -} + std::string operation_name() const override { return "Split long edges"; } +}; } // internal } // Tetrahedral_remeshing diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h index 187178c984a..a6e1b00b119 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -86,7 +87,7 @@ class Adaptive_remesher typedef typename C3t3::Curve_index Curve_index; typedef typename C3t3::Corner_index Corner_index; - typedef Tetrahedral_remeshing_smoother Smoother; + typedef Vertex_smoothing_context SmoothingContext; private: C3t3 m_c3t3; @@ -94,7 +95,7 @@ class Adaptive_remesher const bool m_protect_boundaries; CellSelector m_cell_selector; Visitor& m_visitor; - Smoother m_vertex_smoother;//initialized with initial surface + std::shared_ptr m_smoothing_context;//built from the initial surface C3t3* m_c3t3_pbackup; std::vector m_far_points; @@ -116,14 +117,14 @@ class Adaptive_remesher , m_protect_boundaries(protect_boundaries) , m_cell_selector(cell_selector) , m_visitor(visitor) - , m_vertex_smoother(sizing, cell_selector, protect_boundaries, smooth_constrained_edges) , m_c3t3_pbackup(NULL) , m_tr_pbackup(&tr) { m_c3t3.triangulation().swap(tr); init_c3t3(vcmap, ecmap, fcmap); - m_vertex_smoother.init(m_c3t3); + m_smoothing_context = std::make_shared( + m_c3t3, m_sizing, m_cell_selector, m_protect_boundaries, smooth_constrained_edges); #ifdef CGAL_DUMP_REMESHING_STEPS CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "00-init"); @@ -147,14 +148,14 @@ class Adaptive_remesher , m_protect_boundaries(protect_boundaries) , m_cell_selector(cell_selector) , m_visitor(visitor) - , m_vertex_smoother(sizing, cell_selector, protect_boundaries, smooth_constrained_edges) , m_c3t3_pbackup(&c3t3) , m_tr_pbackup(NULL) { m_c3t3.swap(c3t3); init_c3t3(vcmap, ecmap, fcmap); - m_vertex_smoother.init(m_c3t3); + m_smoothing_context = std::make_shared( + m_c3t3, m_sizing, m_cell_selector, m_protect_boundaries, smooth_constrained_edges); #ifdef CGAL_DUMP_REMESHING_STEPS CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "00-init"); @@ -171,8 +172,10 @@ class Adaptive_remesher void split() { CGAL_assertion(check_vertex_dimensions()); - split_long_edges(m_c3t3, m_sizing, m_protect_boundaries, - m_cell_selector, m_visitor); + typedef Edge_split_operation EdgeSplitOp; + EdgeSplitOp split_op(m_sizing, m_cell_selector, m_protect_boundaries, m_visitor); + Elementary_operation_execution_sequential executor; + executor.execute(split_op, m_c3t3); #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL_assertion(tr().tds().is_valid(true)); @@ -195,8 +198,10 @@ class Adaptive_remesher void collapse() { CGAL_assertion(check_vertex_dimensions()); - collapse_short_edges(m_c3t3, m_sizing, m_protect_boundaries, - m_cell_selector, m_visitor); + typedef Edge_collapse_operation EdgeCollapseOp; + EdgeCollapseOp collapse_op(m_sizing, m_cell_selector, m_protect_boundaries, m_visitor); + Elementary_operation_execution_sequential executor; + executor.execute(collapse_op, m_c3t3); #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL_assertion(tr().tds().is_valid(true)); @@ -216,8 +221,22 @@ class Adaptive_remesher void flip() { - flip_edges(m_c3t3, m_protect_boundaries, - m_cell_selector, m_visitor); + typedef Internal_edge_flip_operation InternalFlipOp; + typedef Boundary_edge_flip_operation BoundaryFlipOp; + + // one incident-cells cache shared by both passes, as in the former flip_edges() + typename InternalFlipOp::Incident_cells_map inc_cells; + + InternalFlipOp internal_flip_op(m_cell_selector, m_visitor, inc_cells); + Elementary_operation_execution_sequential internal_executor; + internal_executor.execute(internal_flip_op, m_c3t3); + + if (!m_protect_boundaries) + { + BoundaryFlipOp boundary_flip_op(m_cell_selector, m_visitor, inc_cells); + Elementary_operation_execution_sequential boundary_executor; + boundary_executor.execute(boundary_flip_op, m_c3t3); + } #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL_assertion(tr().tds().is_valid(true)); @@ -237,7 +256,29 @@ class Adaptive_remesher void smooth() { - m_vertex_smoother.smooth_vertices(m_c3t3); + m_smoothing_context->refresh(m_c3t3); + + // Order matches the former Tetrahedral_remeshing_smoother::smooth_vertices(): + // complex (1D) edges, then surface (2D) vertices, then internal (3D) vertices. + if (!m_protect_boundaries) + { + if (m_smoothing_context->m_smooth_constrained_edges) + { + Complex_edge_vertex_smooth_operation op(m_smoothing_context); + Elementary_operation_execution_sequential executor; + executor.execute(op, m_c3t3); + } + { + Surface_vertex_smooth_operation op(m_smoothing_context); + Elementary_operation_execution_sequential executor; + executor.execute(op, m_c3t3); + } + } + { + Internal_vertex_smooth_operation op(m_smoothing_context); + Elementary_operation_execution_sequential executor; + executor.execute(op, m_c3t3); + } #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL_assertion(tr().tds().is_valid(true)); @@ -633,7 +674,7 @@ class Adaptive_remesher #endif } - m_vertex_smoother.start_flip_smooth_steps(m_c3t3); + m_smoothing_context->start_flip_smooth_steps(m_c3t3); while (it_nb < max_it + nb_extra_iterations) { ++it_nb; 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..66cc5e84fc4 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 @@ -1906,7 +1906,6 @@ void get_edge_info(const typename C3t3::Edge& edge, } } - template void remove_from_bimap(const typename EdgesBimap::left_map::key_type& e, EdgesBimap& edges)