HW06 Liupei Lu

src/angle_defect.cpp

@@ -1,9 +1,24 @@
 #include "../include/angle_defect.h"
+#include "../include/internal_angles.h"
+#include <igl/squared_edge_lengths.h>
+#include <cmath>
 
 void angle_defect(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::VectorXd & D)
+    const Eigen::MatrixXd &V,
+    const Eigen::MatrixXi &F,
+    Eigen::VectorXd &D)
 {
-  D = Eigen::VectorXd::Zero(V.rows());
+  //D = Eigen::VectorXd::Zero(V.rows());
+  Eigen::MatrixXd L, A;
+  igl::squared_edge_lengths(V, F, L);
+  internal_angles(L, A);
+
+  //2pi - A
+  D.resize(V.rows());
+  for (int i = 0; i < F.rows(); i++)
+  {
+    D(F(i, 0)) = 2.00 * M_PI - A(i, 0);
+    D(F(i, 1)) = 2.00 * M_PI - A(i, 1);
+    D(F(i, 2)) = 2.00 * M_PI - A(i, 2);
+  }
 }

src/internal_angles.cpp

@@ -1,8 +1,22 @@
 #include "../include/internal_angles.h"
+#include <cmath>
 
 void internal_angles(
-  const Eigen::MatrixXd & l_sqr,
-  Eigen::MatrixXd & A)
+    const Eigen::MatrixXd &l_sqr,
+    Eigen::MatrixXd &A)
 {
+  A.resize(l_sqr.rows(), 3);
   // Add with your code
+  for (int i = 0; i < l_sqr.rows(); i++)
+  {
+    //get edge squares
+    double a = l_sqr(i, 0); //1,2
+    double b = l_sqr(i, 1); //2,0
+    double c = l_sqr(i, 2); //0,1
+
+    //by cosine law
+    A(i, 0) = acos((b + c - a) / (2.00 * sqrt(b) * sqrt(c)));
+    A(i, 1) = acos((a + c - b) / (2.00 * sqrt(a) * sqrt(c)));
+    A(i, 2) = acos((a + b - c) / (2.00 * sqrt(a) * sqrt(b)));
+  }
 }

src/mean_curvature.cpp

@@ -1,10 +1,38 @@
 #include "../include/mean_curvature.h"
+#include <igl/cotmatrix.h>
+#include <igl/invert_diag.h>
+#include <igl/massmatrix.h>
+#include <igl/per_vertex_normals.h>
 
 void mean_curvature(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::VectorXd & H)
+    const Eigen::MatrixXd &V,
+    const Eigen::MatrixXi &F,
+    Eigen::VectorXd &H)
 {
-  // Replace with your code
-  H = Eigen::VectorXd::Zero(V.rows());
+  //H = Eigen::VectorXd::Zero(V.rows());
+  Eigen::SparseMatrix<double> L, M, M_in;
+
+  igl::cotmatrix(V, F, L);
+  igl::massmatrix(V, F, igl::MassMatrixType::MASSMATRIX_TYPE_VORONOI, M);
+  igl::invert_diag(M, M_in);
+
+  //calculate Hn(mean curvature normal)
+  Eigen::MatrixXd Hn, N;
+  Hn = M_in * L * V;
+
+  //get per vertex normal
+  igl::per_vertex_normals(V, F, N);
+
+  H.resize(V.rows());
+
+  //match normal
+  for (int i = 0; i < V.rows(); i++)
+  {
+    double p = Hn.row(i).dot(N.row(i));
+    if (p >= 0)
+      H(i) = Hn.row(i).norm();
+    //reverse if negative
+    else
+      H(i) = -1.00 * Hn.row(i).norm();
+  }
 }

src/principal_curvatures.cpp

@@ -1,16 +1,89 @@
 #include "../include/principal_curvatures.h"
