Jiayi(Eris) Zhang HW6

src/angle_defect.cpp

@@ -1,9 +1,29 @@
 #include "../include/angle_defect.h"
+#include "internal_angles.h"
+#include <igl/squared_edge_lengths.h>
+#include <math.h>
+
+using namespace Eigen;
 
 void angle_defect(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::VectorXd & D)
 {
-  D = Eigen::VectorXd::Zero(V.rows());
+
+  // get pi
+  double pi = acos(-1);
+
+  D = VectorXd::Ones(V.rows()) * pi * 2;
+
+  Eigen::MatrixXd l_sqr, A;
+  igl::squared_edge_lengths(V, F, l_sqr);
+  internal_angles(l_sqr, A); // get internal angles
+
+  for (int i = 0; i < F.rows(); i++) {
+    for (int j= 0; j < 3; j++) {
+      D(F(i,j)) = D(F(i,j)) - A(i, j);
+    }
+  }
+
 }

src/internal_angles.cpp

@@ -1,8 +1,23 @@
 #include "../include/internal_angles.h"
+#include <math.h> 
 
 void internal_angles(
   const Eigen::MatrixXd & l_sqr,
   Eigen::MatrixXd & A)
 {
-  // Add with your code
+  int num = l_sqr.rows();
+  A.resize(num, 3);
+
+  for (int i = 0; i < num; i++) {
+    for (int j = 0; j < 3; j++) {
+      // side: b, c   angle to: a
+      double a = l_sqr(i, j % 3);
+      double b = l_sqr(i, (j + 1) % 3);
+      double c = l_sqr(i, (j + 2) % 3);
+
+      double cosAngle = ((b + c - a) * 1.0 / (2 * sqrt(b * c)));
+      A(i, j) = acos(cosAngle);
+    }
+  }
+
 }

src/mean_curvature.cpp

@@ -1,10 +1,46 @@
 #include "../include/mean_curvature.h"
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/invert_diag.h>
+#include <igl/per_vertex_normals.h>
+
+using namespace Eigen;
 
 void mean_curvature(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::VectorXd & H)
 {
-  // Replace with your code
+
   H = Eigen::VectorXd::Zero(V.rows());
+
+  // cotangent matrix
+  SparseMatrix<double> L;
+	igl::cotmatrix(V, F, L);
+
+  // mass matrix
+  SparseMatrix<double> M, M_inv;
+  igl::massmatrix(V, F, igl::MassMatrixType::MASSMATRIX_TYPE_DEFAULT, M);
+  igl::invert_diag(M, M_inv);
+
+  // normals: to determine the sign
+  MatrixXd N;
+  igl::per_vertex_normals(V, F, N);
+
+  // curvature normals
+  MatrixXd HN = M_inv * L * V;
+
+  // check and update H
+  for (int i = 0; i < V.rows(); i++) {
+    double cur = HN.row(i).dot(N.row(i));
+    int norm = HN.row(i).norm();
+    if (cur > 0) {
+      H(i) = -norm;
+    } 
+    else {
+      H(i) = norm;
+    }
+  }
+
+
 }

src/principal_curvatures.cpp

@@ -1,4 +1,12 @@
 #include "../include/principal_curvatures.h"
+#include <igl/adjacency_matrix.h>
+#include <igl/per_vertex_normals.h>
+#include <set>
+#include <igl/pinv.h>
+
+
+using namespace Eigen;
+using namespace std;
 
 void principal_curvatures(
   const Eigen::MatrixXd & V,
@@ -13,4 +21,84 @@ void principal_curvatures(
   K2 = Eigen::VectorXd::Zero(V.rows());
   D1 = Eigen::MatrixXd::Zero(V.rows(),3);
   D2 = Eigen::MatrixXd::Zero(V.rows(),3);
+
+  // compute the adjacency matrix: prepare for the two rings
+  SparseMatrix<int> A;
+  igl::adjacency_matrix(F, A);
+
+  // compute the per vertex normals for the mesh
+  MatrixXd N;
+  igl::per_vertex_normals(V, F, N);
+
+  // iterate over all vertices and solve for the quadratic surface
+  for (int i = 0; i < V.rows(); i++) {
+
+    // construct P matrix
+    MatrixXd P;
+    set<int> rings;
+    for (SparseMatrix<int>::InnerIterator iter1(A, i); iter1; ++iter1) {
+      int cur1 = iter1.row();
+      for (SparseMatrix<int>::InnerIterator iter2(A, cur1); iter2; ++iter2) {
+          int cur2 = iter2.row();
+          rings.insert(cur2);
+      }
+    }
+    // discard i
+    rings.erase(i);
+    P.resize(rings.size(), 3);
+    int j = 0;
+    for (set<int>::iterator iter = rings.begin(); iter != rings.end(); ++iter, j++) {
+      int cur = *iter;
+      P.row(j) = V.row(cur) - V.row(i);
+    }
+
+    // PCA
+    SelfAdjointEigenSolver<Matrix3d> PCA(P.transpose() * P);
+    MatrixXd evs = PCA.eigenvectors();
+		Vector3d u = evs.col(2);
+    Vector3d v = evs.col(1);
+    Vector3d w = evs.col(0);
+
+    // uv terms
+    MatrixXd UV(P.rows(), 5);
+    UV.col(0) = P * u;
+    UV.col(1) = P * v;
+    UV.col(2) = (UV.col(0).array().square()).matrix();
+    UV.col(3) = (UV.col(0).array() * UV.col(1).array()).matrix();
+    UV.col(4) = (UV.col(1).array().square()).matrix();
+
+    // solve for coefficients
+    MatrixXd UV_inv;
+		igl::pinv(UV, UV_inv);
+		VectorXd a = UV_inv * (P * w);
+
+    // construct shape operators
+    double E = 1 + a(0) * a(0);
+    double F = a(0) * a(1);
+    double G = 1 + a(1) * a(1);
+    double denom = sqrt(a(0) * a(0) + 1 + a(1) * a(1));
+    double e = 2.0 * a(2) / denom;
+    double f = a(3) * 1.0 / denom;
+    double g = 2.0 * a(4) / denom;
+
+    MatrixXd MatL(2, 2), MatU(2, 2);
+    MatL << e, f, f, g;
+    MatU << E, F, F, G;
+    MatrixXd S = -MatL * MatU.inverse();
+
+    SelfAdjointEigenSolver<MatrixXd> ei(S);
+
+    // update
+    double k1 = ei.eigenvalues()(1);
+    double k2 = ei.eigenvalues()(0);
+    K1(i) = k1;
+    K2(i) = k2;
+
+    Vector3d d1 = ei.eigenvectors()(1, 1) * u + ei.eigenvectors()(1, 0) * v;
+    Vector3d d2 = ei.eigenvectors()(0, 1) * u + ei.eigenvectors()(0, 0) * v;
+    D1.row(i) = d1;
+    D2.row(i) = d2;
+
+  }
+
 }