version 1.4.1

- removed matrix caching

- updated description

- minor refactoring
This commit is contained in:
davidboening 2023-07-12 09:39:22 +02:00
parent 4e4e950cd7
commit f14dca3693
2 changed files with 54 additions and 54 deletions

View File

@ -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<CMeshO>::GetPerMeshAttribute<std::shared_ptr<HeatMethodSolver>>(mesh, std::string("HeatMethodSolver"));
// vcg::tri::Allocator<CMeshO>::DeletePerMeshAttribute<std::shared_ptr<HeatMethodSolver>>(mesh, handle);
// }
std::shared_ptr<HeatMethodSolver> HMSolver;
// delete state if the mesh has changed
bool meshHasChanged = false; // Update with appropriate function
if (meshHasChanged) {
auto handle = vcg::tri::Allocator<CMeshO>::GetPerMeshAttribute<std::shared_ptr<HeatMethodSolver>>(mesh, std::string("HeatMethodSolver"));
vcg::tri::Allocator<CMeshO>::DeletePerMeshAttribute<std::shared_ptr<HeatMethodSolver>>(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<CMeshO>::VertexFace(mesh);
vcg::tri::UpdateTopology<CMeshO>::FaceFace(mesh);
vcg::tri::UpdateNormal<CMeshO>::PerFaceNormalized(mesh);
// state variables
Eigen::SparseMatrix<double> massMatrix;
Eigen::SparseMatrix<double> 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>(HeatMethodSolver(std::move(massMatrix), std::move(cotanMatrix), averageEdgeLength, m));
HMSolver = std::make_shared<HeatMethodSolver>(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<CMeshO>::GetPerMeshAttribute<std::shared_ptr<HeatMethodSolver>>(mesh, std::string("HeatMethodSolver"));
HMSolverHandle() = HMSolver;
}
else {
cb(10, "Recovering State...");
cb(5, "Attempting to Recover State...");
// recover solver
HMSolver = vcg::tri::Allocator<CMeshO>::GetPerMeshAttribute<std::shared_ptr<HeatMethodSolver>>(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<CMeshO>::GetPerMeshAttribute<std::shared_ptr<HeatMethodSolver>>(mesh, std::string("HeatMethodSolver"));
// if m has changed rebuild everything
if (HMSolverHandle()->get_m() != m){
HMSolverHandle() = std::make_shared<HeatMethodSolver>(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...");

View File

@ -256,34 +256,47 @@ inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::Matrix
class HeatMethodSolver {
public:
HeatMethodSolver(Eigen::SparseMatrix<double> &&massMatrix, Eigen::SparseMatrix<double> &&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<CMeshO>::VertexFace(mesh);
vcg::tri::UpdateTopology<CMeshO>::FaceFace(mesh);
vcg::tri::UpdateNormal<CMeshO>::PerFaceNormalized(mesh);
// state variables
Eigen::SparseMatrix<double> massMatrix;
Eigen::SparseMatrix<double> 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<Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>(new Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>);
solver2 = std::shared_ptr<Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>(new Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>);
// 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<Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>> solver1;
std::shared_ptr<Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>> solver2;
Eigen::SparseMatrix<double> massMatrix;
Eigen::SparseMatrix<double> cotanMatrix;
double averageEdgeLength;
double m;
bool factorizationFailed;
bool solverFailed;
};