added heat geodesic (not tested)

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})

src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.cpp

@@ -0,0 +1,247 @@
+#include "trimesh_heatmethod.h"
+
+#include "filter_heatgeodesic.h"
+
+#include <common/plugins/interfaces/meshlab_plugin_logger.h>
+
+/**
+ * @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<std::string, QVariant> 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<CMeshO>::VertexFace(mesh);
+        vcg::tri::UpdateTopology<CMeshO>::FaceFace(mesh);
+        vcg::tri::UpdateNormal<CMeshO>::PerFaceNormalized(mesh);
+
+        cb(100*2/10, "Building Linear System (Heatflow)...");
+        Eigen::SparseMatrix<double> mass(mesh.VN(), mesh.VN());
+        buildMassMatrix(mesh, mass);
+
+        Eigen::SparseMatrix<double> 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<double> system1(mesh.VN(), mesh.VN());
+        system1 = mass - timestep * cotanOperator;
+
+        Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> 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<double> system2(mesh.VN(), mesh.VN());
+        system2 = cotanOperator; //+ 1e-6 * Eigen::Matrix<double,-1,-1>::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<std::string, QVariant>();
+}
+
+void computeHeatGeodesicFromSelection(MeshDocument &md, vcg::CallBackPos *cb, float m){
+
+}
+
+MESHLAB_PLUGIN_NAME_EXPORTER(FilterSamplePlugin)

src/meshlabplugins/filter_heatgeodesic/filter_heatgeodesic.h

@@ -0,0 +1,39 @@
+#ifndef FILTERSAMPLE_PLUGIN_H
+#define FILTERSAMPLE_PLUGIN_H
+
+#include <common/plugins/interfaces/filter_plugin.h>
+
+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<std::string, QVariant> 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

src/meshlabplugins/filter_heatgeodesic/trimesh_heatmethod.h

@@ -0,0 +1,250 @@
+#ifndef TRIMESH_HEAT_METHOD
+#define TRIMESH_HEAT_METHOD
+
+#include <Eigen/Sparse>
+#include <Eigen/SparseCholesky>
+
+#include <vcg/complex/algorithms/update/quality.h>
+
+#include <common/plugins/interfaces/filter_plugin.h>
+
+
+// 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<double> &mass);
+inline void buildCotanMatrix(CMeshO &mesh, Eigen::SparseMatrix<double> &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<double> &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<CMeshO::FaceType*> faces;
+        std::vector<int> indices;
+        vcg::face::VFStarVF<CMeshO::FaceType>(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<double> &cotanOperator){
+    // initialize a hashtable from vertex pointers to ids
+    std::unordered_map<CMeshO::VertexType*, int> 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<CMeshO::FaceType> pos(fp, vp);
+        vcg::face::Pos<CMeshO::FaceType> 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<CMeshO::VertexType*, int> 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<CMeshO::FaceType*, int> 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<CMeshO::FaceType*> faces;
+        std::vector<int> indices;
+        vcg::face::VFStarVF<CMeshO::FaceType>(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