From f14dca369360d0ade59eeb2bda42887876827811 Mon Sep 17 00:00:00 2001 From: davidboening Date: Wed, 12 Jul 2023 09:39:22 +0200 Subject: [PATCH] version 1.4.1 - removed matrix caching - updated description - minor refactoring --- .../filter_heatgeodesic.cpp | 62 ++++++++----------- .../filter_heatgeodesic/trimesh_heatmethod.h | 46 +++++++++----- 2 files changed, 54 insertions(+), 54 deletions(-) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp index a2da10cf0..ec435b462 100644 --- a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp +++ b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp @@ -72,7 +72,12 @@ QString FilterHeatGeodesicPlugin::pythonFilterName(ActionIDType f) const { switch(filterId) { case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - return QString("Computes an approximated geodesic distance from the selected vertices to all others. This algorithm implements the heat method."); + return QString( + "Computes an approximated geodesic distance from the selected vertices to all others. " + "This algorithm implements the heat method. " + "If the mesh has changed, update m to a different value to reflect these changes as" + " this rebuilds the existing cache." + ); default : assert(0); return QString("Unknown Filter"); @@ -198,51 +203,34 @@ inline void FilterHeatGeodesicPlugin::computeHeatGeodesicFromSelection(CMeshO& m return; } + // TODO detect if the mesh has changed and if so delete state + // NOW changing the value of m rebuilds the cache + // + // if (meshHasChanged) { + // auto handle = vcg::tri::Allocator::GetPerMeshAttribute>(mesh, std::string("HeatMethodSolver")); + // vcg::tri::Allocator::DeletePerMeshAttribute>(mesh, handle); + // } + std::shared_ptr HMSolver; - - // delete state if the mesh has changed - bool meshHasChanged = false; // Update with appropriate function - if (meshHasChanged) { - auto handle = vcg::tri::Allocator::GetPerMeshAttribute>(mesh, std::string("HeatMethodSolver")); - vcg::tri::Allocator::DeletePerMeshAttribute>(mesh, handle); - } - - // recover state if it exists + // if no state can be recovered, initialize it if (!vcg::tri::HasPerMeshAttribute(mesh, "HeatMethodSolver")){ - cb(10, "Updating Topology..."); - // update topology and face normals - vcg::tri::UpdateTopology::VertexFace(mesh); - vcg::tri::UpdateTopology::FaceFace(mesh); - vcg::tri::UpdateNormal::PerFaceNormalized(mesh); - - // state variables - Eigen::SparseMatrix massMatrix; - Eigen::SparseMatrix cotanMatrix; - double averageEdgeLength; - - buildMassMatrix(mesh, massMatrix); - // this step is very slow (95% of total time) hence we pass the callback - // to update progress and avoid freezes - cb(20, "Building Cotan Matrix..."); - buildCotanLowerTriMatrix(mesh, cotanMatrix); - averageEdgeLength = computeAverageEdgeLength(mesh); - - cb(70, "Matrix Factorization..."); - HMSolver = std::make_shared(HeatMethodSolver(std::move(massMatrix), std::move(cotanMatrix), averageEdgeLength, m)); + HMSolver = std::make_shared(HeatMethodSolver(mesh, m, cb)); if (HMSolver->factorization_failed()) - log("Warning: factorization has failed. The mesh is non-manifold or badly conditioned (e.g. angles ~ 0deg) or has disconnected components"); + log("Warning: factorization has failed. The mesh is non-manifold or badly conditioned (e.g. angles ~ 0deg) or disconnected components"); auto HMSolverHandle = vcg::tri::Allocator::GetPerMeshAttribute>(mesh, std::string("HeatMethodSolver")); HMSolverHandle() = HMSolver; } else { - cb(10, "Recovering State..."); + cb(5, "Attempting to Recover State..."); // recover solver - HMSolver = vcg::tri::Allocator::GetPerMeshAttribute>(mesh, std::string("HeatMethodSolver"))(); - // if m has changed rebuild factorizations from existing matrices - if (HMSolver->get_m() != m){ - cb(70, "Rebuilding Factorization..."); - HMSolver->rebuildFactorization(m); + auto HMSolverHandle = vcg::tri::Allocator::GetPerMeshAttribute>(mesh, std::string("HeatMethodSolver")); + // if m has changed rebuild everything + if (HMSolverHandle()->get_m() != m){ + HMSolverHandle() = std::make_shared(HeatMethodSolver(mesh, m, cb)); + if (HMSolverHandle()->factorization_failed()) + log("Warning: factorization has failed. The changed mesh is non-manifold or badly conditioned (e.g. angles ~ 0deg) or disconnected components"); } + HMSolver = HMSolverHandle(); } cb(85, "Computing Geodesic Distance..."); diff --git a/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h b/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h index 0887341a4..348c628ca 100644 --- a/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h +++ b/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h @@ -256,34 +256,47 @@ inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::Matrix class HeatMethodSolver { public: - HeatMethodSolver(Eigen::SparseMatrix &&massMatrix, Eigen::SparseMatrix &&cotanMatrix, double averageEdgeLength, double m) - : massMatrix(std::move(massMatrix)), cotanMatrix(std::move(cotanMatrix)), averageEdgeLength(averageEdgeLength), m(m) + HeatMethodSolver(CMeshO& mesh, float m, vcg::CallBackPos* cb = NULL) + : m(m) { + if(cb != NULL) cb(10, "Updating Topology..."); + // update topology and face normals + vcg::tri::UpdateTopology::VertexFace(mesh); + vcg::tri::UpdateTopology::FaceFace(mesh); + vcg::tri::UpdateNormal::PerFaceNormalized(mesh); + + // state variables + Eigen::SparseMatrix massMatrix; + Eigen::SparseMatrix cotanMatrix; + double averageEdgeLength; + + buildMassMatrix(mesh, massMatrix); + // this step is very slow (95% of total time) hence we pass the callback + // to update progress and avoid freezes + if(cb != NULL) cb(20, "Building Cotan Matrix..."); + buildCotanLowerTriMatrix(mesh, cotanMatrix); + averageEdgeLength = computeAverageEdgeLength(mesh); + // initialize pointers solver1 = std::shared_ptr>>(new Eigen::SimplicialLDLT>); solver2 = std::shared_ptr>>(new Eigen::SimplicialLDLT>); // build factorization double timestep = m * averageEdgeLength * averageEdgeLength; + if(cb != NULL) cb(50, "Matrix Factorization..."); solver1->compute(massMatrix - timestep * cotanMatrix); + if(cb != NULL) cb(75, "Matrix Factorization..."); solver2->compute(cotanMatrix); if((solver1->info() != Eigen::Success) || (solver2->info() != Eigen::Success)) factorizationFailed = true; else factorizationFailed = false; } + HeatMethodSolver() = delete; + HeatMethodSolver& operator=(const HeatMethodSolver&) = delete; - void rebuildFactorization(double new_m){ - m = new_m; - double timestep = m * averageEdgeLength * averageEdgeLength; - solver1->compute(massMatrix - timestep * cotanMatrix); - solver2->compute(cotanMatrix); - if((solver1->info() != Eigen::Success) || (solver2->info() != Eigen::Success)) - factorizationFailed = true; - else - factorizationFailed = false; - } + Eigen::VectorXd solve(CMeshO &mesh, Eigen::VectorXd &sourcePoints){ Eigen::VectorXd heatflow = solver1->solve(sourcePoints); // (VN) Eigen::MatrixX3d heatGradient = computeVertexGradient(mesh, heatflow); // (FN, 3) @@ -294,22 +307,21 @@ public: solverFailed = true; else solverFailed = false; - // shift to impose dist(d) = 0 \forall d \in sourcePoints + // shift to impose dist(source) = 0 geodesicDistance.array() -= geodesicDistance.minCoeff(); return geodesicDistance; } + double get_m(){return m;} + bool factorization_failed(){return factorizationFailed;} + bool solver_failed(){return solverFailed;} private: // Eigen::SimplicialLDLT has no copy/move constructor std::shared_ptr>> solver1; std::shared_ptr>> solver2; - Eigen::SparseMatrix massMatrix; - Eigen::SparseMatrix cotanMatrix; - double averageEdgeLength; double m; - bool factorizationFailed; bool solverFailed; };