From 8ef1e76548e7e5cb7b869239b4e5e7f85ab4ded0 Mon Sep 17 00:00:00 2001 From: davidboening Date: Sat, 17 Jun 2023 22:27:09 +0200 Subject: [PATCH] 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