Merge aee5ed8 into 48face2

pmp-library · Jan 28, 2020 · f9c8036 · f9c8036
2 parents 48face2 + aee5ed8
commit f9c8036
Show file tree

Hide file tree

Showing 13 changed files with 315 additions and 139 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -15,6 +15,7 @@ This project aims to adhere to [Semantic Versioning](https://semver.org/spec/v2.
 ### Added
 
 - Add constructors using initializer lists to Matrix/Vector classes
+- Add assignment from and cast from Eigen matrices and vectors
 
 ### Fixed
 

diff --git a/docs/bibliography.bib b/docs/bibliography.bib
@@ -247,3 +247,11 @@ @article{kazhdan_2012
 	journal = {Computer Graphics Forum}
 }
 
+@article{alexa_2011_laplace,
+ author = {Alexa, Marc and Wardetzky, Max},
+ title = {Discrete Laplacians on General Polygonal Meshes},
+ year = {2011},
+ volume = {30},
+ number = {4},
+ journal = {ACM Transactions on Graphics}
+}
diff --git a/external/pmp-data b/external/pmp-data
diff --git a/src/pmp/MatVec.h b/src/pmp/MatVec.h
@@ -15,6 +15,7 @@
 #include <assert.h>
 #include <limits>
 #include <initializer_list>
+#include <Eigen/Dense>
 
 //=============================================================================
 
@@ -166,6 +167,38 @@ class Matrix
             data_[i] = static_cast<Scalar>(m[i]);
     }
 
+    //! assign from Eigen
+    template <typename Derived>
+    Matrix<Scalar, M, N>& operator=(const Eigen::MatrixBase<Derived>& m)
+    {
+        // don't distinguish between row and column vectors
+        if (m.rows()==1 || m.cols()==1)
+        {
+            assert( m.size()==size() );
+            for (int i = 0; i < size(); ++i)
+                (*this)[i] = m[i];
+        }
+        else
+        {
+            assert(m.rows()==rows() && m.cols()==cols());
+            for (int i = 0; i < rows(); ++i)
+                for (int j = 0; j < cols(); ++j)
+                    (*this)(i,j) = m(i,j);
+        }
+        return *this;
+    }
+
+    //! cast to Eigen
+    template <typename OtherScalar>
+    operator Eigen::Matrix<OtherScalar, M, N>() const
+    {
+        Eigen::Matrix<OtherScalar, M, N> m;
+        for (int i = 0; i < rows(); ++i)
+            for (int j = 0; j < cols(); ++j)
+                m(i,j) = static_cast<OtherScalar>((*this)(i,j));
+        return m;
+    }
+
     //! return identity matrix (only for square matrices, N==M)
     static Matrix<Scalar, M, N> identity();
 
@@ -204,7 +237,12 @@ class Matrix
     Scalar* data() { return data_; }
 
     //! normalize matrix/vector by dividing through Frobenius/Euclidean norm
-    void normalize() { *this /= norm(*this); }
+    void normalize() 
+    { 
+        Scalar n = norm(*this);
+        n = (n > std::numeric_limits<Scalar>::min()) ? 1.0 / n : 0.0;
+        *this *= n;
+    }
 
     //! divide matrix by scalar
     Matrix<Scalar, M, N>& operator/=(const Scalar s)

diff --git a/src/pmp/algorithms/SurfaceFairing.cpp b/src/pmp/algorithms/SurfaceFairing.cpp
@@ -140,15 +140,14 @@ void SurfaceFairing::fair(unsigned int k)
     const unsigned int n = vertices.size();
     SparseMatrix A(n, n);
     Eigen::MatrixXd B(n, 3);
+    dvec3 b;
 
     std::map<Vertex, double> row;
     std::vector<Triplet> triplets;
 
     for (unsigned int i = 0; i < n; ++i)
     {
-        B(i, 0) = 0.0;
-        B(i, 1) = 0.0;
-        B(i, 2) = 0.0;
+        b = dvec3(0.0);
 
         setup_matrix_row(vertices[i], vweight_, eweight_, k, row);
 
@@ -163,11 +162,11 @@ void SurfaceFairing::fair(unsigned int k)
             }
             else
             {
-                B(i, 0) -= w * points_[v][0];
-                B(i, 1) -= w * points_[v][1];
-                B(i, 2) -= w * points_[v][2];
+                b -= w * static_cast<dvec3>(points_[v]);
             }
         }
+
+        B.row(i) = (Eigen::Vector3d) b;
     }
 
     A.setFromTriplets(triplets.begin(), triplets.end());
@@ -183,10 +182,7 @@ void SurfaceFairing::fair(unsigned int k)
     else
     {
         for (unsigned int i = 0; i < n; ++i)
-        {
-            auto v = vertices[i];
-            points_[v] = Point(X(idx_[v], 0), X(idx_[v], 1), X(idx_[v], 2));
-        }
+            points_[vertices[i]] = X.row(i);
     }
 }
 

