From 8ef1e76548e7e5cb7b869239b4e5e7f85ab4ded0 Mon Sep 17 00:00:00 2001 From: davidboening Date: Sat, 17 Jun 2023 22:27:09 +0200 Subject: [PATCH 01/12] added heat geodesic (not tested) --- .../filter_heatgeodesic/CMakeLists.txt | 5 + .../filter_heatgeodesic.cpp | 247 +++++++++++++++++ .../filter_heatgeodesic/filter_heatgeodesic.h | 39 +++ .../filter_heatgeodesic/trimesh_heatmethod.h | 250 ++++++++++++++++++ 4 files changed, 541 insertions(+) create mode 100644 src/meshlabplugins/filter_heatgeodesic/CMakeLists.txt create mode 100644 src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp create mode 100644 src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h create mode 100644 src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h diff --git a/src/meshlabplugins/filter_heatgeodesic/CMakeLists.txt b/src/meshlabplugins/filter_heatgeodesic/CMakeLists.txt new file mode 100644 index 000000000..732bcaacc --- /dev/null +++ b/src/meshlabplugins/filter_heatgeodesic/CMakeLists.txt @@ -0,0 +1,5 @@ +set(SOURCES filter_heatgeodesic.cpp) + +set(HEADERS filter_heatgeodesic.h trimesh_heatmethod.h) + +add_meshlab_plugin(filter_heatgeodesic ${SOURCES} ${HEADERS}) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp new file mode 100644 index 000000000..ececbd4b9 --- /dev/null +++ b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp @@ -0,0 +1,247 @@ +#include "trimesh_heatmethod.h" + +#include "filter_heatgeodesic.h" + +#include + +/** + * @brief Constructor usually performs only two simple tasks of filling the two lists + * - typeList: with all the possible id of the filtering actions + * - actionList with the corresponding actions. If you want to add icons to + * your filtering actions you can do here by construction the QActions accordingly + */ +FilterHeatGeodesicPlugin::FilterHeatGeodesicPlugin() +{ + typeList = {FP_COMPUTE_HEATGEODESIC_FROM_SELECTION}; + + for(ActionIDType tt : types()) + actionList.push_back(new QAction(filterName(tt), this)); +} + +FilterHeatGeodesicPlugin::~FilterHeatGeodesicPlugin() +{ +} + +QString FilterHeatGeodesicPlugin::pluginName() const +{ + return "FilterHeatGeodesic"; +} + +/** + * @brief ST() must return the very short string describing each filtering action + * (this string is used also to define the menu entry) + * @param filterId: the id of the filter + * @return the name of the filter + */ +QString FilterHeatGeodesicPlugin::filterName(ActionIDType filterId) const +{ + switch(filterId) { + case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : + return "Compute Geodesic From Selection (Heat)"; + default : + assert(0); + return QString(); + } +} + +/** + * @brief FilterHeatGeodesicPlugin::pythonFilterName if you want that your filter should have a different + * name on pymeshlab, use this function to return its python name. + * @param f + * @return + */ +QString FilterHeatGeodesicPlugin::pythonFilterName(ActionIDType f) const +{ + switch(f) { + case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : + return "compute_heat_geodesic_from_selection"; + default : + assert(0); + return QString(); + } +} + + +/** + * @brief // Info() must return the longer string describing each filtering action + * (this string is used in the About plugin dialog) + * @param filterId: the id of the filter + * @return an info string of the filter + */ + QString FilterHeatGeodesicPlugin::filterInfo(ActionIDType filterId) const +{ + switch(filterId) { + case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : + return "Computes an approximated geodesic distance from the selected vertices to all others. This algorithm implements the heat method."; + default : + assert(0); + return "Unknown Filter"; + } +} + + /** + * @brief The FilterClass describes in which generic class of filters it fits. + * This choice affect the submenu in which each filter will be placed + * More than a single class can be chosen. + * @param a: the action of the filter + * @return the class od the filter + */ +FilterHeatGeodesicPlugin::FilterClass FilterHeatGeodesicPlugin::getClass(const QAction *a) const +{ + switch(ID(a)) { + case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : + return FilterPlugin::Measure; + default : + assert(0); + return FilterPlugin::Generic; + } +} + +/** + * @brief FilterSamplePlugin::filterArity + * @return + */ +FilterPlugin::FilterArity FilterHeatGeodesicPlugin::filterArity(const QAction*) const +{ + return SINGLE_MESH; +} + +/** + * @brief FilterSamplePlugin::getPreConditions + * @return + */ +int FilterHeatGeodesicPlugin::getPreConditions(const QAction*) const +{ + // NOTE: note sure about these + return MeshModel::MM_VERTCOORD | MeshModel::MM_VERTQUALITY | + MeshModel::MM_FACEFACETOPO | MeshModel::MM_VERTFACETOPO | + MeshModel::MM_FACEQUALITY | MeshModel::MM_FACENORMAL | + MeshModel::MM_VERTFLAGSELECT; +} + +/** + * @brief FilterSamplePlugin::postCondition + * @return + */ +int FilterHeatGeodesicPlugin::postCondition(const QAction*) const +{ + // NOTE: note sure about these + return MeshModel::MM_VERTQUALITY | MeshModel::MM_FACEQUALITY; +} + +/** + * @brief This function define the needed parameters for each filter. Return true if the filter has some parameters + * it is called every time, so you can set the default value of parameters according to the mesh + * For each parameter you need to define, + * - the name of the parameter, + * - the default value + * - the string shown in the dialog + * - a possibly long string describing the meaning of that parameter (shown as a popup help in the dialog) + * @param action + * @param m + * @param parlst + */ +RichParameterList FilterHeatGeodesicPlugin::initParameterList(const QAction *action,const MeshModel &m) +{ + RichParameterList parlst; + switch(ID(action)) { + case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : + parlst.addParam(RichFloat("m", 1, "m", "Multiplier used in backward Euler timestep.")); + break; + default : + assert(0); + } + return parlst; +} + +/** + * @brief The Real Core Function doing the actual mesh processing. + * @param action + * @param md: an object containing all the meshes and rasters of MeshLab + * @param par: the set of parameters of each filter + * @param cb: callback object to tell MeshLab the percentage of execution of the filter + * @return true if the filter has been applied correctly, false otherwise + */ +std::map FilterHeatGeodesicPlugin::applyFilter(const QAction * action, const RichParameterList & parameters, MeshDocument &md, unsigned int& /*postConditionMask*/, vcg::CallBackPos *cb) +{ + switch(ID(action)) { + case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : + { + CMeshO &mesh = md.mm()->cm; + cb(100*0/10, "Computing Boundary Conditions..."); + Eigen::VectorXd initialConditions(mesh.VN()); + for(int i=0; i < mesh.vert.size(); i++){ + initialConditions[i] = mesh.vert[i].IsS() ? 1 : 0; + } + cb(100*1/10, "Updating Topology and Computing Face Normals..."); + vcg::tri::UpdateTopology::VertexFace(mesh); + vcg::tri::UpdateTopology::FaceFace(mesh); + vcg::tri::UpdateNormal::PerFaceNormalized(mesh); + + 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(0); + } + } + break; + default : + wrongActionCalled(action); + } + return std::map(); +} + +void computeHeatGeodesicFromSelection(MeshDocument &md, vcg::CallBackPos *cb, float m){ + +} + +MESHLAB_PLUGIN_NAME_EXPORTER(FilterSamplePlugin) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h new file mode 100644 index 000000000..2fb67ad6d --- /dev/null +++ b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h @@ -0,0 +1,39 @@ +#ifndef FILTERSAMPLE_PLUGIN_H +#define FILTERSAMPLE_PLUGIN_H + +#include + +class FilterHeatGeodesicPlugin : public QObject, public FilterPlugin +{ + Q_OBJECT + MESHLAB_PLUGIN_IID_EXPORTER(FILTER_PLUGIN_IID) + Q_INTERFACES(FilterPlugin) + +public: + enum { FP_COMPUTE_HEATGEODESIC_FROM_SELECTION } ; + + FilterHeatGeodesicPlugin(); + virtual ~FilterHeatGeodesicPlugin(); + + QString pluginName() const; + + QString filterName(ActionIDType filter) const; + QString pythonFilterName(ActionIDType f) const; + QString filterInfo(ActionIDType filter) const; + FilterClass getClass(const QAction* a) const; + FilterArity filterArity(const QAction*) const; + int getPreConditions(const QAction *) const; + int postCondition(const QAction* ) const; + RichParameterList initParameterList(const QAction*, const MeshModel &/*m*/); + std::map applyFilter( + const QAction* action, + const RichParameterList & parameters, + MeshDocument &md, + unsigned int& postConditionMask, + vcg::CallBackPos * cb); + +private: + void computeHeatGeodesicFromSelection(MeshDocument& 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 new file mode 100644 index 000000000..3f3eab0e0 --- /dev/null +++ b/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h @@ -0,0 +1,250 @@ +#ifndef TRIMESH_HEAT_METHOD +#define TRIMESH_HEAT_METHOD + +#include +#include + +#include + +#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); + +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 void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix &mass){ + // compute area of all faces + for (CMeshO::FaceIterator fi = mesh.face.begin(); fi != mesh.face.end(); ++fi) + { + vcg::Point3f p0 = fi->V(0)->P(); + vcg::Point3f p1 = fi->V(1)->P(); + vcg::Point3f p2 = fi->V(2)->P(); + double e0 = toEigen(p1 - p0).norm(); + double e1 = toEigen(p2 - p0).norm(); + double e2 = toEigen(p2 - p1).norm(); + double s = (e0 + e1 + e2) / 2; + double area = std::sqrt(s * (s - e0) * (s - e1) * (s - e2)); + // we store face areas in quality field to avoid a hash table + // this will also be useful for the gradient computation + fi->Q() = area; + } + // compute area of the dual cell for each vertex + for (int i = 0; i < mesh.VN(); ++i){ + CMeshO::VertexType *vp = &mesh.vert[i]; + + std::vector faces; + std::vector indices; + vcg::face::VFStarVF(vp, faces, indices); + + double area = 0; + for (int j = 0; j < faces.size(); ++j) + { + area += faces[j]->Q(); + } + area /= 3; + mass.coeffRef(i, i) = area; + } +} + + +inline void buildCotanMatrix(CMeshO &mesh, Eigen::SparseMatrix &cotanOperator){ + // 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 + + // 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 + + double cotan_l = cotan(elf, eln); + double cotan_r = cotan(ern, erf); + + // add to the matrix + cotanOperator.coeffRef(vertex_ids[vp], vertex_ids[vo]) = (cotan_l + cotan_r)/2; + + // move to the next edge + pos.FlipF();pos.FlipE(); + } while (pos != start); + } + + // compute diagonal entries + for (int i = 0; i < mesh.VN(); ++i){ + cotanOperator.coeffRef(i, i) = -cotanOperator.row(i).sum(); + } + +} + + +inline double computeAverageEdgeLength(CMeshO &mesh){ + // compute total length of all edges + double total_length = 0; + for (CMeshO::FaceIterator fi = mesh.face.begin(); fi != mesh.face.end(); ++fi) + { + vcg::Point3f p0 = fi->V(0)->P(); + vcg::Point3f p1 = fi->V(1)->P(); + vcg::Point3f p2 = fi->V(2)->P(); + double e2 = toEigen(p1 - p0).norm(); + double e1 = toEigen(p2 - p0).norm(); + double e0 = toEigen(p2 - p1).norm(); + total_length += (e0 + e1 + e2) / 2; + } + return total_length / (3./2. * mesh.FN()); +} + + +inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorXd &heat){ + Eigen::MatrixX3d heatGradientField(mesh.FN(), 3); + // 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; + } + // compute gradient of heat function at each vertex + for (int i = 0; i < mesh.FN(); ++i){ + CMeshO::FaceType *fp = &mesh.face[i]; + + vcg::Point3f p0 = fp->V(0)->P(); + vcg::Point3f p1 = fp->V(1)->P(); + vcg::Point3f p2 = fp->V(2)->P(); + + // normal unit vector + Eigen::Vector3d n = toEigen(fp->N()); + 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 + Eigen::Vector3d e0 = toEigen(p2 - p1); + e0 /= e0.norm(); + Eigen::Vector3d e1 = toEigen(p0 - p2); + e1 /= e1.norm(); + Eigen::Vector3d e2 = toEigen(p1 - p0); + e2 /= e2.norm(); + // gradient unit vectors + Eigen::Vector3d g0 = n.cross(e0); //v0 grad + Eigen::Vector3d g1 = n.cross(e1); //v1 grad + 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)]) + ) / (2 * faceArea); + + heatGradientField.row(i) = total_grad; + } + return heatGradientField; +} + + +inline Eigen::MatrixX3d normalizeVectorField(const Eigen::MatrixX3d &field){ + Eigen::MatrixX3d normalizedField(field.rows(), 3); + normalizedField.setZero(); + // normalize vector field at each vertex + for (int i = 0; i < field.rows(); ++i){ + Eigen::Vector3d v = field.row(i); + normalizedField.row(i) = v / v.norm(); + } + return normalizedField; +} + + +inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::MatrixX3d &field){ + Eigen::VectorXd divergence(mesh.VN()); + divergence.setZero(); + // initialize a hashtable from face pointers to ids + std::unordered_map face_ids; + for (int i = 0; i < mesh.FN(); ++i){ + face_ids[&mesh.face[i]] = i; + } + + // compute divergence of vector field at each vertex + for (int i = 0; i < mesh.VN(); ++i){ + CMeshO::VertexType *vp = &mesh.vert[i]; + + std::vector faces; + std::vector indices; + vcg::face::VFStarVF(vp, faces, indices); + for (int j = 0; j < faces.size(); ++j) + { + CMeshO::FaceType *fp = faces[j]; + int index = indices[j]; + vcg::Point3f p0 = fp->V(0)->P(); + vcg::Point3f p1 = fp->V(1)->P(); + vcg::Point3f p2 = fp->V(2)->P(); + // (ORDERING) edge vectors + Eigen::Vector3d el, er, eo; //left, right, opposite to vp + if (index == 0){ + el = toEigen(p2 - p0); //-e1 + er = toEigen(p1 - p0); //e2 + eo = toEigen(p2 - p1); //e0 + } else if (index == 1){ + el = toEigen(p0 - p1); //-e2 + er = toEigen(p2 - p1); //e0 + eo = toEigen(p0 - p2); //e1 + } 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); + // normalize edge vectors after cotangent computation + el /= el.norm(); + er /= er.norm(); + // add divergence contribution of given face + Eigen::Vector3d x = field.row(face_ids[fp]); + divergence(i) += (cotl * er.dot(x) + cotr * el.dot(x)) / 2; + } + } + return divergence; +} + + +#endif From dc9385c2461fcc711a5d1e7ee89e8e31ebb7999e Mon Sep 17 00:00:00 2001 From: davidboening Date: Sat, 17 Jun 2023 22:27:51 +0200 Subject: [PATCH 02/12] added heatgeodesic to cmake --- src/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 11600523e..f13281016 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -153,6 +153,7 @@ if(NOT DEFINED MESHLAB_PLUGINS) # it may be already defined in parent directory meshlabplugins/filter_dirt meshlabplugins/filter_fractal meshlabplugins/filter_func + meshlabplugins/filter_heatgeodesic # added meshlabplugins/filter_img_patch_param meshlabplugins/filter_icp meshlabplugins/filter_io_nxs From 66a60362778de1f11acf44f0e41dacc2ce61e0b4 Mon Sep 17 00:00:00 2001 From: davidboening Date: Mon, 19 Jun 2023 11:09:10 +0200 Subject: [PATCH 03/12] version 1.0 fixed requirements, tested on various meshes. --- .../filter_heatgeodesic.cpp | 60 ++++++++++++------- 1 file changed, 38 insertions(+), 22 deletions(-) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp index ececbd4b9..406f235bd 100644 --- a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp +++ b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp @@ -37,7 +37,7 @@ QString FilterHeatGeodesicPlugin::filterName(ActionIDType filterId) const { switch(filterId) { case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - return "Compute Geodesic From Selection (Heat)"; + return QString("Compute Geodesic From Selection (Heat)"); default : assert(0); return QString(); @@ -54,7 +54,7 @@ QString FilterHeatGeodesicPlugin::pythonFilterName(ActionIDType f) const { switch(f) { case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - return "compute_heat_geodesic_from_selection"; + return "compute_approximate_geodesic_from_selection"; default : assert(0); return QString(); @@ -72,10 +72,10 @@ QString FilterHeatGeodesicPlugin::pythonFilterName(ActionIDType f) const { switch(filterId) { case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - return "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."); default : assert(0); - return "Unknown Filter"; + return QString("Unknown Filter"); } } @@ -90,7 +90,8 @@ FilterHeatGeodesicPlugin::FilterClass FilterHeatGeodesicPlugin::getClass(const Q { switch(ID(a)) { case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - return FilterPlugin::Measure; + return FilterPlugin::Measure; // Note: not sure + // return FilterGeodesic::FilterClass(FilterPlugin::VertexColoring + FilterPlugin::Quality); default : assert(0); return FilterPlugin::Generic; @@ -113,10 +114,7 @@ FilterPlugin::FilterArity FilterHeatGeodesicPlugin::filterArity(const QAction*) int FilterHeatGeodesicPlugin::getPreConditions(const QAction*) const { // NOTE: note sure about these - return MeshModel::MM_VERTCOORD | MeshModel::MM_VERTQUALITY | - MeshModel::MM_FACEFACETOPO | MeshModel::MM_VERTFACETOPO | - MeshModel::MM_FACEQUALITY | MeshModel::MM_FACENORMAL | - MeshModel::MM_VERTFLAGSELECT; + return MeshModel::MM_FACEFACETOPO | MeshModel::MM_VERTFACETOPO; } /** @@ -126,7 +124,7 @@ int FilterHeatGeodesicPlugin::getPreConditions(const QAction*) const int FilterHeatGeodesicPlugin::postCondition(const QAction*) const { // NOTE: note sure about these - return MeshModel::MM_VERTQUALITY | MeshModel::MM_FACEQUALITY; + return MeshModel::MM_VERTQUALITY | MeshModel::MM_VERTCOLOR; } /** @@ -164,20 +162,40 @@ RichParameterList FilterHeatGeodesicPlugin::initParameterList(const QAction *act */ std::map FilterHeatGeodesicPlugin::applyFilter(const QAction * action, const RichParameterList & parameters, MeshDocument &md, unsigned int& /*postConditionMask*/, vcg::CallBackPos *cb) { + MeshModel &mm=*(md.mm()); switch(ID(action)) { case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : { - CMeshO &mesh = md.mm()->cm; - cb(100*0/10, "Computing Boundary Conditions..."); - Eigen::VectorXd initialConditions(mesh.VN()); - for(int i=0; i < mesh.vert.size(); i++){ - initialConditions[i] = mesh.vert[i].IsS() ? 1 : 0; - } - cb(100*1/10, "Updating Topology and Computing Face Normals..."); + mm.updateDataMask(MeshModel::MM_FACEFACETOPO); + mm.updateDataMask(MeshModel::MM_VERTFACETOPO); + mm.updateDataMask(MeshModel::MM_FACEQUALITY); + mm.updateDataMask(MeshModel::MM_FACENORMAL); + 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); @@ -230,8 +248,10 @@ std::map FilterHeatGeodesicPlugin::applyFilter(const QAct // set geodesic distance as quality for(int i=0; i < mesh.vert.size(); i++){ - mesh.vert[i].Q() = geodesicDistance(0); + mesh.vert[i].Q() = geodesicDistance(i); } + + vcg::tri::UpdateColor::PerVertexQualityRamp(mesh); } break; default : @@ -240,8 +260,4 @@ std::map FilterHeatGeodesicPlugin::applyFilter(const QAct return std::map(); } -void computeHeatGeodesicFromSelection(MeshDocument &md, vcg::CallBackPos *cb, float m){ - -} - MESHLAB_PLUGIN_NAME_EXPORTER(FilterSamplePlugin) From a3879881cc36713ea572e03bf9a5eec9f90e5996 Mon Sep 17 00:00:00 2001 From: davidboening Date: Mon, 19 Jun 2023 11:39:10 +0200 Subject: [PATCH 04/12] enabled actions to github --- src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp index 406f235bd..9151247a8 100644 --- a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp +++ b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp @@ -196,6 +196,7 @@ std::map FilterHeatGeodesicPlugin::applyFilter(const QAct break; } + cb(100*2/10, "Building Linear System (Heatflow)..."); Eigen::SparseMatrix mass(mesh.VN(), mesh.VN()); buildMassMatrix(mesh, mass); From f8f213ea4e8a90fd80c5759cc6cd51a14a3d2746 Mon Sep 17 00:00:00 2001 From: davidboening Date: Wed, 21 Jun 2023 18:35:03 +0200 Subject: [PATCH 05/12] 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(); From 9d7664aa6ed10becae96b1dc7ee456ce98794d14 Mon Sep 17 00:00:00 2001 From: davidboening Date: Thu, 22 Jun 2023 12:28:53 +0200 Subject: [PATCH 06/12] version 1.2 two levels of caching (factorization, matrices), better encapsulation --- .../filter_heatgeodesic.cpp | 120 +++++++----------- .../filter_heatgeodesic/trimesh_heatmethod.h | 62 +++++++++ 2 files changed, 105 insertions(+), 77 deletions(-) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp index c780e7d07..1d6544786 100644 --- a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp +++ b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp @@ -184,106 +184,72 @@ std::map FilterHeatGeodesicPlugin::applyFilter(const QAct } inline void FilterHeatGeodesicPlugin::computeHeatGeodesicFromSelection(CMeshO& mesh, vcg::CallBackPos* cb, float m){ + vcg::tri::Allocator::CompactEveryVector(mesh); + // build boundary conditions - Eigen::VectorXd boundaryConditions(mesh.VN()); + Eigen::VectorXd sourcePoints(mesh.VN()); int selection_count = 0; for(int i=0; i < mesh.VN(); i++){ - boundaryConditions(i) = mesh.vert[i].IsS() ? (++selection_count, 1) : 0; + sourcePoints(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); + std::shared_ptr HMSolver; - // 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; + // 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 (!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 + if (!vcg::tri::HasPerMeshAttribute(mesh, "HeatMethodSolver")){ + + // 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 buildCotanLowerTriMatrix(mesh, cotanMatrix, cb); - avg_edge_len = computeAverageEdgeLength(mesh); + averageEdgeLength = computeAverageEdgeLength(mesh); - HeatMethodData &HMData = vcg::tri::Allocator::GetPerMeshAttribute(mesh, std::string("HeatMethodData"))(); - HMData = HeatMethodData(massMatrix, cotanMatrix, avg_edge_len); + cb(91, "Matrix Factorization..."); + HMSolver = std::make_shared(HeatMethodSolver(std::move(massMatrix), std::move(cotanMatrix), averageEdgeLength, m)); + if (HMSolver->factorization_failed()) + log("Warning: factorization has failed. The mesh is non-manifold or badly conditioned (angles ~ 0deg)"); + auto HMSolverHandle = vcg::tri::Allocator::GetPerMeshAttribute>(mesh, std::string("HeatMethodSolver")); + HMSolverHandle() = HMSolver; } else { - HeatMethodData &HMData = vcg::tri::Allocator::GetPerMeshAttribute(mesh, std::string("HeatMethodData"))(); - std::tie (massMatrix, cotanMatrix, avg_edge_len) = HMData; + // 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(91, "Rebuilding Factorization..."); + HMSolver->rebuildFactorization(m); + } } - // cholesky type solver - Eigen::SimplicialLDLT> solver; + cb(95, "Computing Geodesic Distance..."); + Eigen::VectorXd geodesicDistance = HMSolver->solve(mesh, sourcePoints); + if (HMSolver->solver_failed()) + log("Warning: solver has failed."); - // 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..."); + cb(99, "Updating Mesh..."); // set geodesic distance as quality and color - for(int i=0; i < mesh.vert.size(); i++){ + for(int i=0; i < mesh.VN(); i++){ mesh.vert[i].Q() = geodesicDistance(i); } vcg::tri::UpdateColor::PerVertexQualityRamp(mesh); diff --git a/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h b/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h index 37ebe5d2e..4de5d1901 100644 --- a/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h +++ b/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h @@ -8,6 +8,8 @@ #include +#include + // LOCAL FUNCTIONS @@ -252,4 +254,64 @@ 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) + { + // initialize pointers + solver1 = std::shared_ptr>>(new Eigen::SimplicialLDLT>); + solver2 = std::shared_ptr>>(new Eigen::SimplicialLDLT>); + + // build factorization + 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; + } + 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) + Eigen::MatrixX3d unitVectorField = normalizeVectorField(-heatGradient); // (FN, 3) + Eigen::VectorXd divergence = computeVertexDivergence(mesh, unitVectorField); // (VN) + Eigen::VectorXd geodesicDistance = solver2->solve(divergence); // (VN) + if((solver1->info() != Eigen::Success) || (solver2->info() != Eigen::Success)) + solverFailed = true; + else + solverFailed = false; + // shift to impose dist(d) = 0 \forall d \in sourcePoints + 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; +}; + #endif From 4e4e950cd79731a28b5859b5679938e24a41be68 Mon Sep 17 00:00:00 2001 From: davidboening Date: Fri, 23 Jun 2023 18:53:40 +0200 Subject: [PATCH 07/12] version 1.3 fixed bottleneck (cotan weights) using Eigen::Triplet changed progress bar changed hashmap to vcg::tri::Index other minor changes --- .../filter_heatgeodesic.cpp | 17 ++-- .../filter_heatgeodesic/trimesh_heatmethod.h | 92 +++++++++---------- 2 files changed, 56 insertions(+), 53 deletions(-) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp index 1d6544786..a2da10cf0 100644 --- a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp +++ b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp @@ -186,7 +186,8 @@ std::map FilterHeatGeodesicPlugin::applyFilter(const QAct inline void FilterHeatGeodesicPlugin::computeHeatGeodesicFromSelection(CMeshO& mesh, vcg::CallBackPos* cb, float m){ vcg::tri::Allocator::CompactEveryVector(mesh); - // build boundary conditions + cb(1, "Checking Selection..."); + // build and check source vertices Eigen::VectorXd sourcePoints(mesh.VN()); int selection_count = 0; for(int i=0; i < mesh.VN(); i++){ @@ -208,7 +209,7 @@ inline void FilterHeatGeodesicPlugin::computeHeatGeodesicFromSelection(CMeshO& m // recover state if it exists 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); @@ -222,27 +223,29 @@ inline void FilterHeatGeodesicPlugin::computeHeatGeodesicFromSelection(CMeshO& m 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); + cb(20, "Building Cotan Matrix..."); + buildCotanLowerTriMatrix(mesh, cotanMatrix); averageEdgeLength = computeAverageEdgeLength(mesh); - cb(91, "Matrix Factorization..."); + cb(70, "Matrix Factorization..."); HMSolver = std::make_shared(HeatMethodSolver(std::move(massMatrix), std::move(cotanMatrix), averageEdgeLength, m)); if (HMSolver->factorization_failed()) - log("Warning: factorization has failed. The mesh is non-manifold or badly conditioned (angles ~ 0deg)"); + log("Warning: factorization has failed. The mesh is non-manifold or badly conditioned (e.g. angles ~ 0deg) or has disconnected components"); auto HMSolverHandle = vcg::tri::Allocator::GetPerMeshAttribute>(mesh, std::string("HeatMethodSolver")); HMSolverHandle() = HMSolver; } else { + cb(10, "Recovering 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(91, "Rebuilding Factorization..."); + cb(70, "Rebuilding Factorization..."); HMSolver->rebuildFactorization(m); } } - cb(95, "Computing Geodesic Distance..."); + cb(85, "Computing Geodesic Distance..."); Eigen::VectorXd geodesicDistance = HMSolver->solve(mesh, sourcePoints); if (HMSolver->solver_failed()) log("Warning: solver has failed."); diff --git a/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h b/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h index 4de5d1901..0887341a4 100644 --- a/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h +++ b/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h @@ -27,10 +27,18 @@ namespace { }; } -// GLOBAL FUNCTIONS +// simple struct to track progress bar status +struct ProgressBarCallbackT { + vcg::CallBackPos* callback; + int current_progress; + int allocated_progress; +}; + inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix &mass){ mass.resize(mesh.VN(), mesh.VN()); + mass.reserve(Eigen::VectorXi::Constant(mesh.VN(),1)); + // compute area of all faces for (CMeshO::FaceIterator fi = mesh.face.begin(); fi != mesh.face.end(); ++fi) { @@ -60,29 +68,22 @@ inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix &mass){ area += faces[j]->Q(); } area /= 3; - mass.coeffRef(i, i) = area; + mass.insert(i, i) = area; } + + mass.makeCompressed(); } -inline void buildCotanLowerTriMatrix(CMeshO &mesh, Eigen::SparseMatrix &cotanOperator, vcg::CallBackPos* cb = NULL){ +inline void buildCotanLowerTriMatrix(CMeshO &mesh, Eigen::SparseMatrix &cotanOperator){ 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; - } - // progress bar variables - int update_size = mesh.FN() / 90; // 90 updates (90% weight of the progress bar) - int progress = 0; + typedef Eigen::Triplet T; + std::vector tripletList; + tripletList.reserve(3*mesh.VN() + 3*mesh.FN()); // 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]; CMeshO::VertexType* v0 = fi->V(0); @@ -102,28 +103,31 @@ inline void buildCotanLowerTriMatrix(CMeshO &mesh, Eigen::SparseMatrix & 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]; + int i0 = vcg::tri::Index(mesh,v0); + int i1 = vcg::tri::Index(mesh,v1); + int i2 = vcg::tri::Index(mesh,v2); // save only lower triangular part if (i0 > i1) - cotanOperator.coeffRef(i0, i1) += alpha2; + tripletList.push_back(T(i0,i1,alpha2)); else - cotanOperator.coeffRef(i1, i0) += alpha2; + tripletList.push_back(T(i1,i0,alpha2)); if (i0 > i2) - cotanOperator.coeffRef(i0, i2) += alpha1; + tripletList.push_back(T(i0,i2,alpha1)); else - cotanOperator.coeffRef(i2, i0) += alpha1; + tripletList.push_back(T(i2,i0,alpha1)); if (i1 > i2) - cotanOperator.coeffRef(i1, i2) += alpha0; + tripletList.push_back(T(i1,i2,alpha0)); else - cotanOperator.coeffRef(i2, i1) += alpha0; + tripletList.push_back(T(i2,i1,alpha0)); - cotanOperator.coeffRef(i0, i0) -= (alpha1 + alpha2); - cotanOperator.coeffRef(i1, i1) -= (alpha0 + alpha2); - cotanOperator.coeffRef(i2, i2) -= (alpha0 + alpha1); + tripletList.push_back(T(i0,i0,-(alpha1 + alpha2))); + tripletList.push_back(T(i1,i1,-(alpha0 + alpha2))); + tripletList.push_back(T(i2,i2,-(alpha0 + alpha1))); } + + cotanOperator.setFromTriplets(tripletList.begin(), tripletList.end()); + cotanOperator.makeCompressed(); } @@ -146,18 +150,21 @@ inline double computeAverageEdgeLength(CMeshO &mesh){ inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorXd &heat){ Eigen::MatrixX3d heatGradientField(mesh.FN(), 3); - // 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; - } // compute gradient of heat function at each vertex for (int i = 0; i < mesh.FN(); ++i){ CMeshO::FaceType *fp = &mesh.face[i]; - vcg::Point3f p0 = fp->V(0)->P(); - vcg::Point3f p1 = fp->V(1)->P(); - vcg::Point3f p2 = fp->V(2)->P(); + CMeshO::VertexType* v0 = fp->V(0); + CMeshO::VertexType* v1 = fp->V(1); + CMeshO::VertexType* v2 = fp->V(2); + + vcg::Point3f p0 = v0->P(); + vcg::Point3f p1 = v1->P(); + vcg::Point3f p2 = v2->P(); + + int i0 = vcg::tri::Index(mesh,v0); + int i1 = vcg::tri::Index(mesh,v1); + int i2 = vcg::tri::Index(mesh,v2); // normal unit vector Eigen::Vector3d n = toEigen(fp->N()); @@ -177,11 +184,7 @@ inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorX Eigen::Vector3d g2 = n.cross(e2); //v2 grad // add vertex gradient contributions - 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); + Eigen::Vector3d tri_grad = (g0 * heat(i0) + g1 * heat(i1) + g2 * heat(i2)) / (2 * faceArea); heatGradientField.row(i) = tri_grad; } @@ -204,11 +207,6 @@ inline Eigen::MatrixX3d normalizeVectorField(const Eigen::MatrixX3d &field){ inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::MatrixX3d &field){ Eigen::VectorXd divergence(mesh.VN()); divergence.setZero(); - // initialize a hashtable from face pointers to ids - std::unordered_map face_ids; - for (int i = 0; i < mesh.FN(); ++i){ - face_ids[&mesh.face[i]] = i; - } // compute divergence of vector field at each vertex for (int i = 0; i < mesh.VN(); ++i){ @@ -221,9 +219,11 @@ inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::Matrix { CMeshO::FaceType *fp = faces[j]; int index = indices[j]; + vcg::Point3f p0 = fp->V(0)->P(); vcg::Point3f p1 = fp->V(1)->P(); vcg::Point3f p2 = fp->V(2)->P(); + // (ORDERING) edge vectors Eigen::Vector3d el, er, eo; //left, right, opposite to vp if (index == 0){ @@ -246,7 +246,7 @@ inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::Matrix el /= el.norm(); er /= er.norm(); // add divergence contribution of given face - Eigen::Vector3d x = field.row(face_ids[fp]); + Eigen::Vector3d x = field.row(vcg::tri::Index(mesh,fp)); divergence(i) += (cotl * er.dot(x) + cotr * el.dot(x)) / 2; } } From f14dca369360d0ade59eeb2bda42887876827811 Mon Sep 17 00:00:00 2001 From: davidboening Date: Wed, 12 Jul 2023 09:39:22 +0200 Subject: [PATCH 08/12] 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; }; From 3faad7bd73adcedbf0b95bea2100cbb64fca9b3f Mon Sep 17 00:00:00 2001 From: davidboening Date: Thu, 20 Jul 2023 23:01:18 +0200 Subject: [PATCH 09/12] candidate final - removed heatgeodesic filter - merged heatgeodesic filter into geodesic filter - moved algorithmic part to vcglib --- src/CMakeLists.txt | 3 +- .../filter_geodesic/filter_geodesic.cpp | 93 ++++- .../filter_geodesic/filter_geodesic.h | 4 +- .../filter_heatgeodesic/CMakeLists.txt | 5 - .../filter_heatgeodesic.cpp | 249 ------------- .../filter_heatgeodesic/filter_heatgeodesic.h | 39 --- .../filter_heatgeodesic/trimesh_heatmethod.h | 329 ------------------ 7 files changed, 81 insertions(+), 641 deletions(-) delete mode 100644 src/meshlabplugins/filter_heatgeodesic/CMakeLists.txt delete mode 100644 src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp delete mode 100644 src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h delete mode 100644 src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f13281016..87d461de4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -152,8 +152,7 @@ if(NOT DEFINED MESHLAB_PLUGINS) # it may be already defined in parent directory meshlabplugins/filter_developability meshlabplugins/filter_dirt meshlabplugins/filter_fractal - meshlabplugins/filter_func - meshlabplugins/filter_heatgeodesic # added + meshlabplugins/filter_func meshlabplugins/filter_img_patch_param meshlabplugins/filter_icp meshlabplugins/filter_io_nxs diff --git a/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp b/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp index a5149e09b..34493709b 100644 --- a/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp +++ b/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp @@ -44,7 +44,8 @@ FilterGeodesic::FilterGeodesic() typeList = { FP_QUALITY_BORDER_GEODESIC, FP_QUALITY_POINT_GEODESIC, - FP_QUALITY_SELECTED_GEODESIC + FP_QUALITY_SELECTED_GEODESIC, + FP_QUALITY_SELECTED_GEODESIC_HEAT }; for(ActionIDType tt : types()) @@ -67,6 +68,8 @@ QString FilterGeodesic::filterName(ActionIDType filter) const return QString("Colorize by geodesic distance from a given point"); case FP_QUALITY_SELECTED_GEODESIC: return QString("Colorize by geodesic distance from the selected points"); + case FP_QUALITY_SELECTED_GEODESIC_HEAT: + return QString("Colorize by approximated geodesic distance from the selected points"); default: assert(0); return QString(); } } @@ -79,6 +82,8 @@ QString FilterGeodesic::pythonFilterName(ActionIDType f) const return QString("compute_scalar_by_geodesic_distance_from_given_point_per_vertex"); case FP_QUALITY_SELECTED_GEODESIC: return QString("compute_scalar_by_geodesic_distance_from_selection_per_vertex"); + case FP_QUALITY_SELECTED_GEODESIC_HEAT: + return QString("compute_scalar_by_heat_geodesic_distance_from_selection_per_vertex"); default: assert(0); return QString(); } } @@ -87,9 +92,10 @@ QString FilterGeodesic::filterInfo(ActionIDType filterId) const { switch(filterId) { - case FP_QUALITY_BORDER_GEODESIC : return tr("Store in the quality field the geodesic distance from borders and color the mesh accordingly."); - case FP_QUALITY_POINT_GEODESIC : return tr("Store in the quality field the geodesic distance from a given point on the mesh surface and color the mesh accordingly."); - case FP_QUALITY_SELECTED_GEODESIC : return tr("Store in the quality field the geodesic distance from the selected points on the mesh surface and color the mesh accordingly."); + case FP_QUALITY_BORDER_GEODESIC : return tr("Store in the quality field the geodesic distance from borders and color the mesh accordingly."); + case FP_QUALITY_POINT_GEODESIC : return tr("Store in the quality field the geodesic distance from a given point on the mesh surface and color the mesh accordingly."); + case FP_QUALITY_SELECTED_GEODESIC : return tr("Store in the quality field the geodesic distance from the selected points on the mesh surface and color the mesh accordingly."); + case FP_QUALITY_SELECTED_GEODESIC_HEAT : return tr("Store in the quality field the approximated geodesic distance, computed via heat method, from the selected points on the mesh surface and color the mesh accordingly. This method is very sensitive to triangulation and currently does not support disconnected components."); default : assert(0); } return QString("error!"); @@ -99,9 +105,10 @@ FilterGeodesic::FilterClass FilterGeodesic::getClass(const QAction *a) const { switch(ID(a)) { - case FP_QUALITY_BORDER_GEODESIC : - case FP_QUALITY_SELECTED_GEODESIC : - case FP_QUALITY_POINT_GEODESIC : return FilterGeodesic::FilterClass(FilterPlugin::VertexColoring + FilterPlugin::Quality); + case FP_QUALITY_BORDER_GEODESIC : + case FP_QUALITY_SELECTED_GEODESIC : + case FP_QUALITY_POINT_GEODESIC : return FilterGeodesic::FilterClass(FilterPlugin::VertexColoring + FilterPlugin::Quality); + case FP_QUALITY_SELECTED_GEODESIC_HEAT : return FilterGeodesic::FilterClass(FilterPlugin::VertexColoring + FilterPlugin::Quality); default : assert(0); } return FilterPlugin::Generic; @@ -111,9 +118,10 @@ int FilterGeodesic::getRequirements(const QAction *action) { switch(ID(action)) { - case FP_QUALITY_BORDER_GEODESIC : - case FP_QUALITY_SELECTED_GEODESIC: - case FP_QUALITY_POINT_GEODESIC : return MeshModel::MM_VERTFACETOPO; + case FP_QUALITY_BORDER_GEODESIC : + case FP_QUALITY_SELECTED_GEODESIC : + case FP_QUALITY_POINT_GEODESIC : return MeshModel::MM_VERTFACETOPO; + case FP_QUALITY_SELECTED_GEODESIC_HEAT : return MeshModel::MM_VERTFACETOPO + MeshModel::MM_FACEFACETOPO; default: assert(0); } return 0; @@ -236,7 +244,56 @@ std::map FilterGeodesic::applyFilter(const QAction *filte else log("Warning: no vertices are selected! aborting geodesic computation."); } - break; + break; + case FP_QUALITY_SELECTED_GEODESIC_HEAT: + { + m.updateDataMask(MeshModel::MM_FACEFACETOPO); + m.updateDataMask(MeshModel::MM_VERTFACETOPO); + m.updateDataMask(MeshModel::MM_FACEQUALITY); + m.updateDataMask(MeshModel::MM_FACENORMAL); + m.updateDataMask(MeshModel::MM_VERTQUALITY); + m.updateDataMask(MeshModel::MM_VERTCOLOR); + + std::vector seedVec; + ForEachVertex(m.cm, [&seedVec] (CMeshO::VertexType & v) { + if (v.IsS()) seedVec.push_back(&v); + }); + + float param_m = par.getFloat("m"); + + if (seedVec.size() > 0){ + // cache factorizations to reduce runtime on successive computations + std::pair cache; + if (!vcg::tri::HasPerMeshAttribute(m.cm, "GeodesicHeatCache")){ + cache = std::make_pair(param_m, tri::GeodesicHeat::BuildCache(m.cm, param_m)); + // save cache for next compute + auto GeodesicHeatCacheHandle = tri::Allocator::GetPerMeshAttribute>(m.cm, std::string("GeodesicHeatCache")); + GeodesicHeatCacheHandle() = cache; + } + else { + // recover cache + auto GeodesicHeatCacheHandle = vcg::tri::Allocator::GetPerMeshAttribute>(m.cm, std::string("GeodesicHeatCache")); + // if m has changed rebuild everything + if (std::get<0>(GeodesicHeatCacheHandle()) != param_m){ + GeodesicHeatCacheHandle() = std::make_pair(param_m, tri::GeodesicHeat::BuildCache(m.cm, param_m)); + } + cache = GeodesicHeatCacheHandle(); + } + + if (tri::GeodesicHeat::ComputeFromCache(m.cm, seedVec, std::get<1>(cache))){ + tri::UpdateColor::PerVertexQualityRamp(m.cm); + } + else{ + log("Warning: heat method has failed. The mesh is most likely badly conditioned (e.g. angles ~ 0deg) or has disconnected components"); + // delete cache as its most likely useless after failure + auto GeodesicHeatCacheHandle = vcg::tri::Allocator::GetPerMeshAttribute>(m.cm, std::string("GeodesicHeatCache")); + tri::Allocator::DeletePerMeshAttribute>(m.cm, GeodesicHeatCacheHandle); + } + } + else + log("Warning: no vertices are selected! aborting geodesic computation."); + } + break; default: wrongActionCalled(filter); break; @@ -255,7 +312,10 @@ RichParameterList FilterGeodesic::initParameterList(const QAction *action, const break; case FP_QUALITY_SELECTED_GEODESIC : parlst.addParam(RichPercentage("maxDistance",m.cm.bbox.Diag(),0,m.cm.bbox.Diag()*2,"Max Distance","If not zero it indicates a cut off value to be used during geodesic distance computation.")); - break; + break; + case FP_QUALITY_SELECTED_GEODESIC_HEAT : + parlst.addParam(RichFloat("m", 1, "m", "Multiplier used in backward Euler timestep.")); + break; default: break; // do not add any parameter for the other filters } return parlst; @@ -265,10 +325,11 @@ int FilterGeodesic::postCondition(const QAction * filter) const { switch (ID(filter)) { - case FP_QUALITY_BORDER_GEODESIC : - case FP_QUALITY_SELECTED_GEODESIC : - case FP_QUALITY_POINT_GEODESIC : return MeshModel::MM_VERTCOLOR + MeshModel::MM_VERTQUALITY; - default : return MeshModel::MM_ALL; + case FP_QUALITY_BORDER_GEODESIC : + case FP_QUALITY_SELECTED_GEODESIC : + case FP_QUALITY_POINT_GEODESIC : return MeshModel::MM_VERTCOLOR + MeshModel::MM_VERTQUALITY; + case FP_QUALITY_SELECTED_GEODESIC_HEAT : return MeshModel::MM_VERTCOLOR + MeshModel::MM_VERTQUALITY; + default : return MeshModel::MM_ALL; } } MESHLAB_PLUGIN_NAME_EXPORTER(FilterGeodesic) diff --git a/src/meshlabplugins/filter_geodesic/filter_geodesic.h b/src/meshlabplugins/filter_geodesic/filter_geodesic.h index 6fe6d24c2..498b43aba 100644 --- a/src/meshlabplugins/filter_geodesic/filter_geodesic.h +++ b/src/meshlabplugins/filter_geodesic/filter_geodesic.h @@ -26,6 +26,7 @@ #include #include #include +#include class FilterGeodesic : public QObject, public FilterPlugin @@ -42,7 +43,8 @@ class FilterGeodesic : public QObject, public FilterPlugin enum { FP_QUALITY_BORDER_GEODESIC, FP_QUALITY_POINT_GEODESIC, - FP_QUALITY_SELECTED_GEODESIC + FP_QUALITY_SELECTED_GEODESIC, + FP_QUALITY_SELECTED_GEODESIC_HEAT } ; diff --git a/src/meshlabplugins/filter_heatgeodesic/CMakeLists.txt b/src/meshlabplugins/filter_heatgeodesic/CMakeLists.txt deleted file mode 100644 index 732bcaacc..000000000 --- a/src/meshlabplugins/filter_heatgeodesic/CMakeLists.txt +++ /dev/null @@ -1,5 +0,0 @@ -set(SOURCES filter_heatgeodesic.cpp) - -set(HEADERS filter_heatgeodesic.h trimesh_heatmethod.h) - -add_meshlab_plugin(filter_heatgeodesic ${SOURCES} ${HEADERS}) diff --git a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp b/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp deleted file mode 100644 index ec435b462..000000000 --- a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp +++ /dev/null @@ -1,249 +0,0 @@ -#include "trimesh_heatmethod.h" - -#include "filter_heatgeodesic.h" - -#include - -/** - * @brief Constructor usually performs only two simple tasks of filling the two lists - * - typeList: with all the possible id of the filtering actions - * - actionList with the corresponding actions. If you want to add icons to - * your filtering actions you can do here by construction the QActions accordingly - */ -FilterHeatGeodesicPlugin::FilterHeatGeodesicPlugin() -{ - typeList = {FP_COMPUTE_HEATGEODESIC_FROM_SELECTION}; - - for(ActionIDType tt : types()) - actionList.push_back(new QAction(filterName(tt), this)); -} - -FilterHeatGeodesicPlugin::~FilterHeatGeodesicPlugin() -{ -} - -QString FilterHeatGeodesicPlugin::pluginName() const -{ - return "FilterHeatGeodesic"; -} - -/** - * @brief ST() must return the very short string describing each filtering action - * (this string is used also to define the menu entry) - * @param filterId: the id of the filter - * @return the name of the filter - */ -QString FilterHeatGeodesicPlugin::filterName(ActionIDType filterId) const -{ - switch(filterId) { - case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - return QString("Compute Geodesic From Selection (Heat)"); - default : - assert(0); - return QString(); - } -} - -/** - * @brief FilterHeatGeodesicPlugin::pythonFilterName if you want that your filter should have a different - * name on pymeshlab, use this function to return its python name. - * @param f - * @return - */ -QString FilterHeatGeodesicPlugin::pythonFilterName(ActionIDType f) const -{ - switch(f) { - case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - return "compute_approximate_geodesic_from_selection"; - default : - assert(0); - return QString(); - } -} - - -/** - * @brief // Info() must return the longer string describing each filtering action - * (this string is used in the About plugin dialog) - * @param filterId: the id of the filter - * @return an info string of the filter - */ - QString FilterHeatGeodesicPlugin::filterInfo(ActionIDType filterId) 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. " - "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"); - } -} - - /** - * @brief The FilterClass describes in which generic class of filters it fits. - * This choice affect the submenu in which each filter will be placed - * More than a single class can be chosen. - * @param a: the action of the filter - * @return the class od the filter - */ -FilterHeatGeodesicPlugin::FilterClass FilterHeatGeodesicPlugin::getClass(const QAction *a) const -{ - switch(ID(a)) { - case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - return FilterPlugin::Measure; // Note: not sure - // return FilterGeodesic::FilterClass(FilterPlugin::VertexColoring + FilterPlugin::Quality); - default : - assert(0); - return FilterPlugin::Generic; - } -} - -/** - * @brief FilterSamplePlugin::filterArity - * @return - */ -FilterPlugin::FilterArity FilterHeatGeodesicPlugin::filterArity(const QAction*) const -{ - return SINGLE_MESH; -} - -/** - * @brief FilterSamplePlugin::getPreConditions - * @return - */ -int FilterHeatGeodesicPlugin::getPreConditions(const QAction*) const -{ - // NOTE: note sure about these - return MeshModel::MM_FACEFACETOPO | MeshModel::MM_VERTFACETOPO; -} - -/** - * @brief FilterSamplePlugin::postCondition - * @return - */ -int FilterHeatGeodesicPlugin::postCondition(const QAction*) const -{ - // NOTE: note sure about these - return MeshModel::MM_VERTQUALITY | MeshModel::MM_VERTCOLOR; -} - -/** - * @brief This function define the needed parameters for each filter. Return true if the filter has some parameters - * it is called every time, so you can set the default value of parameters according to the mesh - * For each parameter you need to define, - * - the name of the parameter, - * - the default value - * - the string shown in the dialog - * - a possibly long string describing the meaning of that parameter (shown as a popup help in the dialog) - * @param action - * @param m - * @param parlst - */ -RichParameterList FilterHeatGeodesicPlugin::initParameterList(const QAction *action,const MeshModel &m) -{ - RichParameterList parlst; - switch(ID(action)) { - case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - parlst.addParam(RichFloat("m", 1, "m", "Multiplier used in backward Euler timestep.")); - break; - default : - assert(0); - } - return parlst; -} - -/** - * @brief The Real Core Function doing the actual mesh processing. - * @param action - * @param md: an object containing all the meshes and rasters of MeshLab - * @param par: the set of parameters of each filter - * @param cb: callback object to tell MeshLab the percentage of execution of the filter - * @return true if the filter has been applied correctly, false otherwise - */ -std::map FilterHeatGeodesicPlugin::applyFilter(const QAction * action, const RichParameterList & parameters, MeshDocument &md, unsigned int& /*postConditionMask*/, vcg::CallBackPos *cb) -{ - MeshModel &mm=*(md.mm()); - switch(ID(action)) { - case FP_COMPUTE_HEATGEODESIC_FROM_SELECTION : - { - mm.updateDataMask(MeshModel::MM_FACEFACETOPO); - mm.updateDataMask(MeshModel::MM_VERTFACETOPO); - mm.updateDataMask(MeshModel::MM_FACEQUALITY); - mm.updateDataMask(MeshModel::MM_FACENORMAL); - mm.updateDataMask(MeshModel::MM_VERTQUALITY); - mm.updateDataMask(MeshModel::MM_VERTCOLOR); - CMeshO &mesh = mm.cm; - // TODO: compact all vectors - computeHeatGeodesicFromSelection(mesh, cb, parameters.getFloat("m")); - } - break; - default : - wrongActionCalled(action); - } - return std::map(); -} - -inline void FilterHeatGeodesicPlugin::computeHeatGeodesicFromSelection(CMeshO& mesh, vcg::CallBackPos* cb, float m){ - vcg::tri::Allocator::CompactEveryVector(mesh); - - cb(1, "Checking Selection..."); - // build and check source vertices - Eigen::VectorXd sourcePoints(mesh.VN()); - int selection_count = 0; - for(int i=0; i < mesh.VN(); i++){ - sourcePoints(i) = mesh.vert[i].IsS() ? (++selection_count, 1) : 0; - } - if (selection_count < 1){ - log("Warning: no vertices are selected! aborting computation."); - 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; - // if no state can be recovered, initialize it - if (!vcg::tri::HasPerMeshAttribute(mesh, "HeatMethodSolver")){ - 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 disconnected components"); - auto HMSolverHandle = vcg::tri::Allocator::GetPerMeshAttribute>(mesh, std::string("HeatMethodSolver")); - HMSolverHandle() = HMSolver; - } - else { - cb(5, "Attempting to Recover State..."); - // recover solver - 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..."); - Eigen::VectorXd geodesicDistance = HMSolver->solve(mesh, sourcePoints); - if (HMSolver->solver_failed()) - log("Warning: solver has failed."); - - cb(99, "Updating Mesh..."); - // set geodesic distance as quality and color - for(int i=0; i < mesh.VN(); 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 deleted file mode 100644 index b421550ab..000000000 --- a/src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h +++ /dev/null @@ -1,39 +0,0 @@ -#ifndef FILTERSAMPLE_PLUGIN_H -#define FILTERSAMPLE_PLUGIN_H - -#include - -class FilterHeatGeodesicPlugin : public QObject, public FilterPlugin -{ - Q_OBJECT - MESHLAB_PLUGIN_IID_EXPORTER(FILTER_PLUGIN_IID) - Q_INTERFACES(FilterPlugin) - -public: - enum { FP_COMPUTE_HEATGEODESIC_FROM_SELECTION } ; - - FilterHeatGeodesicPlugin(); - virtual ~FilterHeatGeodesicPlugin(); - - QString pluginName() const; - - QString filterName(ActionIDType filter) const; - QString pythonFilterName(ActionIDType f) const; - QString filterInfo(ActionIDType filter) const; - FilterClass getClass(const QAction* a) const; - FilterArity filterArity(const QAction*) const; - int getPreConditions(const QAction *) const; - int postCondition(const QAction* ) const; - RichParameterList initParameterList(const QAction*, const MeshModel &/*m*/); - std::map applyFilter( - const QAction* action, - const RichParameterList & parameters, - MeshDocument &md, - unsigned int& postConditionMask, - vcg::CallBackPos * cb); - -private: - 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 deleted file mode 100644 index 348c628ca..000000000 --- a/src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h +++ /dev/null @@ -1,329 +0,0 @@ -#ifndef TRIMESH_HEAT_METHOD -#define TRIMESH_HEAT_METHOD - -#include -#include - -#include - -#include - -#include - - -// LOCAL FUNCTIONS - -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(); - }; -} - -// simple struct to track progress bar status -struct ProgressBarCallbackT { - vcg::CallBackPos* callback; - int current_progress; - int allocated_progress; -}; - - -inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix &mass){ - mass.resize(mesh.VN(), mesh.VN()); - mass.reserve(Eigen::VectorXi::Constant(mesh.VN(),1)); - - // compute area of all faces - for (CMeshO::FaceIterator fi = mesh.face.begin(); fi != mesh.face.end(); ++fi) - { - vcg::Point3f p0 = fi->V(0)->P(); - vcg::Point3f p1 = fi->V(1)->P(); - vcg::Point3f p2 = fi->V(2)->P(); - double e0 = toEigen(p1 - p0).norm(); - double e1 = toEigen(p2 - p0).norm(); - double e2 = toEigen(p2 - p1).norm(); - double s = (e0 + e1 + e2) / 2; - double area = std::sqrt(s * (s - e0) * (s - e1) * (s - e2)); - // we store face areas in quality field to avoid a hash table - // this will also be useful for the gradient computation - fi->Q() = area; - } - // compute area of the dual cell for each vertex - for (int i = 0; i < mesh.VN(); ++i){ - CMeshO::VertexType *vp = &mesh.vert[i]; - - std::vector faces; - std::vector indices; - vcg::face::VFStarVF(vp, faces, indices); - - double area = 0; - for (int j = 0; j < faces.size(); ++j) - { - area += faces[j]->Q(); - } - area /= 3; - mass.insert(i, i) = area; - } - - mass.makeCompressed(); -} - - -inline void buildCotanLowerTriMatrix(CMeshO &mesh, Eigen::SparseMatrix &cotanOperator){ - cotanOperator.resize(mesh.VN(), mesh.VN()); - - typedef Eigen::Triplet T; - std::vector tripletList; - tripletList.reserve(3*mesh.VN() + 3*mesh.FN()); - - // compute cotan weights - for (int i = 0; i < mesh.FN(); ++i){ - CMeshO::FaceType* fi = &mesh.face[i]; - - CMeshO::VertexType* v0 = fi->V(0); - CMeshO::VertexType* v1 = fi->V(1); - CMeshO::VertexType* v2 = fi->V(2); - - vcg::Point3f p0 = v0->P(); - vcg::Point3f p1 = v1->P(); - vcg::Point3f p2 = v2->P(); - - 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 = vcg::tri::Index(mesh,v0); - int i1 = vcg::tri::Index(mesh,v1); - int i2 = vcg::tri::Index(mesh,v2); - - // save only lower triangular part - if (i0 > i1) - tripletList.push_back(T(i0,i1,alpha2)); - else - tripletList.push_back(T(i1,i0,alpha2)); - if (i0 > i2) - tripletList.push_back(T(i0,i2,alpha1)); - else - tripletList.push_back(T(i2,i0,alpha1)); - if (i1 > i2) - tripletList.push_back(T(i1,i2,alpha0)); - else - tripletList.push_back(T(i2,i1,alpha0)); - - tripletList.push_back(T(i0,i0,-(alpha1 + alpha2))); - tripletList.push_back(T(i1,i1,-(alpha0 + alpha2))); - tripletList.push_back(T(i2,i2,-(alpha0 + alpha1))); - } - - cotanOperator.setFromTriplets(tripletList.begin(), tripletList.end()); - cotanOperator.makeCompressed(); -} - - -inline double computeAverageEdgeLength(CMeshO &mesh){ - // compute total length of all edges - double total_length = 0; - for (CMeshO::FaceIterator fi = mesh.face.begin(); fi != mesh.face.end(); ++fi) - { - vcg::Point3f p0 = fi->V(0)->P(); - vcg::Point3f p1 = fi->V(1)->P(); - vcg::Point3f p2 = fi->V(2)->P(); - double e2 = toEigen(p1 - p0).norm(); - double e1 = toEigen(p2 - p0).norm(); - double e0 = toEigen(p2 - p1).norm(); - total_length += (e0 + e1 + e2) / 2; - } - return total_length / (3./2. * mesh.FN()); -} - - -inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorXd &heat){ - Eigen::MatrixX3d heatGradientField(mesh.FN(), 3); - // compute gradient of heat function at each vertex - for (int i = 0; i < mesh.FN(); ++i){ - CMeshO::FaceType *fp = &mesh.face[i]; - - CMeshO::VertexType* v0 = fp->V(0); - CMeshO::VertexType* v1 = fp->V(1); - CMeshO::VertexType* v2 = fp->V(2); - - vcg::Point3f p0 = v0->P(); - vcg::Point3f p1 = v1->P(); - vcg::Point3f p2 = v2->P(); - - int i0 = vcg::tri::Index(mesh,v0); - int i1 = vcg::tri::Index(mesh,v1); - int i2 = vcg::tri::Index(mesh,v2); - - // normal unit vector - Eigen::Vector3d n = toEigen(fp->N()); - n /= n.norm(); - // face area - double faceArea = fp->Q(); - // edge unit vectors (counter-clockwise) - Eigen::Vector3d e0 = toEigen(p2 - p1); - e0 /= e0.norm(); - Eigen::Vector3d e1 = toEigen(p0 - p2); - e1 /= e1.norm(); - Eigen::Vector3d e2 = toEigen(p1 - p0); - e2 /= e2.norm(); - // gradient unit vectors - Eigen::Vector3d g0 = n.cross(e0); //v0 grad - Eigen::Vector3d g1 = n.cross(e1); //v1 grad - Eigen::Vector3d g2 = n.cross(e2); //v2 grad - - // add vertex gradient contributions - Eigen::Vector3d tri_grad = (g0 * heat(i0) + g1 * heat(i1) + g2 * heat(i2)) / (2 * faceArea); - - heatGradientField.row(i) = tri_grad; - } - return heatGradientField; -} - - -inline Eigen::MatrixX3d normalizeVectorField(const Eigen::MatrixX3d &field){ - Eigen::MatrixX3d normalizedField(field.rows(), 3); - normalizedField.setZero(); - // normalize vector field at each vertex - for (int i = 0; i < field.rows(); ++i){ - Eigen::Vector3d v = field.row(i); - normalizedField.row(i) = v / v.norm(); - } - return normalizedField; -} - - -inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::MatrixX3d &field){ - Eigen::VectorXd divergence(mesh.VN()); - divergence.setZero(); - - // compute divergence of vector field at each vertex - for (int i = 0; i < mesh.VN(); ++i){ - CMeshO::VertexType *vp = &mesh.vert[i]; - - std::vector faces; - std::vector indices; - vcg::face::VFStarVF(vp, faces, indices); - for (int j = 0; j < faces.size(); ++j) - { - CMeshO::FaceType *fp = faces[j]; - int index = indices[j]; - - vcg::Point3f p0 = fp->V(0)->P(); - vcg::Point3f p1 = fp->V(1)->P(); - vcg::Point3f p2 = fp->V(2)->P(); - - // (ORDERING) edge vectors - Eigen::Vector3d el, er, eo; //left, right, opposite to vp - if (index == 0){ - el = toEigen(p2 - p0); //-e1 - er = toEigen(p1 - p0); //e2 - eo = toEigen(p2 - p1); //e0 - } else if (index == 1){ - el = toEigen(p0 - p1); //-e2 - er = toEigen(p2 - p1); //e0 - eo = toEigen(p0 - p2); //e1 - } 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); - double cotr = cotan(-er, eo); - // normalize edge vectors after cotangent computation - el /= el.norm(); - er /= er.norm(); - // add divergence contribution of given face - Eigen::Vector3d x = field.row(vcg::tri::Index(mesh,fp)); - divergence(i) += (cotl * er.dot(x) + cotr * el.dot(x)) / 2; - } - } - return divergence; -} - - -class HeatMethodSolver { -public: - 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; - - Eigen::VectorXd solve(CMeshO &mesh, Eigen::VectorXd &sourcePoints){ - Eigen::VectorXd heatflow = solver1->solve(sourcePoints); // (VN) - Eigen::MatrixX3d heatGradient = computeVertexGradient(mesh, heatflow); // (FN, 3) - Eigen::MatrixX3d unitVectorField = normalizeVectorField(-heatGradient); // (FN, 3) - Eigen::VectorXd divergence = computeVertexDivergence(mesh, unitVectorField); // (VN) - Eigen::VectorXd geodesicDistance = solver2->solve(divergence); // (VN) - if((solver1->info() != Eigen::Success) || (solver2->info() != Eigen::Success)) - solverFailed = true; - else - solverFailed = false; - // 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; - double m; - bool factorizationFailed; - bool solverFailed; -}; - -#endif From 7d9b26dfe8c84cbd675c5e6cf6efc6b91eee3057 Mon Sep 17 00:00:00 2001 From: davidboening Date: Fri, 21 Jul 2023 10:40:05 +0200 Subject: [PATCH 10/12] Updated descriptions, added progress bar. --- .../filter_geodesic/filter_geodesic.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp b/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp index 34493709b..5a73e102d 100644 --- a/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp +++ b/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp @@ -95,7 +95,7 @@ QString FilterGeodesic::filterInfo(ActionIDType filterId) const case FP_QUALITY_BORDER_GEODESIC : return tr("Store in the quality field the geodesic distance from borders and color the mesh accordingly."); case FP_QUALITY_POINT_GEODESIC : return tr("Store in the quality field the geodesic distance from a given point on the mesh surface and color the mesh accordingly."); case FP_QUALITY_SELECTED_GEODESIC : return tr("Store in the quality field the geodesic distance from the selected points on the mesh surface and color the mesh accordingly."); - case FP_QUALITY_SELECTED_GEODESIC_HEAT : return tr("Store in the quality field the approximated geodesic distance, computed via heat method, from the selected points on the mesh surface and color the mesh accordingly. This method is very sensitive to triangulation and currently does not support disconnected components."); + case FP_QUALITY_SELECTED_GEODESIC_HEAT : return tr("Store in the quality field the approximated geodesic distance, computed via heat method (Crane et al.), from the selected points on the mesh surface and color the mesh accordingly. As this implementation does not use intrinsic triangulation it is very sensitive to trinagulation. First run takes longer as factorization has to be build. "); default : assert(0); } return QString("error!"); @@ -107,7 +107,7 @@ FilterGeodesic::FilterClass FilterGeodesic::getClass(const QAction *a) const { case FP_QUALITY_BORDER_GEODESIC : case FP_QUALITY_SELECTED_GEODESIC : - case FP_QUALITY_POINT_GEODESIC : return FilterGeodesic::FilterClass(FilterPlugin::VertexColoring + FilterPlugin::Quality); + case FP_QUALITY_POINT_GEODESIC : case FP_QUALITY_SELECTED_GEODESIC_HEAT : return FilterGeodesic::FilterClass(FilterPlugin::VertexColoring + FilterPlugin::Quality); default : assert(0); } @@ -127,7 +127,7 @@ int FilterGeodesic::getRequirements(const QAction *action) return 0; } -std::map FilterGeodesic::applyFilter(const QAction *filter, const RichParameterList & par, MeshDocument &md, unsigned int& /*postConditionMask*/, vcg::CallBackPos * /*cb*/) +std::map FilterGeodesic::applyFilter(const QAction *filter, const RichParameterList & par, MeshDocument &md, unsigned int& /*postConditionMask*/, vcg::CallBackPos *cb) { MeshModel &m=*(md.mm()); CMeshO::VertexIterator vi; @@ -265,12 +265,14 @@ std::map FilterGeodesic::applyFilter(const QAction *filte // cache factorizations to reduce runtime on successive computations std::pair cache; if (!vcg::tri::HasPerMeshAttribute(m.cm, "GeodesicHeatCache")){ + cb(10, "Building Cache: Computing Factorizations."); cache = std::make_pair(param_m, tri::GeodesicHeat::BuildCache(m.cm, param_m)); // save cache for next compute auto GeodesicHeatCacheHandle = tri::Allocator::GetPerMeshAttribute>(m.cm, std::string("GeodesicHeatCache")); GeodesicHeatCacheHandle() = cache; } else { + cb(70, "Recovering Cache."); // recover cache auto GeodesicHeatCacheHandle = vcg::tri::Allocator::GetPerMeshAttribute>(m.cm, std::string("GeodesicHeatCache")); // if m has changed rebuild everything @@ -279,6 +281,7 @@ std::map FilterGeodesic::applyFilter(const QAction *filte } cache = GeodesicHeatCacheHandle(); } + cb(80, "Computing Geodesic Distance."); if (tri::GeodesicHeat::ComputeFromCache(m.cm, seedVec, std::get<1>(cache))){ tri::UpdateColor::PerVertexQualityRamp(m.cm); @@ -314,7 +317,7 @@ RichParameterList FilterGeodesic::initParameterList(const QAction *action, const parlst.addParam(RichPercentage("maxDistance",m.cm.bbox.Diag(),0,m.cm.bbox.Diag()*2,"Max Distance","If not zero it indicates a cut off value to be used during geodesic distance computation.")); break; case FP_QUALITY_SELECTED_GEODESIC_HEAT : - parlst.addParam(RichFloat("m", 1, "m", "Multiplier used in backward Euler timestep.")); + parlst.addParam(RichFloat("m", 1.0, tr("Euler Step"), tr("Multiplier used in backward Euler timestep. Changing this value will reset the cache."))); break; default: break; // do not add any parameter for the other filters } @@ -328,7 +331,7 @@ int FilterGeodesic::postCondition(const QAction * filter) const case FP_QUALITY_BORDER_GEODESIC : case FP_QUALITY_SELECTED_GEODESIC : case FP_QUALITY_POINT_GEODESIC : return MeshModel::MM_VERTCOLOR + MeshModel::MM_VERTQUALITY; - case FP_QUALITY_SELECTED_GEODESIC_HEAT : return MeshModel::MM_VERTCOLOR + MeshModel::MM_VERTQUALITY; + case FP_QUALITY_SELECTED_GEODESIC_HEAT : return MeshModel::MM_VERTCOLOR + MeshModel::MM_VERTQUALITY + MeshModel::MM_FACEQUALITY; default : return MeshModel::MM_ALL; } } From ae6212cb8f910c121e3bc8f1d8064ba49fdb44ff Mon Sep 17 00:00:00 2001 From: davidboening Date: Fri, 21 Jul 2023 19:33:00 +0200 Subject: [PATCH 11/12] reflecting vcglib changes --- .../filter_geodesic/filter_geodesic.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp b/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp index 5a73e102d..dd600f699 100644 --- a/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp +++ b/src/meshlabplugins/filter_geodesic/filter_geodesic.cpp @@ -263,25 +263,26 @@ std::map FilterGeodesic::applyFilter(const QAction *filte if (seedVec.size() > 0){ // cache factorizations to reduce runtime on successive computations - std::pair cache; + std::pair::GeodesicHeatCache> cache; if (!vcg::tri::HasPerMeshAttribute(m.cm, "GeodesicHeatCache")){ - cb(10, "Building Cache: Computing Factorizations."); + cb(20, "Building Cache: Computing Factorizations..."); cache = std::make_pair(param_m, tri::GeodesicHeat::BuildCache(m.cm, param_m)); // save cache for next compute - auto GeodesicHeatCacheHandle = tri::Allocator::GetPerMeshAttribute>(m.cm, std::string("GeodesicHeatCache")); + auto GeodesicHeatCacheHandle = tri::Allocator::GetPerMeshAttribute::GeodesicHeatCache>>(m.cm, std::string("GeodesicHeatCache")); GeodesicHeatCacheHandle() = cache; } else { - cb(70, "Recovering Cache."); + cb(10, "Recovering Cache..."); // recover cache - auto GeodesicHeatCacheHandle = vcg::tri::Allocator::GetPerMeshAttribute>(m.cm, std::string("GeodesicHeatCache")); + auto GeodesicHeatCacheHandle = vcg::tri::Allocator::GetPerMeshAttribute::GeodesicHeatCache>>(m.cm, std::string("GeodesicHeatCache")); // if m has changed rebuild everything if (std::get<0>(GeodesicHeatCacheHandle()) != param_m){ + cb(20, "Parameter Changed: Rebuilding Factorizations..."); GeodesicHeatCacheHandle() = std::make_pair(param_m, tri::GeodesicHeat::BuildCache(m.cm, param_m)); } cache = GeodesicHeatCacheHandle(); } - cb(80, "Computing Geodesic Distance."); + cb(80, "Computing Geodesic Distance..."); if (tri::GeodesicHeat::ComputeFromCache(m.cm, seedVec, std::get<1>(cache))){ tri::UpdateColor::PerVertexQualityRamp(m.cm); @@ -289,8 +290,8 @@ std::map FilterGeodesic::applyFilter(const QAction *filte else{ log("Warning: heat method has failed. The mesh is most likely badly conditioned (e.g. angles ~ 0deg) or has disconnected components"); // delete cache as its most likely useless after failure - auto GeodesicHeatCacheHandle = vcg::tri::Allocator::GetPerMeshAttribute>(m.cm, std::string("GeodesicHeatCache")); - tri::Allocator::DeletePerMeshAttribute>(m.cm, GeodesicHeatCacheHandle); + auto GeodesicHeatCacheHandle = vcg::tri::Allocator::GetPerMeshAttribute::GeodesicHeatCache>>(m.cm, std::string("GeodesicHeatCache")); + tri::Allocator::DeletePerMeshAttribute::GeodesicHeatCache>>(m.cm, GeodesicHeatCacheHandle); } } else From d7166cfdf5c54824f81759a17a15f033535c2c37 Mon Sep 17 00:00:00 2001 From: Alessandro Muntoni Date: Mon, 31 Jul 2023 10:43:30 +0200 Subject: [PATCH 12/12] update vcg submodule --- src/vcglib | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/vcglib b/src/vcglib index 87edc09b9..80b1d21bd 160000 --- a/src/vcglib +++ b/src/vcglib @@ -1 +1 @@ -Subproject commit 87edc09b9077b9e9b5fc474544ea76b17b665399 +Subproject commit 80b1d21bdab5031c2ffdba971ff8b0b7c7a76911