+#include <igl/adjacency_matrix.h>
+#include <igl/sort.h>
+#include <igl/pinv.h>
+#include <Eigen/SVD>
+#include <igl/eigs.h>
 
 void principal_curvatures(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::MatrixXd & D1,
-  Eigen::MatrixXd & D2,
-  Eigen::VectorXd & K1,
-  Eigen::VectorXd & K2)
+    const Eigen::MatrixXd &V,
+    const Eigen::MatrixXi &F,
+    Eigen::MatrixXd &D1,
+    Eigen::MatrixXd &D2,
+    Eigen::VectorXd &K1,
+    Eigen::VectorXd &K2)
 {
-  // Replace with your code
-  K1 = Eigen::VectorXd::Zero(V.rows());
-  K2 = Eigen::VectorXd::Zero(V.rows());
-  D1 = Eigen::MatrixXd::Zero(V.rows(),3);
-  D2 = Eigen::MatrixXd::Zero(V.rows(),3);
+    // Replace with your code
+    K1 = Eigen::VectorXd::Zero(V.rows());
+    K2 = Eigen::VectorXd::Zero(V.rows());
+    D1 = Eigen::MatrixXd::Zero(V.rows(), 3);
+    D2 = Eigen::MatrixXd::Zero(V.rows(), 3);
+
+    Eigen::SparseMatrix<double> A, A2;
+    igl::adjacency_matrix(F, A);
+    //Compute 2 ring
+    A2 = A * A;
+
+    for (int i = 0; i < V.rows(); i++)
+    {
+        Eigen::MatrixXd P;
+        int k = A2.col(i).nonZeros();
+        P.resize(k, 3);
+        int index = 0;
+
+        //grab all other vertices to P in 2 rings
+        for (Eigen::SparseMatrix<double>::InnerIterator it(A2, i); it; ++it)
+        {
+            P.row(index) = V.row(it.row()) - V.row(i);
+            index++;
+        }
+
+        //SVD decomposition of P
+        Eigen::JacobiSVD<Eigen::MatrixXd> svd(P.transpose() * P, Eigen::ComputeThinU | Eigen::ComputeThinV);
+        Eigen::MatrixXd U = svd.matrixU();
+
+        //get two most principle directions and one least principle direction
+        Eigen::VectorXd u = U.col(0) / U.col(0).norm();
+        Eigen::VectorXd v = U.col(1) / U.col(1).norm();
+        Eigen::VectorXd w = U.col(2) / U.col(2).norm();
+
+        //maxtrix for w
+        Eigen::MatrixXd C(k, 5), Ci;
+        Eigen::VectorXd v1, v2, v3, v4, v5;
+        v1 = P * u;
+        v2 = P * v;
+        v3 = (P * u).cwiseProduct(P * u);
+        v4 = (P * u).cwiseProduct(P * v);
+        v5 = (P * v).cwiseProduct(P * v);
+        C << v1, v2, v3, v4, v5;
+
+        //solve for a1-a5
+        igl::pinv(C, Ci);
+        Eigen::VectorXd a = Ci * P * w;
+
+        //get S coefficients
+        double E, F, G, e, f, g;
+        E = 1.0 + a(0) * a(0);
+        F = a(0) * a(1);
+        G = 1 + a(1) * a(1);
+        e = (2 * a(2)) / sqrt(a(0) * a(0) + 1 + a(1) * a(1));
+        f = (a(3)) / sqrt(a(0) * a(0) + 1 + a(1) * a(1));
+        g = (2 * a(4)) / sqrt(a(0) * a(0) + 1 + a(1) * a(1));
+
+        Eigen::Matrix2d first, second;
+        first << e, f,
+                f, g;
+
+        second << E, F,
+                F, G;
+
+        //perform enigen decomposition
+        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> solver;
+        solver.compute(-first * second.inverse());
+
+        K1(i) = solver.eigenvalues()(0);
+        K2(i) = solver.eigenvalues()(1);
+        D1.row(i) = solver.eigenvectors()(0, 0) * u + solver.eigenvectors()(0, 1) * v;
+        D2.row(i) = solver.eigenvectors()(1, 0) * u + solver.eigenvectors()(1, 1) * v;
+    }
 }