From 4f8fb3ef064d662380f14976ea348f5387028732 Mon Sep 17 00:00:00 2001
From: Feilong Song <feilong.song@mail.utoronto.ca>
Date: Wed, 14 Nov 2018 00:02:34 -0500
Subject: [PATCH] hw6

---
 src/angle_defect.cpp         | 14 ++++++-
 src/internal_angles.cpp      |  9 ++++
 src/mean_curvature.cpp       | 21 +++++++++-
 src/principal_curvatures.cpp | 79 +++++++++++++++++++++++++++++++++---
 4 files changed, 116 insertions(+), 7 deletions(-)

diff --git a/src/angle_defect.cpp b/src/angle_defect.cpp
index 78589fb..30c4075 100644
--- a/src/angle_defect.cpp
+++ b/src/angle_defect.cpp
@@ -1,9 +1,21 @@
 #include "../include/angle_defect.h"
+#include <igl/squared_edge_lengths.h>
+#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, A;
+  igl::squared_edge_lengths(V, F, l_sqr);
+  internal_angles(l_sqr, A);
+
+  D = Eigen::VectorXd::Constant(V.rows(), 2 * M_PI);
+
+  for (int i = 0; i < A.rows(); i++) {
+    for (int j = 0; j < 3; j++) {
+      D(F(i, j)) -= A(i, j);
+    }
+  }
 }
diff --git a/src/internal_angles.cpp b/src/internal_angles.cpp
index 3e14c75..b39d3b6 100644
--- a/src/internal_angles.cpp
+++ b/src/internal_angles.cpp
@@ -5,4 +5,13 @@ void internal_angles(
   Eigen::MatrixXd & A)
 {
   // Add with your code
+  A = Eigen::MatrixXd::Zero(l_sqr.rows(), l_sqr.cols());
+  for (int m = 0; m < A.rows(); m++) {
+    for (int i = 0; i < 3; i++) {
+      double l1 = l_sqr(m, i);
+      double l2 = l_sqr(m, (i + 1) % 3);
+      double l3 = l_sqr(m, (i + 2) % 3);
+      A(m, i) = acos((l3 + l2 - l1) / (2.0 * sqrt(l3 * l2)));
+    }
+  }
 }
diff --git a/src/mean_curvature.cpp b/src/mean_curvature.cpp
index 41fce89..57a3a2d 100644
--- a/src/mean_curvature.cpp
+++ b/src/mean_curvature.cpp
@@ -1,4 +1,9 @@
 #include "../include/mean_curvature.h"
+#include <Eigen/Sparse>
+#include <igl/massmatrix.h>
+#include <igl/cotmatrix.h>
+#include <igl/invert_diag.h>
+#include <igl/per_vertex_normals.h>
 
 void mean_curvature(
   const Eigen::MatrixXd & V,
@@ -6,5 +11,19 @@ void mean_curvature(
   Eigen::VectorXd & H)
 {
   // Replace with your code
-  H = Eigen::VectorXd::Zero(V.rows());
+  Eigen::SparseMatrix<double> L, M, Mi;
+  igl::cotmatrix(V, F, L);
+  igl::massmatrix(V, F, igl::MassMatrixType::MASSMATRIX_TYPE_DEFAULT, M);
+  igl::invert_diag(M, Mi);
+
+  Eigen::MatrixXd N;
+  igl::per_vertex_normals(V, F, N);
+
+  Eigen::MatrixXd MLV = Mi * L * V;
+
+  H.resize(MLV.rows());
+  for (int i = 0; i < H.size(); i++) {
+    int sign = (MLV.row(i).dot(N.row(i)) >= 0) ? 1 : -1;
+    H(i) = sign * MLV.row(i).norm();
+  }
 }
diff --git a/src/principal_curvatures.cpp b/src/principal_curvatures.cpp
index 6a51d2e..c24f0fb 100644
--- a/src/principal_curvatures.cpp
+++ b/src/principal_curvatures.cpp
@@ -1,4 +1,10 @@
 #include "../include/principal_curvatures.h"
+#include <Eigen/Eigenvalues>
+#include <Eigen/Sparse>
+#include <Eigen/Geometry>
+#include <igl/adjacency_matrix.h>
+#include <set>
+#include <igl/pinv.h>
 
 void principal_curvatures(
   const Eigen::MatrixXd & V,
@@ -6,11 +12,74 @@ void principal_curvatures(
   Eigen::MatrixXd & D1,
   Eigen::MatrixXd & D2,
   Eigen::VectorXd & K1,
-  Eigen::VectorXd & K2)
-{
+  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);
-}
+  D1 = Eigen::MatrixXd::Zero(V.rows(), 3);
+  D2 = Eigen::MatrixXd::Zero(V.rows(), 3);
+
+  Eigen::SparseMatrix<int> adj;
+  igl::adjacency_matrix(F, adj);
+
+  for (int i = 0; i < V.rows(); i++) {
+    // Gather two ring vertices
+    std::set<int> twoRing;
+    for(Eigen::SparseMatrix<int>::InnerIterator it(adj, i); it; ++it) {
+      int oneRing = it.row();
+      twoRing.insert(oneRing);
+      for(Eigen::SparseMatrix<int>::InnerIterator it1(adj, oneRing); it1; ++it1) {
+        int twoRingi = it1.row();
+        twoRing.insert(twoRingi);
+      }
+    }
+    // Construct P
+    Eigen::MatrixXd P(twoRing.size(), 3);
+    int j = 0;
+    for (auto idx : twoRing) {
+      P.row(j++) = V.row(idx) - V.row(i);
+    }
+
+    // Compute plane passing through V, eigen decomposition
+    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigenSolver(P.transpose() * P);
+    Eigen::VectorXd u = eigenSolver.eigenvectors().col(2);
+    Eigen::VectorXd v = eigenSolver.eigenvectors().col(1);
+    Eigen::VectorXd w = eigenSolver.eigenvectors().col(0);
+    Eigen::VectorXd B = P * w;
+    Eigen::MatrixXd S(twoRing.size(), 2);
+    S.col(0) = P * u;
+    S.col(1) = P * v;
+
+    // Solve for known coefficients
+    Eigen::MatrixXd ax(P.rows(), 5);
+    ax << S, S.col(0).cwiseProduct(S.col(0)), S.col(0).cwiseProduct(S.col(1)), S.col(1).cwiseProduct(S.col(1));
+    Eigen::MatrixXd X;
+    igl::pinv(ax, X);
+    Eigen::VectorXd A = X * B;
+
+    // Shape operator
+    double e, f, g, E, Fd, G;
+    E = 1 + A(0) * A(0);
+    Fd = 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 = (1 * 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::MatrixXd left(2, 2);
+    left << e, f,
+            f, g;
+    Eigen::MatrixXd right(2, 2);
+    right << E, Fd,
+             Fd, G;
+    Eigen::MatrixXd shapeOp(2, 2);
+    shapeOp = - left * right.inverse();
+
+    // Eigen solver
+    Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> eigenSolverSO(shapeOp);
+    K1(i) = eigenSolverSO.eigenvalues()(0);
+    K2(i) = eigenSolverSO.eigenvalues()(1);
+    D1.row(i) = eigenSolverSO.eigenvectors()(1, 1) * u + eigenSolverSO.eigenvectors()(1, 0) * v;
+    D2.row(i) = eigenSolverSO.eigenvectors()(0, 1) * u + eigenSolverSO.eigenvectors()(0, 0) * v;
+  }
+
+}
\ No newline at end of file