Tianhao Hu

src/angle_defect.cpp

@@ -1,9 +1,22 @@
 #include "../include/angle_defect.h"
-
+#include <Eigen/Sparse>
+#include <igl/squared_edge_lengths.h>
+#include "../include/internal_angles.h"
 void angle_defect(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::VectorXd & D)
 {
-  D = Eigen::VectorXd::Zero(V.rows());
+  Eigen::MatrixXd l_sqr;
+  igl::squared_edge_lengths(V, F, l_sqr);
+
+  Eigen::MatrixXd A;
+  internal_angles(l_sqr, A);
+
+  D.resize(V.rows());
+  D.setConstant(2*M_PI);
+
+  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,13 @@ void internal_angles(
   const Eigen::MatrixXd & l_sqr,
   Eigen::MatrixXd & A)
 {
-  // Add with your code
+  A = Eigen::MatrixXd(l_sqr.rows(), 3);
+  for (int m = 0; m < A.rows(); m++)
+    for (int i = 0; i < 3; i++)
+    {
+      int j = (i + 1) % 3, k = (i + 2) % 3;
+      double cosine = (l_sqr(m, i) - l_sqr(m, j) - l_sqr(m, k)) /
+        (-2 * sqrt(l_sqr(m, j)*l_sqr(m, k)));
+      A(m, i) = acos(cosine);
+    }
 }

src/mean_curvature.cpp

@@ -1,10 +1,32 @@
 #include "../include/mean_curvature.h"
+#include <Eigen/Sparse>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/invert_diag.h>
+#include <igl/per_vertex_normals.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());
+  Eigen::SparseMatrix<double> L; L.setZero();
+  igl::cotmatrix(V, F, L);
+
+  Eigen::SparseMatrix<double> M; M.setZero();
+  igl::massmatrix(V, F, igl::MassMatrixType::MASSMATRIX_TYPE_DEFAULT, M);
+  igl::invert_diag(M, M);
+
+  Eigen::MatrixXd N; N.setZero();
+  igl::per_vertex_normals(V, F, N);
+
+  Eigen::MatrixXd Hn = M*L*V;
+  H.resize(Hn.rows());
+
+  for (int i = 0; i < Hn.rows(); i++)
+  {
+    double dot = Hn.row(i).dot(N.row(i));
+    H[i] = Hn.row(i).norm()*(dot /abs(dot));
+  }
+
 }

src/principal_curvatures.cpp

@@ -1,4 +1,7 @@
 #include "../include/principal_curvatures.h"
+#include <igl/adjacency_matrix.h>
+#include <igl/pinv.h>
+#include <igl/sort.h>
 
 void principal_curvatures(
   const Eigen::MatrixXd & V,
@@ -13,4 +16,60 @@ void principal_curvatures(
   K2 = Eigen::VectorXd::Zero(V.rows());
   D1 = Eigen::MatrixXd::Zero(V.rows(),3);
   D2 = Eigen::MatrixXd::Zero(V.rows(),3);
+
+  Eigen::SparseMatrix<int> adj;
+  igl::adjacency_matrix(F, adj);
+  Eigen::SparseMatrix<int> adj2 = adj*adj;
+
+  for (int i = 0; i < V.rows(); ++i)
+  {
+    Eigen::MatrixXd P(adj2.col(i).nonZeros(), 3);
+    int p = 0;
+    for (Eigen::SparseMatrix<int>::InnerIterator it(adj2, i); it; ++it)
+    {
+      P.row(p) = V.row(it.index()) - V.row(i);
+      p++;
+    }
+
+    Eigen::JacobiSVD<Eigen::MatrixXd> svd(P, Eigen::ComputeThinU | Eigen::ComputeThinV);
+    Eigen::MatrixXd W = svd.matrixV();
+    Eigen::MatrixXd UV = P*W.leftCols(2);
+    Eigen::MatrixXd B = P*W.rightCols(1);
+
+    Eigen::MatrixXd C(B.rows(), 5);
+    for (int j = 0; j < UV.rows(); ++j)
+        C.row(j) << UV(j, 0), UV(j, 1), UV(j, 0)*UV(j, 0), UV(j, 0)*UV(j, 1), UV(j, 1)*UV(j, 1);
+
+    igl::pinv(C, C);
+    Eigen::VectorXd A = C*B;
+
+    double E = 1 + A(0)*A(0);
+    double _F = A(0)*A(1);
+    double G = 1 + A(1)*A(1);
+    double e = 2 * A(2) / sqrt(A(0)*A(0) + 1 + A(1)*A(1));
+    double f = A(3) / sqrt(A(0)*A(0) + 1 + A(1)*A(1));
+    double g = 2 * A(4) / sqrt(A(0)*A(0) + 1 + A(1)*A(1));
+
+    Eigen::MatrixXd SL(2, 2);
+    Eigen::MatrixXd SR(2, 2);
+    SL << e, f, f, g;
+    SR << E, _F, _F, G;
+    Eigen::MatrixXd S = -SL*SR.inverse();
+
+    Eigen::EigenSolver<Eigen::MatrixXd> eigensolver(S);
+    Eigen::MatrixXd eigen_vecs = eigensolver.eigenvectors().real();
+    Eigen::VectorXd eigen_vals = eigensolver.eigenvalues().real().cwiseAbs();
+    Eigen::VectorXd eigen_vals_sorted;
+    Eigen::VectorXd eigen_index_sorted;
+    igl::sort(eigen_vals, 1, false, eigen_vals_sorted, eigen_index_sorted);
+
+    Eigen::VectorXd d1 = eigen_vecs.col(eigen_index_sorted(0));
+    Eigen::VectorXd d2 = eigen_vecs.col(eigen_index_sorted(1));
+
+    K1(i) = eigen_vals_sorted(0);
+    K2(i) = eigen_vals_sorted(1);
+
+    D1.row(i) = W.leftCols(2)*d1;
+    D2.row(i) = W.leftCols(2)*d2;
+  }
 }