diff --git a/src/pmp/algorithms/SurfaceHoleFilling.cpp b/src/pmp/algorithms/SurfaceHoleFilling.cpp
@@ -456,9 +456,7 @@ void SurfaceHoleFilling::relaxation()
         else
             triplets.emplace_back(i, idx[v], c);
 
-        B(i, 0) = b[0];
-        B(i, 1) = b[1];
-        B(i, 2) = b[2];
+        B.row(i) = (Eigen::Vector3d) b;
     }
 
     // solve least squares system
@@ -477,10 +475,7 @@ void SurfaceHoleFilling::relaxation()
     // copy solution to mesh vertices
     for (int i = 0; i < n; ++i)
     {
-        Vertex v = vertices[i];
-        points_[v][0] = X(i, 0);
-        points_[v][1] = X(i, 1);
-        points_[v][2] = X(i, 2);
+        points_[vertices[i]] = X.row(i);
     }
 
     // clean up

diff --git a/src/pmp/algorithms/SurfaceNormals.cpp b/src/pmp/algorithms/SurfaceNormals.cpp
@@ -15,28 +15,64 @@ namespace pmp {
 
 //=============================================================================
 
+Normal SurfaceNormals::compute_face_normal(const SurfaceMesh& mesh, Face f)
+{
+    Halfedge h = mesh.halfedge(f);
+    Halfedge hend = h;
+
+    auto vpoint = mesh.get_vertex_property<Point>("v:point");
+
+    Point p0 = vpoint[mesh.to_vertex(h)];
+    h = mesh.next_halfedge(h);
+    Point p1 = vpoint[mesh.to_vertex(h)];
+    h = mesh.next_halfedge(h);
+    Point p2 = vpoint[mesh.to_vertex(h)];
+
+    if (mesh.next_halfedge(h) == hend) // face is a triangle
+    {
+        return normalize(cross(p2 -= p1, p0 -= p1));
+    }
+    else // face is a general polygon
+    {
+        Normal n(0, 0, 0);
+
+        // this computes the sum of cross products (area-weighted normals)
+        // of triangles generated by inserting the centroid c:
+        //   sum_i (p_{i} - c) x (p_{i+1} - c)
+        // The point c cancels out, leading to
+        //   sum_i (p_{i} x p_{i+1}
+        // This vector then has to be normalized.
+        for (auto h : mesh.halfedges(f))
+        {
+            n += cross(vpoint[mesh.from_vertex(h)], vpoint[mesh.to_vertex(h)]);
+        }
+
+        return normalize(n);
+    }
+}
+
+//-----------------------------------------------------------------------------
+
 Normal SurfaceNormals::compute_vertex_normal(const SurfaceMesh& mesh, Vertex v)
 {
     Point nn(0, 0, 0);
-    Halfedge h = mesh.halfedge(v);
 
-    if (h.is_valid())
+    if (!mesh.is_isolated(v))
     {
         auto vpoint = mesh.get_vertex_property<Point>("v:point");
-
-        const Halfedge hend = h;
         const Point p0 = vpoint[v];
 
-        Point n, p1, p2;
+        Normal n;
+        Point p1, p2;
         Scalar cosine, angle, denom;
+        bool is_triangle;
 
-        do
+        for (auto h: mesh.halfedges(v))
         {
             if (!mesh.is_boundary(h))
             {
                 p1 = vpoint[mesh.to_vertex(h)];
                 p1 -= p0;
-
                 p2 = vpoint[mesh.from_vertex(mesh.prev_halfedge(h))];
                 p2 -= p0;
 
@@ -51,20 +87,17 @@ Normal SurfaceNormals::compute_vertex_normal(const SurfaceMesh& mesh, Vertex v)
                         cosine = 1.0;
                     angle = acos(cosine);
 
-                    n = cross(p1, p2);
+                    // compute triangle or polygon normal
+                    is_triangle = (mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h))) == h);
+                    n = is_triangle ? 
+                        normalize(cross(p1,p2)) :
+                        compute_face_normal(mesh, mesh.face(h));
 
-                    // check whether normal is != 0
-                    denom = norm(n);
-                    if (denom > std::numeric_limits<Scalar>::min())
-                    {
-                        n *= angle / denom;
-                        nn += n;
-                    }
+                    n *= angle;
+                    nn += n;
                 }
             }
-
-            h = mesh.cw_rotated_halfedge(h);
-        } while (h != hend);
+        }
 
         nn = normalize(nn);
     }
@@ -74,43 +107,6 @@ Normal SurfaceNormals::compute_vertex_normal(const SurfaceMesh& mesh, Vertex v)
 
 //-----------------------------------------------------------------------------
 
