[Meng Wei] final version

src/angle_defect.cpp

@@ -1,9 +1,21 @@
 #include "../include/angle_defect.h"
+#include "../include/internal_angles.h"
+#include <igl/squared_edge_lengths.h>
 
 void angle_defect(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::VectorXd & D)
 {
-  D = Eigen::VectorXd::Zero(V.rows());
+  // Find internal_angles
+  Eigen::MatrixXd l_sqr, A;
+  igl::squared_edge_lengths(V, F, l_sqr);
+  internal_angles(l_sqr, A);
+
+  D = Eigen::VectorXd::Ones(V.rows()) * 3.141592653589793238463 * 2;
+  for(int i = 0; i < F.rows(); i ++) {
+    for(int j = 0; j < 3; j ++) {
+      D(F(i, j)) -= A(i, j);
+    }
+  }
 }

src/internal_angles.cpp

@@ -4,5 +4,23 @@ void internal_angles(
   const Eigen::MatrixXd & l_sqr,
   Eigen::MatrixXd & A)
 {
-  // Add with your code
+  A.resize(l_sqr.rows(), 3);
+  // Law of Cosines:
+  for (int i = 0; i < l_sqr.rows(); i ++) {
+    double a2 = l_sqr(i, 0);
+    double b2 = l_sqr(i, 1);
+    double c2 = l_sqr(i, 2);
+
+    double cos_C = (a2 + b2 - c2) / (2.0 * sqrt(a2) * sqrt(b2));
+    double angle_C = std::acos(cos_C);
+
+    double cos_A = (b2 + c2 - a2) / (2.0 * sqrt(b2) * sqrt(c2));
+    double angle_A = std::acos(cos_A);
+
+    double angle_B = 180 - angle_C - angle_A;
+
+    A(i, 0) = angle_A;
+    A(i, 1) = angle_B;
+    A(i, 2) = angle_C;
+  }
 }

src/mean_curvature.cpp

@@ -1,10 +1,34 @@
 #include "../include/mean_curvature.h"
+#include <igl/cotmatrix.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/massmatrix.h>
+#include <igl/invert_diag.h>
 
 void mean_curvature(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::VectorXd & H)
 {
-  // Replace with your code
-  H = Eigen::VectorXd::Zero(V.rows());
+  H.resize(V.rows());
+  // Compute cot laplacian and mass matrix
+  Eigen::SparseMatrix<double> L, M, inverse_M;
+  igl::cotmatrix(V, F, L);
+  igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+  igl::invert_diag(M, inverse_M);
+
+  // HN = M-1 LV
+  Eigen::MatrixXd HN = inverse_M * L * V;
+
+  // Compute normal vectors
+  Eigen::MatrixXd N;
+  igl::per_vertex_normals(V, F, N);
+
+  // Compute H and its sign
+  for(int i = 0; i < HN.rows(); i ++) {
+    H(i) = HN.row(i).norm();
+    if (HN.row(i).dot( N.row(i) ) < 0) {
+      // Non-consistency
+      H(i) *= -1;
+    }
+  }
 }

src/principal_curvatures.cpp

@@ -1,4 +1,9 @@
 #include "../include/principal_curvatures.h"
+#include <igl/pinv.h>
+#include <igl/slice.h>
+#include <igl/adjacency_matrix.h>
+#include <set>
+#include <Eigen/Eigen>
 
 void principal_curvatures(
   const Eigen::MatrixXd & V,
@@ -13,4 +18,78 @@ void principal_curvatures(
   K2 = Eigen::VectorXd::Zero(V.rows());
   D1 = Eigen::MatrixXd::Zero(V.rows(),3);
   D2 = Eigen::MatrixXd::Zero(V.rows(),3);
+
+  // Grab adjacent vertex
+  Eigen::SparseMatrix<int> adj_matrix;
+  igl::adjacency_matrix(F, adj_matrix);
+
+  for (int i = 0; i < V.rows(); i ++) {
+    // Find adjacent vertices
+    std::set<int> adj_ver;
+    // Idea inspired by Eigen Sparse Matrix documentation
+    for (Eigen::SparseMatrix<int>::InnerIterator it1(adj_matrix, i); it1; ++it1) {
+      int level_1 = it1.row();
+      adj_ver.insert(level_1);
+      for (Eigen::SparseMatrix<int>::InnerIterator it2(adj_matrix, i); it2; ++it2) {
+        int level_2 = it2.row();
+        if (level_2 != i ) {
+          adj_ver.insert(level_2);
+        }
+      }
+    }
+
+    // Find P: (vi - v)
+    Eigen::MatrixXd P(adj_ver.size(), 3);
+    int j = 0;
+    for(auto cur_vi : adj_ver) {
+      P.row(j) = V.row(cur_vi) - V.row(i);
+      j ++;
+    }
+
+    // PCA: by issue #18 and documentation
+    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(P.transpose()*P);
+    // Since R3, two most:
+    Eigen::VectorXd u = es.eigenvectors().col(1);
+    Eigen::VectorXd v = es.eigenvectors().col(2);
+    // one leasT:
+    Eigen::VectorXd w = P * es.eigenvectors().col(0);
+
+    // Construct coefficients
+    Eigen::MatrixXd coeff(P.rows(), 5);
+    coeff.col(0) = P * u;
+    coeff.col(1) = P * v;
+    coeff.col(2) = (P * u).array().square();
+    coeff.col(3) = (P * u).cwiseProduct(P * v);
+    coeff.col(4) = (P * v).array().square();
+
+    // Compute As
+    Eigen::MatrixXd inverse_coeff;
+    igl::pinv(coeff, inverse_coeff);
+    Eigen::VectorXd A = inverse_coeff * w;
+
+    // Compute variables
+    double E = 1 + A(0) * A(0);
+    double F = A(0) * A(1);
+    double G = 1 + A(1) * A(1);
+    double common = std::sqrt(E + G - 1);
+    double e = (2 * A(2)) / common;
+    double f = A(3) / common;
+    double g = (2 * A(4)) / common;
+
+    // Construct S
+    Eigen::Matrix2d left, Right, S;
+    left << e, f,
+            f, g;
+    Right << E, F,
+             F, G;
+    S = - left * Right.inverse();
+
+    // PCA on S
+    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es_S(S);
+    K1(i) = es_S.eigenvalues()(0);
+    K2(i) = es_S.eigenvalues()(1);
+    D1.row(i) = es_S.eigenvectors()(0, 0) * u + es_S.eigenvectors()(0, 1) * v;
+    D2.row(i) = es_S.eigenvectors()(1, 0) * u + es_S.eigenvectors()(0, 1) * v;
+
+  }
 }