Adrian She Submission

src/angle_defect.cpp

@@ -1,9 +1,23 @@
 #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());
+  // Get list of internal angles
+  Eigen::MatrixXd L;
+  igl::squared_edge_lengths(V, F, L);
+  Eigen::MatrixXd A;
+  internal_angles(L, A);
+
+  // From angles compute angle defects
+  D = Eigen::VectorXd::Ones(V.rows()) * 2 * std::acos(-1); // acos(-1) = pi
+  for (int i = 0; i < F.rows(); i++ ){
+	for (int j = 0; j <= 2; j++) {
+ 	D(F(i,j)) = D(F(i,j)) - A(i,j);
+	} 
+  }
 }

src/internal_angles.cpp

@@ -4,5 +4,15 @@ void internal_angles(
   const Eigen::MatrixXd & l_sqr,
   Eigen::MatrixXd & A)
 {
-  // Add with your code
+  A.resize(l_sqr.rows(), 3);
+
+  // Apply the cosine law to compute each internal angle
+  for (int i = 0; i < l_sqr.rows(); i++) {
+	A(i, 0) = (l_sqr(i, 0) - l_sqr(i,1) - l_sqr(i,2)) / (-2.0 * std::sqrt(l_sqr(i,1) * l_sqr(i,2))); 
+	A(i, 1) = (l_sqr(i, 1) - l_sqr(i,2) - l_sqr(i,0)) / (-2.0 * std::sqrt(l_sqr(i,0) * l_sqr(i,2))); 
+	A(i, 2) = (l_sqr(i, 2) - l_sqr(i,0) - l_sqr(i,1)) / (-2.0 * std::sqrt(l_sqr(i,0) * l_sqr(i,1)));
+	A(i, 0) = std::acos(A(i,0));
+ 	A(i, 1) = std::acos(A(i,1));
+	A(i, 2) = std::acos(A(i,2));
+  }
 }

src/mean_curvature.cpp

@@ -1,10 +1,44 @@
 #include "../include/mean_curvature.h"
+#include "igl/cotmatrix.h"
+#include "igl/invert_diag.h"
+#include "igl/massmatrix.h"
+#include "igl/squared_edge_lengths.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
+  // Construct Laplacian and mass matrices
+  Eigen::SparseMatrix<double> L;
+  Eigen::SparseMatrix<double> M;
+  Eigen::SparseMatrix<double> M_inv;
+
+  igl::massmatrix(V, F, igl::MassMatrixType::MASSMATRIX_TYPE_BARYCENTRIC, M);
+  igl::cotmatrix(V, F, L);
+  igl::invert_diag(M, M_inv);
+
+  // Compute mean curvature vectors as 1/2 (M^(-1) L V)
+  Eigen::MatrixXd H_vec(V.rows(), 3);
+  for (int i = 0; i <= 2; i++) {
+	H_vec.col(i) = 0.5 * M_inv * L * V.col(i); 
+  }
+
+  // Produce mean curvature by comparing normals
+  Eigen::MatrixXd N;
+  igl::per_vertex_normals(V, F, N);
+
   H = Eigen::VectorXd::Zero(V.rows());
+  for (int i = 0; i < V.rows(); i++) {
+	if (N(i, 2) > 0) {
+		H(i) = H_vec.row(i).norm();
+	}
+	else if (N(i, 2) < 0) {
+		H(i) = -H_vec.row(i).norm();
+	}
+	else {
+		H(i) = 0;
+	}
+  }
 }

src/principal_curvatures.cpp

@@ -1,4 +1,7 @@
 #include "../include/principal_curvatures.h"
+#include "igl/adjacency_list.h"
+#include "igl/pinv.h"
+#include <Eigen/Eigenvalues>
 
 void principal_curvatures(
   const Eigen::MatrixXd & V,
@@ -8,9 +11,68 @@ void principal_curvatures(
   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);
+
+  std::vector<std::vector<int>> L;
+  igl::adjacency_list(F, L); 
+
+  for (int i = 0; i < V.rows(); i++) {
+	// BFS on vertex v to get the two ring;
+	std::vector<int> V_tr;
+        for (int j = 0; j < L[i].size(); j++) {
+			V_tr.push_back(L[i][j]);
+			for (int k = 0; k < L[j].size(); k++) {
+					V_tr.push_back(L[j][k]);
+			}
+		}
+
+
+ 	// Compute matrix P of offset positions 
+	Eigen::MatrixXd P(V_tr.size(), 3);
+	for (int j = 0; j < V_tr.size(); j++) {
+		P.row(j) = V.row(V_tr[j]) - V.row(i);
+	}	
+
+	// Compute eigenvectors of P^T * P to get principal directions u, v , w
+        Eigen::MatrixXd D = P.transpose() * P;
+	Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(D);
+	Eigen::Vector3d u = eigensolver.eigenvectors().col(2);
+	Eigen::Vector3d v = eigensolver.eigenvectors().col(1);
+	Eigen::Vector3d w = eigensolver.eigenvectors().col(0);
+  	Eigen::MatrixXd A(V_tr.size(), 5);
+
+	// Compute the coordinates of points in P in (u, v, w) basis
+	Eigen::VectorXd b(V_tr.size());
+	for (int i = 0; i < V_tr.size(); i++) {
+		A(i, 0) = P.row(i).dot(u);
+		A(i, 1) = P.row(i).dot(v);
+		b(i) = P.row(i).dot(w);
+	}
+
+	// Construct quadratic height field coefficients
+	A.col(2) = (A.col(0)).cwiseProduct(A.col(0));
+	A.col(3) = (A.col(1)).cwiseProduct(A.col(0));
+	A.col(4) = (A.col(1)).cwiseProduct(A.col(1));
+	Eigen::MatrixXd A_inv;
+	igl::pinv(A, -1, A_inv); // use default tolerance
+	Eigen::VectorXd a = A_inv *  b;
+
+	// compute shape operator
+	Eigen::Matrix2d first_ff;
+	first_ff << 1.0 + a(0) * a(0), a(0) * a(1), a(0) * a(1), 1.0 + a(1) * a(1); 
+	double div = std::sqrt(1 + a(0) * a(0) + a(1) * a(1));
+	Eigen::Matrix2d second_ff;
+	second_ff << 2.0 * a(2) / div, a(3) / div, a(3) / div, 2.0 * a(4) /div;
+	Eigen::Matrix2d S = - second_ff * first_ff.inverse();
+
+	// eigen decompose the shape operator to obtain curvatures and directions
+	eigensolver.compute(S);
+	K1(i) = eigensolver.eigenvalues()[1];
+	K2(i) = eigensolver.eigenvalues()[0];
+	D1.row(i) = eigensolver.eigenvectors()(1,0) * u + eigensolver.eigenvectors()(1,1) * v ;
+	D2.row(i) = eigensolver.eigenvectors()(0,0) * u + eigensolver.eigenvectors()(0,1) * v; 
+ } 
 }