-Normal SurfaceNormals::compute_face_normal(const SurfaceMesh& mesh, Face f)
-{
-    Halfedge h = mesh.halfedge(f);
-    Halfedge hend = h;
-
-    auto vpoint = mesh.get_vertex_property<Point>("v:point");
-
-    Point p0 = vpoint[mesh.to_vertex(h)];
-    h = mesh.next_halfedge(h);
-    Point p1 = vpoint[mesh.to_vertex(h)];
-    h = mesh.next_halfedge(h);
-    Point p2 = vpoint[mesh.to_vertex(h)];
-
-    if (mesh.next_halfedge(h) == hend) // face is a triangle
-    {
-        return normalize(cross(p2 -= p1, p0 -= p1));
-    }
-    else // face is a general polygon
-    {
-        Normal n(0, 0, 0);
-
-        // this computes the sum of cross products (area-weighted normals)
-        // of triangles generated by inserting the centroid c:
-        //   sum_i (p_{i} - c) x (p_{i+1} - c)
-        // The point c cancels out, leading to
-        //   sum_i (p_{i} x p_{i+1}
-        // This vector then has to be normalized.
-        for (auto h : mesh.halfedges(f))
-        {
-            n += cross(vpoint[mesh.from_vertex(h)], vpoint[mesh.to_vertex(h)]);
-        }
-
-        return normalize(n);
-    }
-}
-
-//-----------------------------------------------------------------------------
 
 Normal SurfaceNormals::compute_corner_normal(const SurfaceMesh& mesh,
                                              Halfedge h, Scalar crease_angle)
@@ -138,13 +134,10 @@ Normal SurfaceNormals::compute_corner_normal(const SurfaceMesh& mesh,
 
         Point n, p1, p2;
         Scalar cosine, angle, denom;
+        bool is_triangle;
 
         // compute normal of h's face
-        p1 = vpoint[mesh.to_vertex(mesh.next_halfedge(h))];
-        p1 -= p0;
-        p2 = vpoint[mesh.from_vertex(h)];
-        p2 -= p0;
-        const Point nf = normalize(cross(p1, p2));
+        const Point nf = compute_face_normal(mesh, mesh.face(h));
 
         // average over all incident faces
         do
@@ -156,31 +149,29 @@ Normal SurfaceNormals::compute_corner_normal(const SurfaceMesh& mesh,
                 p2 = vpoint[mesh.from_vertex(h)];
                 p2 -= p0;
 
-                n = cross(p1, p2);
+                // compute triangle or polygon normal
+                is_triangle = (mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h))) == h);
+                n = is_triangle ? 
+                    normalize(cross(p1,p2)) :
+                    compute_face_normal(mesh, mesh.face(h));
 
-                // check whether normal is != 0
-                denom = norm(n);
-                if (denom > std::numeric_limits<Scalar>::min())
+                // check whether normal is withing crease_angle bound
+                if (dot(n, nf) >= cos_crease_angle)
                 {
-                    n /= denom;
 
-                    // check whether normal is withing crease_angle bound
-                    if (dot(n, nf) >= cos_crease_angle)
+                    // check whether we can robustly compute angle
+                    denom = sqrt(dot(p1, p1) * dot(p2, p2));
+                    if (denom > std::numeric_limits<Scalar>::min())
                     {
-                        // check whether we can robustly compute angle
-                        denom = sqrt(dot(p1, p1) * dot(p2, p2));
-                        if (denom > std::numeric_limits<Scalar>::min())
-                        {
-                            cosine = dot(p1, p2) / denom;
-                            if (cosine < -1.0)
-                                cosine = -1.0;
-                            else if (cosine > 1.0)
-                                cosine = 1.0;
-                            angle = acos(cosine);
-
-                            n *= angle;
-                            nn += n;
-                        }
+                        cosine = dot(p1, p2) / denom;
+                        if (cosine < -1.0)
+                            cosine = -1.0;
+                        else if (cosine > 1.0)
+                            cosine = 1.0;
+                        angle = acos(cosine);
+
+                        n *= angle;
+                        nn += n;
                     }
                 }
             }

diff --git a/src/pmp/algorithms/SurfaceNormals.h b/src/pmp/algorithms/SurfaceNormals.h
@@ -54,6 +54,9 @@ class SurfaceNormals
     static Normal compute_vertex_normal(const SurfaceMesh& mesh, Vertex v);
 
     //! \brief Compute the normal vector of face \c f.
+    //! \details Normal is computed as (normalized) sum of per-corner
+    //! cross products of the two incident edges. This corresponds to
+    //! the normalized vector area in \cite alexa_2011_laplace
     static Normal compute_face_normal(const SurfaceMesh& mesh, Face f);
 
     //! \brief Compute the normal vector of the polygon corner specified by the