From 81751462f876cce85cd380191a6ef9ccbdc8951f Mon Sep 17 00:00:00 2001 From: IasonManolas Date: Mon, 6 Jul 2026 18:56:54 +0300 Subject: [PATCH 1/5] Tetrahedral_remeshing: add sequential elementary-operation framework Introduces ElementaryOperation, the interface for a single mesh-editing operation (candidate source, per-candidate execution, name), and ElementaryOperationExecutionSequential, which collects an operation's candidates once and applies each in the order the operation returns them. This is the common scaffolding that the split, collapse, flip and smooth operations are refactored onto in the following commits. It is generic and not yet used by any code path. --- .../internal/elementary_operations.h | 84 +++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/elementary_operations.h diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/elementary_operations.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/elementary_operations.h new file mode 100644 index 00000000000..f3d84e5a452 --- /dev/null +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/elementary_operations.h @@ -0,0 +1,84 @@ +// 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 +#endif + +namespace CGAL { +namespace Tetrahedral_remeshing { +namespace internal { + +template +class ElementaryOperation +{ +public: + using C3t3 = C3t3_; + using Triangulation = typename C3t3::Triangulation; + using ElementType = ElementType_; + using ElementSource = ElementSource_; + + ElementaryOperation() = default; + virtual ~ElementaryOperation() = default; + + virtual ElementSource get_element_source(const C3t3& c3t3) const = 0; + virtual bool execute_operation(const ElementType& e, C3t3& c3t3) = 0; + virtual std::string operation_name() const = 0; +}; + +template +class ElementaryOperationExecutionSequential +{ +public: + using C3t3 = typename Operation::C3t3; + using ElementSource = typename Operation::ElementSource; + + bool execute(Operation& op, C3t3& c3t3) const + { + ElementSource candidates = op.get_element_source(c3t3); + if (candidates.empty()) + return false; + +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + std::size_t nb_done = 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 + std::cout << op.operation_name() << ": " << nb_done << "/" + << candidates.size() << " done." << std::endl; +#endif + return true; + } +}; + +} // namespace internal +} // namespace Tetrahedral_remeshing +} // namespace CGAL + +#endif // CGAL_TETRAHEDRAL_REMESHING_ELEMENTARY_OPERATIONS_H From a2e83867c2ea0f37f783e1710047312d8d6c241c Mon Sep 17 00:00:00 2001 From: IasonManolas Date: Tue, 7 Jul 2026 09:39:41 +0300 Subject: [PATCH 2/5] Tetrahedral_remeshing: refactor edge split onto the elementary-operation framework Replace the split_long_edges() free function with EdgeSplitOperation, an ElementaryOperation whose get_element_source() collects and length- orders the splittable edges (stable_sort, longest first, matching the original bimap ordering) and whose execute_operation() re-validates and splits one edge. Adaptive_remesher::split() now drives it through ElementaryOperationExecutionSequential. Output is unchanged: remeshing a mesh through this path yields a byte-identical result to the previous implementation. Add generic timing/progress reporting to the sequential executor: elapsed time per operation under CGAL_TETRAHEDRAL_REMESHING_VERBOSE and a live per-candidate progress line under CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS, replacing the equivalent per-operation reporting the free function used to do itself. The visitor before_split()/after_split() hooks and the CGAL_TETRAHEDRAL_REMESHING_DEBUG diagnostic dumps are preserved. --- .../internal/elementary_operations.h | 14 +- .../internal/split_long_edges.h | 200 +++++++++--------- .../tetrahedral_adaptive_remeshing_impl.h | 7 +- 3 files changed, 116 insertions(+), 105 deletions(-) diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/elementary_operations.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/elementary_operations.h index f3d84e5a452..a6ca50b3298 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/elementary_operations.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/elementary_operations.h @@ -18,6 +18,7 @@ #include #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE +#include #include #include #endif @@ -58,6 +59,11 @@ class ElementaryOperationExecutionSequential #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE std::size_t nb_done = 0; + 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) { @@ -67,11 +73,17 @@ class ElementaryOperationExecutionSequential ++nb_done; #endif } +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS + std::cout << "\r" << op.operation_name() << "... (" + << ++nb_processed << "/" << candidates.size() << ")"; + std::cout.flush(); +#endif } #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + timer.stop(); std::cout << op.operation_name() << ": " << nb_done << "/" - << candidates.size() << " done." << std::endl; + << candidates.size() << " done (" << timer.time() << " sec)." << std::endl; #endif return true; } 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..6119762eceb 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,129 @@ 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) +template +class EdgeSplitOperation + : public ElementaryOperation>, + 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; + + using Long_edges_with_lengths = std::vector>; + using Base = ElementaryOperation, Long_edges_with_lengths>; + using ElementType = typename Base::ElementType; + using ElementSource = typename Base::ElementSource; + +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: + EdgeSplitOperation(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_element_source(const C3t3& c3t3) const override { - auto [splittable, boundary] = can_be_split(e, c3t3, protect_boundaries, cell_selector); - if (!splittable) - continue; + Long_edges_with_lengths long_edges; + 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.push_back(std::make_pair(*sqlen, make_vertex_pair(e))); + } - 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.begin(), long_edges.end(), + [](const std::pair& a, const std::pair& b) { + return a.first > b.first; + }); - 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) + ofs << "2 " << point(le.second.first->point()) + << " " << point(le.second.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.size() << " 0 0" << std::endl; #endif - while(!long_edges.empty()) + return long_edges; + } + + bool execute_operation(const ElementType& 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.second; 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..cece5ef993f 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 @@ -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 EdgeSplitOperation EdgeSplitOp; + EdgeSplitOp split_op(m_sizing, m_cell_selector, m_protect_boundaries, m_visitor); + ElementaryOperationExecutionSequential executor; + executor.execute(split_op, m_c3t3); #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL_assertion(tr().tds().is_valid(true)); From eac2b14e083ff43d30dd5d2da5ee0e8215d77ab5 Mon Sep 17 00:00:00 2001 From: IasonManolas Date: Tue, 7 Jul 2026 10:38:41 +0300 Subject: [PATCH 3/5] Tetrahedral_remeshing: refactor edge flip onto the elementary-operation framework Replace the flip_edges() free function with two ElementaryOperations sharing EdgeFlipOperationBase: - InternalEdgeFlipOperation mirrors flip_all_edges(): reset the cell caches, collect the internal edges, and run find_best_flip() on each. - BoundaryEdgeFlipOperation mirrors flipBoundaryEdges(): compute the per-vertex boundary valences once in get_element_source(), then flip each boundary edge on the surface when it lowers the valence cost, updating the valences. Adaptive_remesher::flip() owns a single incident-cells cache and passes it to both operations by reference, running the internal pass, then the boundary pass when boundaries are not protected -- reproducing the map sharing and pass order of the former flip_edges(). The find_best_flip()/flip_on_surface() machinery and the CGAL_TETRAHEDRAL_REMESHING_DEBUG orientation check are unchanged. Remeshing a mesh through this path yields a byte-identical result to the previous implementation. --- .../internal/flip_edges.h | 299 +++++++++++++----- .../tetrahedral_adaptive_remeshing_impl.h | 18 +- 2 files changed, 238 insertions(+), 79 deletions(-) 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 ed22929e608..8629abc5434 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 @@ -1899,97 +1896,245 @@ std::size_t flipBoundaryEdges( return nb_success; } -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 EdgeFlipOperationBase { - 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: + typedef typename C3t3::Triangulation Tr; + 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 boost::container::small_vector Cells_vector; + typedef std::unordered_map Incident_cells_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; + + EdgeFlipOperationBase(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 InternalEdgeFlipOperation + : public EdgeFlipOperationBase, + public ElementaryOperation, + std::vector>> +{ + using BaseClass = EdgeFlipOperationBase; + 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 = ElementaryOperation>; + using ElementType = typename Base::ElementType; + using ElementSource = typename Base::ElementSource; + + InternalEdgeFlipOperation(CellSelector& cell_selector, Visitor& visitor, Incident_cells_map& incident_cells) + : BaseClass(cell_selector, visitor, incident_cells) {} + + ElementSource get_element_source(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 ElementType& 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 BoundaryEdgeFlipOperation + : public EdgeFlipOperationBase, + public ElementaryOperation, + std::vector>> +{ + using BaseClass = EdgeFlipOperationBase; + 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 = ElementaryOperation>; + using ElementType = typename Base::ElementType; + using ElementSource = typename Base::ElementSource; + + BoundaryEdgeFlipOperation(CellSelector& cell_selector, Visitor& visitor, Incident_cells_map& incident_cells) + : BaseClass(cell_selector, visitor, incident_cells) {} + + ElementSource get_element_source(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 ElementType& 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; + + const Facet& f0 = boundary_facets[0]; + const Facet& f1 = boundary_facets[1]; + + // 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); + + if (!tr.tds().is_edge(vh2, vh3)) // most-likely to happen early exit + { + const Surface_patch_index surfi = c3t3.surface_patch_index(boundary_facets[0]); + + int v0 = m_boundary_vertices_valences.at(vh0)[surfi]; + int v1 = m_boundary_vertices_valences.at(vh1)[surfi]; + int v2 = m_boundary_vertices_valences.at(vh2)[surfi]; + int v3 = m_boundary_vertices_valences.at(vh3)[surfi]; + + if (v0 < 2 || v1 < 2 || v2 < 2 || v3 < 2) + return false; + + int m0 = (m_boundary_vertices_valences.at(vh0).size() > 1 ? 4 : 6); + int m1 = (m_boundary_vertices_valences.at(vh1).size() > 1 ? 4 : 6); + int m2 = (m_boundary_vertices_valences.at(vh2).size() > 1 ? 4 : 6); + int m3 = (m_boundary_vertices_valences.at(vh3).size() > 1 ? 4 : 6); + + int initial_cost = (v0 - m0)*(v0 - m0) + + (v1 - m1)*(v1 - m1) + + (v2 - m2)*(v2 - m2) + + (v3 - m3)*(v3 - m3); + v0--; + v1--; + v2++; + v3++; + + int final_cost = (v0 - m0)*(v0 - m0) + + (v1 - m1)*(v1 - m1) + + (v2 - m2)*(v2 - m2) + + (v3 - m3)*(v3 - m3); + if (initial_cost > final_cost) + { + const Sliver_removal_result db = flip_on_surface(c3t3, edge, vh2, vh3, + inc_cells, + MIN_ANGLE_BASED, + m_visitor); + if (db == VALID_FLIP) + { + Cell_handle c_new; + int li, lj, lk; + tr.tds().is_facet(vh2, vh3, vh0, c_new, li, lj, lk); + c3t3.add_to_complex(c_new, (6 - li - lj - lk), surfi); + + tr.tds().is_facet(vh2, vh3, vh1, c_new, li, lj, lk); + c3t3.add_to_complex(c_new, (6 - li - lj - lk), surfi); + + m_boundary_vertices_valences[vh0][surfi]--; + m_boundary_vertices_valences[vh1][surfi]--; + m_boundary_vertices_valences[vh2][surfi]++; + m_boundary_vertices_valences[vh3][surfi]++; + + return true; + } + } + } + + return false; + } + + 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/tetrahedral_adaptive_remeshing_impl.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h index cece5ef993f..896e3c1e982 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 @@ -219,8 +219,22 @@ class Adaptive_remesher void flip() { - flip_edges(m_c3t3, m_protect_boundaries, - m_cell_selector, m_visitor); + typedef InternalEdgeFlipOperation InternalFlipOp; + typedef BoundaryEdgeFlipOperation 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); + ElementaryOperationExecutionSequential 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); + ElementaryOperationExecutionSequential boundary_executor; + boundary_executor.execute(boundary_flip_op, m_c3t3); + } #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL_assertion(tr().tds().is_valid(true)); From 538324201b82c25e4becb65c6dfcdcd9296ec7e3 Mon Sep 17 00:00:00 2001 From: IasonManolas Date: Tue, 7 Jul 2026 18:04:10 +0300 Subject: [PATCH 4/5] Tetrahedral_remeshing: refactor vertex smoothing onto the elementary-operation framework Replace Tetrahedral_remeshing_smoother with Vertex_smoothing_context (smoothing state + setup helpers), VertexSmoothOperationBase (per-vertex helpers), and one ElementaryOperation per pass: ComplexEdgeVertexSmoothOperation, SurfaceVertexSmoothOperation and InternalVertexSmoothOperation. Each get_element_source() accumulates the per-vertex moves and each execute_operation() moves one vertex, mirroring the former smooth_edges_in_complex/on_surfaces/internal_vertices. Adaptive_remesher keeps the context as a persistent member and runs the three operations in the previous order (1D, 2D, 3D) each smoothing step. Remeshing yields a byte-identical result to the previous implementation. --- .../internal/smooth_vertices.h | 1089 ++++++++--------- .../tetrahedral_adaptive_remeshing_impl.h | 38 +- 2 files changed, 544 insertions(+), 583 deletions(-) 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..b96cadc6fc3 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,26 +42,36 @@ 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; +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; - 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; + struct Move + { + Vector_3 move; + int neighbors; + FT mass; + }; - 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; + std::unordered_map m_vertex_id; + std::vector m_free_vertices; + bool m_flip_smooth_steps = false; + std::vector m_moves; + + using Incident_cells_vector = boost::container::small_vector; + std::vector m_inc_cells; using Segment_vec = std::vector; using Segment_iter = typename Segment_vec::iterator; @@ -68,71 +79,68 @@ class Tetrahedral_remeshing_smoother using AABB_segment_traits = CGAL::AABB_traits_3; using AABB_segment_tree = CGAL::AABB_tree; -private: - typedef CGAL::Tetrahedral_remeshing::internal::FMLS FMLS; - std::vector subdomain_FMLS; - std::unordered_map> subdomain_FMLS_indices; + 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; Triangle_vec m_aabb_triangles; AABB_triangle_tree m_triangles_aabb_tree; + FT m_aabb_epsilon; + Segment_vec m_aabb_segments; AABB_segment_tree m_segments_aabb_tree; - FT m_aabb_epsilon; - const SizingFunction& m_sizing; + 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; + + FT m_total_move = 0; + const CellSelector& m_cell_selector; + const SizingFunction& m_sizing; const bool m_protect_boundaries; const bool m_smooth_constrained_edges; - // the 2 following variables become useful and valid - // just before flip/smooth steps, when no vertices get inserted - // nor removed anymore - std::unordered_map m_vertex_id; - std::vector m_free_vertices{}; - bool m_flip_smooth_steps{false}; - - struct Move - { - Vector_3 move; - int neighbors; - FT mass; - }; + using FMLS = CGAL::Tetrahedral_remeshing::internal::FMLS; + std::vector subdomain_FMLS; + std::unordered_map> subdomain_FMLS_indices; -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) + 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 - createMLSSurfaces(subdomain_FMLS, - subdomain_FMLS_indices, - vertices_normals, - vertices_surface_indices, - c3t3); + createMLSSurfaces(subdomain_FMLS, subdomain_FMLS_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,14 +153,7 @@ 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); @@ -167,8 +168,12 @@ class Tetrahedral_remeshing_smoother return m_free_vertices[vid]; } - // this function can be used iff m_vertex_id - // has already been initialized +private: + bool is_selected(const Cell_handle c) const + { + return get(m_cell_selector, c); + } + void reset_free_vertices(const Tr& tr) { // when flip-smooth steps start, @@ -240,29 +245,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 +294,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 +325,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 +397,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,11 +419,11 @@ 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 (typename Vertices_normals_map::iterator vnm_it = m_vertices_normals.begin(); + vnm_it != m_vertices_normals.end(); ++vnm_it) { //value type is map - for (typename VertexNormalsMap::mapped_type::iterator it = vnm_it->second.begin(); + for (typename Vertices_normals_map::mapped_type::iterator it = vnm_it->second.begin(); it != vnm_it->second.end(); ++it) { Vector_3& n = it->second; @@ -462,44 +459,94 @@ 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); + } + } + } + + 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++; + } + } +}; - const FMLS& fmls = subdomain_FMLS[subdomain_FMLS_indices.at(si)]; - int it_nb = 0; - const int max_it_nb = 5; - const double epsilon = fmls.getPNScale() / 1000.; - const double sq_eps = CGAL::square(epsilon); +template +class VertexSmoothOperationBase +{ +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; - do - { - point = result; + 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; - fmls.fastProjectionCPU(point, result, res_normal); + 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; - 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 Context = Vertex_smoothing_context; - return Point_3(result[0], result[1], result[2]); +public: + std::shared_ptr m_context{nullptr}; + bool m_smooth_constrained_edges; + + VertexSmoothOperationBase(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; + } + + 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; } - Dihedral_angle_cosine max_cosine(const Tr& tr, - const boost::container::small_vector& cells) + 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. + Dihedral_angle_cosine max_cos_dh = cosine_of_90_degrees(); for (Cell_handle c : cells) { if(!is_selected(c)) @@ -519,12 +566,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(); const typename Tr::Geom_traits::Point_3 pv = point(backup); bool valid_orientation = false; @@ -532,40 +579,35 @@ 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 - ? max_cosine(tr, inc_cells) - : Dihedral_angle_cosine(CGAL::ZERO, 0., 1.);//Dummy unused value + 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.); bool valid_try = true; do { v->set_point(typename Tr::Point(pv + frac * move)); - valid_try = true; valid_orientation = true; angles_improved = true; for (const typename Tr::Cell_handle& ci : inc_cells) { - if (CGAL::POSITIVE != CGAL::orientation(point(ci->vertex(0)->point()), - point(ci->vertex(1)->point()), - point(ci->vertex(2)->point()), - point(ci->vertex(3)->point()))) + if (CGAL::POSITIVE != CGAL::orientation(point(ci->vertex(0)->point()), point(ci->vertex(1)->point()), + point(ci->vertex(2)->point()), point(ci->vertex(3)->point()))) { frac = 0.5 * frac; valid_try = false; 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)) { Dihedral_angle_cosine max_cos_ci = max_cos_dihedral_angle(tr, ci, false); if (curr_max_cos < max_cos_ci) { - // keep move only if new cosine is smaller than previous one - // i.e. if angle is larger frac = 0.5 * frac; valid_try = false; angles_improved = false; @@ -574,98 +616,108 @@ class Tetrahedral_remeshing_smoother } } } - } - while(!valid_try && frac > 0.1); - - // if move failed, cancel move - bool valid_move = valid_orientation && angles_improved; + } while(!valid_try && frac > 0.1); + bool valid_move = valid_orientation && angles_improved; if(!valid_move) v->set_point(backup); - #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE else total_move += CGAL::approximate_sqrt(CGAL::squared_distance(pv, point(v->point()))); #endif - return valid_move; } +}; + +template +class ComplexEdgeVertexSmoothOperation + : public VertexSmoothOperationBase, + public ElementaryOperation +{ +public: + using BaseClass = VertexSmoothOperationBase; + using Vertex_handle = typename C3t3::Triangulation::Vertex_handle; + using Base = ElementaryOperation; + using ElementType = Vertex_handle; + using ElementSource = typename Base::ElementSource; - 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); + using Surface_patch_index = typename C3t3::Surface_patch_index; - 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); - } - } - } + using BaseClass::m_context; - 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++; - } - } + 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; +public: + ComplexEdgeVertexSmoothOperation(std::shared_ptr context) + : BaseClass(context) {} - // 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 + ElementSource get_element_source(const C3t3& c3t3) const override { - 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; + perform_global_preprocessing(c3t3); + return c3t3.triangulation().finite_vertex_handles(); } - template - std::size_t smooth_edges_in_complex(C3t3& c3t3, + bool execute_operation(const ElementType& v, C3t3& c3t3) override + { + 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; + + 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; + #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(); - 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); + m_context->m_moves.assign(nbv, typename BaseClass::Context::Move{CGAL::NULL_VECTOR, 0, 0.}); - //collect neighbors for (const Edge& e : c3t3.edges_in_complex()) { const Vertex_handle vh0 = e.first->vertex(e.second); @@ -674,485 +726,372 @@ 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) { - moves[i0].move += density * Vector_3(p0, p1); - moves[i0].mass += density; - ++moves[i0].neighbors; + m_context->m_moves[i0].move += density * Vector_3(p0, p1); + m_context->m_moves[i0].mass += density; + ++m_context->m_moves[i0].neighbors; } if (vh1_moving) { - moves[i1].move += density * Vector_3(p1, p0); - moves[i1].mass += density; - ++moves[i1].neighbors; + m_context->m_moves[i1].move += density * Vector_3(p1, p0); + m_context->m_moves[i1].mass += density; + ++m_context->m_moves[i1].neighbors; } } + } +}; + +template +class SurfaceVertexSmoothOperation + : public VertexSmoothOperationBase, + public ElementaryOperation +{ +public: + using BaseClass = VertexSmoothOperationBase; + using Base = ElementaryOperation; + using ElementType = typename C3t3::Triangulation::Vertex_handle; + using ElementSource = typename Base::ElementSource; + + 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(); + const std::size_t nbv = tr.number_of_vertices(); + const typename BaseClass::Context::Move default_move{CGAL::NULL_VECTOR, 0, 0.}; + m_context->m_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)); - - sum_projections += Vector_3(tmp_pos, normal_projection); - tmp_pos = normal_projection; - } -#endif //CGAL_TET_REMESHING_EDGE_SMOOTHING_DISABLE_PROJECTION + const Vertex_handle vh0 = e.first->vertex(e.second); + const Vertex_handle vh1 = e.first->vertex(e.third); - const Point_3 new_pos = current_pos + sum_projections; + const std::size_t& i0 = m_context->vertex_id(vh0); + const std::size_t& i1 = m_context->vertex_id(vh1); -#else // AABB_tree projection + 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); - 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) + { + m_context->m_moves[i0].move += density * Vector_3(p0, p1); + m_context->m_moves[i0].mass += density; + ++m_context->m_moves[i0].neighbors; + } + if (vh1_moving) + { + m_context->m_moves[i1].move += density * Vector_3(p1, p0); + m_context->m_moves[i1].mass += density; + ++m_context->m_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: + SurfaceVertexSmoothOperation(std::shared_ptr context) + : BaseClass(context) {} + + ElementSource get_element_source(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 ElementType& v, C3t3& c3t3) override + { + auto& tr = c3t3.triangulation(); + 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 std::size_t nb_neighbors = m_context->m_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 Vector_3 move = m_context->m_moves[vid].move / m_context->m_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) { const auto op1 = get_intersection_point(proj); const auto op2 = get_intersection_point(proj_opp); - - const FT sqd1 = (op1 == std::nullopt) ? 0. - : CGAL::squared_distance(smoothed_position, op1.value()); - const FT sqd2 = (op2 == std::nullopt) ? 0. - : CGAL::squared_distance(smoothed_position, op2.value()); - - if (sqd1 != 0. && sqd1 < sqd2) + const FT sqd1 = (op1 == std::nullopt) ? 0. : CGAL::squared_distance(smoothed_position, op1.value()); + const FT sqd2 = (op2 == std::nullopt) ? 0. : CGAL::squared_distance(smoothed_position, op2.value()); + if(sqd1 != 0. && sqd1 < sqd2) new_pos = op1.value(); - else if (sqd2 != 0) + else if(sqd2 != 0) new_pos = op2.value(); else new_pos = smoothed_position; } - else //no valid projection + else + { new_pos = smoothed_position; + } } -#endif //CGAL_TET_REMESHING_SMOOTHING_WITH_MLS +#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; - - const Point_3 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 - - 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; + if(mls_proj == std::nullopt) + return false; + new_pos = *mls_proj; +#else + new_pos = m_context->m_segments_aabb_tree.closest_point(current_pos); #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 InternalVertexSmoothOperation + : public VertexSmoothOperationBase, + public ElementaryOperation { - std::size_t nb_done_3d = 0; - auto& tr = c3t3.triangulation(); +public: + using BaseClass = VertexSmoothOperationBase; + using Base = ElementaryOperation; + using ElementType = typename C3t3::Triangulation::Vertex_handle; + using ElementSource = typename Base::ElementSource; + + 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: + InternalVertexSmoothOperation(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(); + + const std::size_t nbv = tr.number_of_vertices(); + const typename BaseClass::Context::Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/}; + m_context->m_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) + { + m_context->m_moves[i0].move += density * Vector_3(p0, p1); + m_context->m_moves[i0].mass += density; + ++m_context->m_moves[i0].neighbors; + } + if (vh1_moving) + { + m_context->m_moves[i1].move += density * Vector_3(p1, p0); + m_context->m_moves[i1].mass += density; + ++m_context->m_moves[i1].neighbors; + } } } } - // iterate over vertices and move - for(Vertex_handle v : tr.finite_vertex_handles()) + ElementSource get_element_source(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 ElementType& 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); - - 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 - ); - } + auto& tr = c3t3.triangulation(); + 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 && m_context->m_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 = m_context->m_moves[vid].move / m_context->m_moves[vid].mass; + const Point_3 new_pos = point(v->point()) + move; + return BaseClass::check_inversion_and_move(v, new_pos, m_context->m_inc_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/tetrahedral_adaptive_remeshing_impl.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h index 896e3c1e982..64da5630e3f 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 @@ -87,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; @@ -95,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; @@ -117,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"); @@ -148,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"); @@ -254,7 +254,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) + { + ComplexEdgeVertexSmoothOperation op(m_smoothing_context); + ElementaryOperationExecutionSequential executor; + executor.execute(op, m_c3t3); + } + { + SurfaceVertexSmoothOperation op(m_smoothing_context); + ElementaryOperationExecutionSequential executor; + executor.execute(op, m_c3t3); + } + } + { + InternalVertexSmoothOperation op(m_smoothing_context); + ElementaryOperationExecutionSequential executor; + executor.execute(op, m_c3t3); + } #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL_assertion(tr().tds().is_valid(true)); @@ -650,7 +672,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; From 6c49bad30fc6f3e55d4ba1df83195648c9f7e7b2 Mon Sep 17 00:00:00 2001 From: IasonManolas Date: Wed, 8 Jul 2026 11:31:48 +0300 Subject: [PATCH 5/5] Tetrahedral_remeshing: refactor edge collapse onto the elementary-operation framework Replace the collapse_short_edges() free function with EdgeCollapseOperation. get_element_source() collects and length-orders (shortest first) the collapsible edges once; execute_operation() collapses one edge, skipping those invalidated by an earlier collapse via a should_skip map. Unlike the former dynamic bimap worklist, the candidate list is static: edges that become short during a pass are not re-collapsed in that pass but in the next remeshing iteration. This changes the result (it is not byte-identical), but Dihedral_Angle.Mean is unchanged within noise (+0.024% on bear). The shared collapse_edge()/collapse() helpers now mark invalidated edges through remove_from_bimap (repurposed to set the skip flag); the now-unused update_bimap is removed. --- .../internal/collapse_short_edges.h | 186 +++++++----------- .../tetrahedral_adaptive_remeshing_impl.h | 6 +- .../internal/tetrahedral_remeshing_helpers.h | 31 +-- 3 files changed, 86 insertions(+), 137 deletions(-) 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 0b5c2f89533..1235e6ce5c7 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 @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -30,6 +31,7 @@ #include #include +#include #include #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE @@ -976,7 +978,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, @@ -1083,7 +1085,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 */, @@ -1230,125 +1232,89 @@ auto can_be_collapsed(const typename C3T3::Edge& e, return Collapsible {true, boundary}; } -template -void collapse_short_edges(C3T3& c3t3, - const Sizing& sizing, - const bool protect_boundaries, - CellSelector cell_selector, - Visitor& visitor) +template +class EdgeCollapseOperation + : public ElementaryOperation> { - 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 = std::vector; + using Base = ElementaryOperation; + using ElementType = typename Base::ElementType; + using ElementSource = typename Base::ElementSource; + +private: + const SizingFunction& m_sizing; + const CellSelector& m_cell_selector; + bool m_protect_boundaries; + Visitor& m_visitor; + + // Edges invalidated by an earlier collapse in this pass. The candidate list + // is collected once, so an edge whose cells were destroyed by a preceding + // collapse is marked here (by collapse_edge, through remove_from_bimap) and + // skipped when the pass reaches it -- this replaces the former dynamic bimap + // worklist. Edges that only *become* short during the pass are not + // re-collapsed here; they are handled in the next remeshing iteration. + boost::unordered_map m_should_skip; - //collect long edges - Boost_bimap short_edges; - for (const Edge& e : tr.finite_edges()) +public: + EdgeCollapseOperation(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_element_source(const C3t3& c3t3) const override { - auto [collapsible, boundary] = can_be_collapsed(e, c3t3, protect_boundaries, cell_selector); - if (!collapsible) - continue; + std::vector> short_edges_with_length; + const Tr& tr = c3t3.triangulation(); - const auto sqlen = is_too_short(e, boundary, sizing, c3t3, cell_selector); - if(sqlen != std::nullopt) - short_edges.insert(short_edge(e, sqlen.value())); - } + for (const Edge& e : tr.finite_edges()) + { + auto [collapsible, boundary] = can_be_collapsed(e, c3t3, m_protect_boundaries, m_cell_selector); + if (!collapsible) + continue; -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - debug::dump_edges(short_edges, "short_edges.polylines.txt"); + const auto sqlen = is_too_short(e, boundary, m_sizing, c3t3, m_cell_selector); + if (sqlen != std::nullopt) + short_edges_with_length.push_back(std::make_pair(e, sqlen.value())); + } - 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 + // shortest first; stable to match the original bimap's ordering + std::stable_sort(short_edges_with_length.begin(), short_edges_with_length.end(), + [](const std::pair& a, const std::pair& b) { + return a.second < b.second; + }); + + Short_edges short_edges; + short_edges.reserve(short_edges_with_length.size()); + for (const auto& ef : short_edges_with_length) + short_edges.push_back(ef.first); + return short_edges; + } - while(!short_edges.empty()) + bool execute_operation(const ElementType& edge, C3t3& c3t3) override { - //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 - - short_edges.right.erase(eit); - - 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)); - -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - const auto p1 = e.first->vertex(e.second)->point(); - const auto p2 = e.first->vertex(e.third)->point(); -#endif - - 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; - - const auto sqlen = is_too_short(eshort, boundary, sizing, c3t3, cell_selector); - update_bimap(eshort, short_edges, sqlen); - } + if (m_should_skip.find(edge) != m_should_skip.end()) + return false; -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - ++nb_collapses; -#endif + const Vertex_handle vh = collapse_edge(edge, c3t3, m_sizing, m_protect_boundaries, + m_cell_selector, m_should_skip, m_visitor); + return (vh != Vertex_handle()); + } -#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; -#endif - } - }//end loop on short_edges -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - short_success.close(); - short_fail.close(); -#endif + std::string operation_name() const override { return "Collapse short edges"; } +}; -#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - timer.stop(); - std::cout << " done (" << nb_collapses << " collapses, in " - << timer.time() << " seconds)." << std::endl; -#endif -} } } } 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 64da5630e3f..02e1846cadd 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 @@ -198,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 EdgeCollapseOperation EdgeCollapseOp; + EdgeCollapseOp collapse_op(m_sizing, m_cell_selector, m_protect_boundaries, m_visitor); + ElementaryOperationExecutionSequential executor; + executor.execute(collapse_op, m_c3t3); #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL_assertion(tr().tds().is_valid(true)); 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 2cb6bef49ef..905cb617e04 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 @@ -1896,32 +1896,13 @@ void get_edge_info(const typename C3t3::Edge& edge, } -template -void remove_from_bimap(const typename EdgesBimap::left_map::key_type& e, - EdgesBimap& edges) +// Mark an edge invalidated by a collapse so the collapse operation skips it +// when the (statically collected) candidate list reaches it. +template +void remove_from_bimap(const EdgeType& e, + ShouldSkipContainer& should_skip) { - typename EdgesBimap::left_map::iterator eit = edges.left.find(e); - if (eit != edges.left.end()) - edges.left.erase(eit); -} - -// if e is in 'edges' -template -void -update_bimap(typename EdgesBimap::left_map::key_type& e, //Edge - EdgesBimap& edges, - const std::optional sqlen) -{ - if(sqlen == std::nullopt) - remove_from_bimap(e, edges); - else - { - typename EdgesBimap::left_map::iterator eit = edges.left.find(e); - if(eit != edges.left.end()) - edges.left.replace_data(eit, sqlen.value()); - else - edges.left.insert(typename EdgesBimap::left_map::value_type(e, sqlen.value())); - } + should_skip[e] = true; } template