From f8f213ea4e8a90fd80c5759cc6cd51a14a3d2746 Mon Sep 17 00:00:00 2001 From: davidboening Date: Wed, 21 Jun 2023 18:35:03 +0200 Subject: [PATCH] version 1.1 Added matrix caching, major speedup in cotan weights computation --- .../filter_heatgeodesic.cpp | 188 ++++++++++-------- .../filter_heatgeodesic/filter_heatgeodesic.h | 2 +- .../filter_heatgeodesic/trimesh_heatmethod.h | 139 ++++++------- 3 files changed, 181 insertions(+), 148 deletions(-) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp index 9151247a8..c780e7d07 100644 --- a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp +++ b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp @@ -173,86 +173,8 @@ std::map FilterHeatGeodesicPlugin::applyFilter(const QAct mm.updateDataMask(MeshModel::MM_VERTQUALITY); mm.updateDataMask(MeshModel::MM_VERTCOLOR); CMeshO &mesh = mm.cm; - - cb(100*0/10, "Updating Topology and Computing Face Normals..."); - vcg::tri::UpdateTopology::VertexFace(mesh); - vcg::tri::UpdateTopology::FaceFace(mesh); - vcg::tri::UpdateNormal::PerFaceNormalized(mesh); - - cb(100*1/10, "Computing Boundary Conditions..."); - Eigen::VectorXd initialConditions(mesh.VN()); - int selection_count = 0; - for(int i=0; i < mesh.vert.size(); i++){ - if (mesh.vert[i].IsS()){ - initialConditions[i] = 1; - ++selection_count; - } - else { - initialConditions[i] = 0; - } - } - if (selection_count < 1){ - log("Warning: no vertices are selected! aborting computation."); - break; - } - - - cb(100*2/10, "Building Linear System (Heatflow)..."); - Eigen::SparseMatrix mass(mesh.VN(), mesh.VN()); - buildMassMatrix(mesh, mass); - - Eigen::SparseMatrix cotanOperator(mesh.VN(), mesh.VN()); - buildCotanMatrix(mesh, cotanOperator); - - double avg_edge_len = computeAverageEdgeLength(mesh); - double timestep = parameters.getFloat("m") * avg_edge_len * avg_edge_len; - Eigen::SparseMatrix system1(mesh.VN(), mesh.VN()); - system1 = mass - timestep * cotanOperator; - - Eigen::SimplicialLDLT> solver; - - cb(100*3/10, "Computing Matrix Factorization (Heatflow)..."); - solver.compute(system1); - if(solver.info() != Eigen::Success) { - log("Error: Factorization Failed (Heatflow)."); - } - cb(100*4/10, "Solving Linear System (Heatflow)..."); - Eigen::VectorXd heatflow = solver.solve(initialConditions); // (VN) - if(solver.info() != Eigen::Success) { - log("Error: Solver Failed (Heatflow)."); - } - cb(100*5/10, "Computing Heat Gradient..."); - Eigen::MatrixX3d heatGradient = computeVertexGradient(mesh, heatflow); // (FN, 3) - - Eigen::MatrixX3d normalizedVectorField = normalizeVectorField(-heatGradient); // (FN, 3) - cb(100*6/10, "Computing Divergence..."); - Eigen::VectorXd divergence = computeVertexDivergence(mesh, normalizedVectorField); // (VN) - - cb(100*7/10, "Building Linear System (Geodesic)..."); - Eigen::SparseMatrix system2(mesh.VN(), mesh.VN()); - system2 = cotanOperator; //+ 1e-6 * Eigen::Matrix::Identity(mesh.VN(), mesh.VN()); - - cb(100*8/10, "Computing Matrix Factorization (Geodesic)..."); - solver.compute(system2); - if(solver.info() != Eigen::Success) { - log("Error: Factorization Failed (Geodesic)."); - } - cb(100*9/10, "Solving Linear System (Geodesic)..."); - Eigen::VectorXd geodesicDistance = solver.solve(divergence); - if(solver.info() != Eigen::Success) { - log("Error: Solver Failed (Geodesic)."); - } - cb(100*10/10, "Saving Geodesic Distance..."); - // invert and shift - geodesicDistance.array() *= -1; // no clue as to why this needs to be here - geodesicDistance.array() -= geodesicDistance.minCoeff(); - - // set geodesic distance as quality - for(int i=0; i < mesh.vert.size(); i++){ - mesh.vert[i].Q() = geodesicDistance(i); - } - - vcg::tri::UpdateColor::PerVertexQualityRamp(mesh); + // TODO: compact all vectors + computeHeatGeodesicFromSelection(mesh, cb, parameters.getFloat("m")); } break; default : @@ -261,4 +183,110 @@ std::map FilterHeatGeodesicPlugin::applyFilter(const QAct return std::map(); } +inline void FilterHeatGeodesicPlugin::computeHeatGeodesicFromSelection(CMeshO& mesh, vcg::CallBackPos* cb, float m){ + // build boundary conditions + Eigen::VectorXd boundaryConditions(mesh.VN()); + int selection_count = 0; + for(int i=0; i < mesh.VN(); i++){ + boundaryConditions(i) = mesh.vert[i].IsS() ? (++selection_count, 1) : 0; + } + if (selection_count < 1){ + log("Warning: no vertices are selected! aborting computation."); + return; + } + + // 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 avg_edge_len; + + typedef std::tuple< + Eigen::SparseMatrix, // system1 fact + Eigen::SparseMatrix, // system2 fact + double // avg edge len + > HeatMethodData; + + // recover state if it exists + if (!vcg::tri::HasPerMeshAttribute(mesh, "HeatMethodData")){ + // current solution saves matrices not factorizations + // if we save factorization it is best to at least also save cotanMatrix + // to avoid an expensive reconstruction when parameter("m") changes + buildMassMatrix(mesh, massMatrix); + // this step is very slow (95% of total time) hence we pass the callback + // to update progress and avoid freezes + buildCotanLowerTriMatrix(mesh, cotanMatrix, cb); + avg_edge_len = computeAverageEdgeLength(mesh); + + HeatMethodData &HMData = vcg::tri::Allocator::GetPerMeshAttribute(mesh, std::string("HeatMethodData"))(); + HMData = HeatMethodData(massMatrix, cotanMatrix, avg_edge_len); + } + else { + HeatMethodData &HMData = vcg::tri::Allocator::GetPerMeshAttribute(mesh, std::string("HeatMethodData"))(); + std::tie (massMatrix, cotanMatrix, avg_edge_len) = HMData; + } + + // cholesky type solver + Eigen::SimplicialLDLT> solver; + + // build system 1 + double timestep = m * avg_edge_len * avg_edge_len; + Eigen::SparseMatrix system1(mesh.VN(), mesh.VN()); + system1 = massMatrix - timestep * cotanMatrix; + + cb(91, "Computing Factorization 1..."); + + // factorize system 1 + solver.compute(system1); + if(solver.info() != Eigen::Success) { + log("Error: Factorization Failed (Heatflow)."); + } + + cb(93, "Solving System 1..."); + // solve system 1 + Eigen::VectorXd heatflow = solver.solve(boundaryConditions); // (VN) + if(solver.info() != Eigen::Success) { + log("Error: Solver Failed (Heatflow)."); + } + + cb(95, "Computing Gradient, VectorField, Divergence..."); + // intermediate steps + Eigen::MatrixX3d heatGradient = computeVertexGradient(mesh, heatflow); // (FN, 3) + Eigen::MatrixX3d normalizedVectorField = normalizeVectorField(-heatGradient); // (FN, 3) + Eigen::VectorXd divergence = computeVertexDivergence(mesh, normalizedVectorField); // (VN) + + // build system 2 + Eigen::SparseMatrix system2(mesh.VN(), mesh.VN()); + system2 = cotanMatrix; //+ 1e-6 * Eigen::Matrix::Identity(mesh.VN(), mesh.VN()); + + cb(96, "Computing Factorization 2..."); + // factorize system 2 + solver.compute(system2); + if(solver.info() != Eigen::Success) { + log("Error: Factorization Failed (Geodesic)."); + } + + cb(97, "Solving System 2..."); + // solve system 2 + Eigen::VectorXd geodesicDistance = solver.solve(divergence); + + if(solver.info() != Eigen::Success) { + log("Error: Solver Failed (Geodesic)."); + } + + // shift to impose boundary conditions (dist(d) = 0 \forall d \in init_cond) + geodesicDistance.array() -= geodesicDistance.minCoeff(); + + cb(99, "Updating Quality and Color..."); + // set geodesic distance as quality and color + for(int i=0; i < mesh.vert.size(); i++){ + mesh.vert[i].Q() = geodesicDistance(i); + } + vcg::tri::UpdateColor::PerVertexQualityRamp(mesh); +} + MESHLAB_PLUGIN_NAME_EXPORTER(FilterSamplePlugin) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h index 2fb67ad6d..b421550ab 100644 --- a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h +++ b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h @@ -33,7 +33,7 @@ public: vcg::CallBackPos * cb); private: - void computeHeatGeodesicFromSelection(MeshDocument& mesh, vcg::CallBackPos* cb, float m); + void computeHeatGeodesicFromSelection(CMeshO& mesh, vcg::CallBackPos* cb, float m); }; #endif diff --git a/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h b/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h index 3f3eab0e0..37ebe5d2e 100644 --- a/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h +++ b/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h @@ -9,30 +9,26 @@ #include -// foward declarations -inline Eigen::Vector3d toEigen(const vcg::Point3f& p); -inline double cotan(const Eigen::Vector3d& v0, const Eigen::Vector3d& v1); -inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix &mass); -inline void buildCotanMatrix(CMeshO &mesh, Eigen::SparseMatrix &cotanOperator); -inline double computeAverageEdgeLength(CMeshO &mesh); -inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorXd &heat); -inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::MatrixX3d &field); -inline Eigen::MatrixX3d normalizeVectorField(const Eigen::MatrixX3d &field); +// LOCAL FUNCTIONS -inline Eigen::Vector3d toEigen(const vcg::Point3f& p) -{ - return Eigen::Vector3d(p.X(), p.Y(), p.Z()); -}; +namespace { + inline Eigen::Vector3d toEigen(const vcg::Point3f& p) + { + return Eigen::Vector3d(p.X(), p.Y(), p.Z()); + }; -inline double cotan(const Eigen::Vector3d& v0, const Eigen::Vector3d& v1) -{ - // cos(theta) / sin(theta) - return v0.dot(v1) / v0.cross(v1).norm(); -}; + inline double cotan(const Eigen::Vector3d& v0, const Eigen::Vector3d& v1) + { + // cos(theta) / sin(theta) + return v0.dot(v1) / v0.cross(v1).norm(); + }; +} +// GLOBAL FUNCTIONS inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix &mass){ + mass.resize(mesh.VN(), mesh.VN()); // compute area of all faces for (CMeshO::FaceIterator fi = mesh.face.begin(); fi != mesh.face.end(); ++fi) { @@ -55,7 +51,7 @@ inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix &mass){ std::vector faces; std::vector indices; vcg::face::VFStarVF(vp, faces, indices); - + double area = 0; for (int j = 0; j < faces.size(); ++j) { @@ -67,55 +63,65 @@ inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix &mass){ } -inline void buildCotanMatrix(CMeshO &mesh, Eigen::SparseMatrix &cotanOperator){ +inline void buildCotanLowerTriMatrix(CMeshO &mesh, Eigen::SparseMatrix &cotanOperator, vcg::CallBackPos* cb = NULL){ + cotanOperator.resize(mesh.VN(), mesh.VN()); // initialize a hashtable from vertex pointers to ids std::unordered_map vertex_ids; for (int i = 0; i < mesh.VN(); ++i){ vertex_ids[&mesh.vert[i]] = i; } - // iterate over all vertices to fill cotan matrix - for (int i = 0; i < mesh.VN(); ++i){ - CMeshO::VertexType *vp = &mesh.vert[i]; - CMeshO::FaceType *fp = vp->VFp(); - vcg::face::Pos pos(fp, vp); - vcg::face::Pos start(fp, vp); - // iterate over all incident edges of vp - do { - // get the vertex opposite to vp - pos.FlipV(); - CMeshO::VertexType *vo = pos.V(); - // move to left vertex - pos.FlipE();pos.FlipV(); - CMeshO::VertexType *vl = pos.V(); - // move back then to right vertex - pos.FlipV();pos.FlipE(); // back to vo - pos.FlipF();pos.FlipE();pos.FlipV(); - CMeshO::VertexType *vr = pos.V(); - pos.FlipV();pos.FlipE();pos.FlipF();pos.FlipV(); // back to vp + // progress bar variables + int update_size = mesh.FN() / 90; // 90 updates (90% weight of the progress bar) + int progress = 0; - // compute cotan of left edges and right edges - Eigen::Vector3d elf = toEigen(vo->P() - vl->P()); // far left edge - Eigen::Vector3d eln = toEigen(vp->P() - vl->P()); // near left edge - Eigen::Vector3d erf = toEigen(vp->P() - vr->P()); // far right edge - Eigen::Vector3d ern = toEigen(vo->P() - vr->P()); // near right edge + // compute cotan weights + for (int i = 0; i < mesh.FN(); ++i){ + // progress bar update + if (cb != NULL && i % update_size == 89){ + cb(++progress, "Computing Cotan Weights..."); + } + CMeshO::FaceType* fi = &mesh.face[i]; - double cotan_l = cotan(elf, eln); - double cotan_r = cotan(ern, erf); + CMeshO::VertexType* v0 = fi->V(0); + CMeshO::VertexType* v1 = fi->V(1); + CMeshO::VertexType* v2 = fi->V(2); - // add to the matrix - cotanOperator.coeffRef(vertex_ids[vp], vertex_ids[vo]) = (cotan_l + cotan_r)/2; + vcg::Point3f p0 = v0->P(); + vcg::Point3f p1 = v1->P(); + vcg::Point3f p2 = v2->P(); - // move to the next edge - pos.FlipF();pos.FlipE(); - } while (pos != start); + Eigen::Vector3d e0 = toEigen(p2 - p1); + Eigen::Vector3d e1 = toEigen(p0 - p2); + Eigen::Vector3d e2 = toEigen(p1 - p0); + + // first edge is inverted to get correct orientation + double alpha0 = cotan(-e1, e2) / 2; + double alpha1 = cotan(-e2, e0) / 2; + double alpha2 = cotan(-e0, e1) / 2; + + int i0 = vertex_ids[v0]; + int i1 = vertex_ids[v1]; + int i2 = vertex_ids[v2]; + + // save only lower triangular part + if (i0 > i1) + cotanOperator.coeffRef(i0, i1) += alpha2; + else + cotanOperator.coeffRef(i1, i0) += alpha2; + if (i0 > i2) + cotanOperator.coeffRef(i0, i2) += alpha1; + else + cotanOperator.coeffRef(i2, i0) += alpha1; + if (i1 > i2) + cotanOperator.coeffRef(i1, i2) += alpha0; + else + cotanOperator.coeffRef(i2, i1) += alpha0; + + cotanOperator.coeffRef(i0, i0) -= (alpha1 + alpha2); + cotanOperator.coeffRef(i1, i1) -= (alpha0 + alpha2); + cotanOperator.coeffRef(i2, i2) -= (alpha0 + alpha1); } - - // compute diagonal entries - for (int i = 0; i < mesh.VN(); ++i){ - cotanOperator.coeffRef(i, i) = -cotanOperator.row(i).sum(); - } - } @@ -156,8 +162,7 @@ inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorX n /= n.norm(); // face area double faceArea = fp->Q(); - // (ORDERING): edge unit vectors (assuming counter-clockwise ordering) - // note if the ordering is clockwise, the gradient will point in the opposite direction + // edge unit vectors (counter-clockwise) Eigen::Vector3d e0 = toEigen(p2 - p1); e0 /= e0.norm(); Eigen::Vector3d e1 = toEigen(p0 - p2); @@ -170,13 +175,13 @@ inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorX Eigen::Vector3d g2 = n.cross(e2); //v2 grad // add vertex gradient contributions - Eigen::Vector3d total_grad = ( - g0 * heat(vertex_ids[fp->V(0)]) + - g1 * heat(vertex_ids[fp->V(1)]) + - g2 * heat(vertex_ids[fp->V(2)]) + Eigen::Vector3d tri_grad = ( + g0 * heat(vertex_ids[fp->V(0)]) + + g1 * heat(vertex_ids[fp->V(1)]) + + g2 * heat(vertex_ids[fp->V(2)]) ) / (2 * faceArea); - heatGradientField.row(i) = total_grad; + heatGradientField.row(i) = tri_grad; } return heatGradientField; } @@ -227,14 +232,14 @@ inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::Matrix el = toEigen(p0 - p1); //-e2 er = toEigen(p2 - p1); //e0 eo = toEigen(p0 - p2); //e1 - } else if (index == 2){ + } else if (index == 2){ el = toEigen(p1 - p2); //-e0 er = toEigen(p0 - p2); //e1 eo = toEigen(p1 - p0); //e2 } // compute left and right cotangents - double cotl = cotan(-el, eo); // -el -> angle between el and eo - double cotr = cotan(er, eo); + double cotl = cotan(-el, -eo); + double cotr = cotan(-er, eo); // normalize edge vectors after cotangent computation el /= el.norm(); er /= er.norm();