Sajjad Nezhadi HW6

src/angle_defect.cpp

@@ -1,9 +1,27 @@
 #include "../include/angle_defect.h"
+#include "../include/internal_angles.h"
+#include <igl/squared_edge_lengths.h>
+#include <iostream>
 
 void angle_defect(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::VectorXd & D)
 {
   D = Eigen::VectorXd::Zero(V.rows());
+  Eigen::MatrixXd K, L;
+
+  igl::squared_edge_lengths(V, F, L);
+  internal_angles(L, K);
+
+   for (int i = 0; i< V.rows(); i++){
+    D(i) = 2*M_PI;
+   }
+
+   for (int i = 0; i < F.rows(); i++) {
+    for(int j= 0; j < 3; j++){
+        D(F(i,j)) -= K(i, j);
+      }
+    }
+    std::cout << D << std::endl;
 }

src/internal_angles.cpp

@@ -1,8 +1,26 @@
 #include "../include/internal_angles.h"
+#include <igl/cotmatrix.h>
+#include <iostream>
 
 void internal_angles(
   const Eigen::MatrixXd & l_sqr,
   Eigen::MatrixXd & A)
 {
-  // Add with your code
+    A.resize(l_sqr.rows(), 3);
+
+    for (int t = 0; t < l_sqr.rows(); t++){
+
+        double a = l_sqr(t, 0);
+        double b = l_sqr(t, 1);
+        double c = l_sqr(t, 2);
+
+        double cos_aij = ((a*a) + (b*b) - (c*c))/(2*a*b);
+        double cos_aik = ((c*c) + (b*b) - (a*a))/(2*c*b);
+        double cos_ajk = ((c*c) + (a*a) - (b*b))/(2*c*a);
+
+        A(t, 0) = acos(cos_aij);
+        A(t, 1) = acos(cos_aik);
+        A(t, 2) = acos(cos_ajk);
+    }
+
 }

src/mean_curvature.cpp

@@ -1,4 +1,9 @@
 #include "../include/mean_curvature.h"
+#include <igl/invert_diag.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <iostream>
 
 void mean_curvature(
   const Eigen::MatrixXd & V,
@@ -7,4 +12,20 @@ void mean_curvature(
 {
   // Replace with your code
   H = Eigen::VectorXd::Zero(V.rows());
+  Eigen::SparseMatrix<double> L;
+  Eigen::SparseMatrix<double> M;
+  Eigen::SparseMatrix<double> Y;
+  Eigen::MatrixXd N;
+
+  igl::cotmatrix(V, F, L); 
+  igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+  igl::invert_diag(M, Y);
+  igl::per_vertex_normals(V, F, N);
+  std::cout << Y.cols() << std::endl;
+  std::cout << L.rows() << std::endl;
+  std::cout << V.rows() << std::endl;
+
+  Eigen::MatrixXd temp = Y * L * V * N.transpose();
+  for (int i = 0; i < V.rows(); i++)
+    H(i) = -temp(i,i);
 }

src/principal_curvatures.cpp

@@ -1,4 +1,7 @@
 #include "../include/principal_curvatures.h"
+#include <igl/adjacency_matrix.h>
+#include <igl/pinv.h>
+#include <iostream>
 
 void principal_curvatures(
   const Eigen::MatrixXd & V,
@@ -13,4 +16,56 @@ 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> A;
+  igl::adjacency_matrix(F, A);
+  A = A * A;
+
+  for(int i=0; i<V.rows(); i++){
+    Eigen::SparseVector<int> Ring = A.col(i);
+    Eigen::MatrixXd P(Ring.nonZeros(), 3);
+
+    int k = 0;
+    for (Eigen::SparseVector<double>::InnerIterator it(Ring); it; ++it)
+        P.row(k++) = V.row(it.index()) - V.row(i);
+
+    Eigen::JacobiSVD<Eigen::MatrixXd> svd(P.transpose() * P, Eigen::ComputeFullU);
+    Eigen::MatrixXd SD = svd.matrixU();
+    std::cout << SD.cols() << std::endl;
+
+    Eigen::VectorXd t_1 = P * SD.col(0);
+    Eigen::VectorXd t_2 = P * SD.col(1);
+    Eigen::VectorXd t_3 = t_1.array().colwise()*t_1.array();
+    Eigen::VectorXd t_4 = t_1.array().colwise()*t_2.array();
+    Eigen::VectorXd t_5 = t_2.array().colwise()*t_2.array();
+    Eigen::MatrixXd t(P.rows(), 5);
+    t<< t_1, t_2, t_3, t_4, t_5;
+
+    Eigen::MatrixXd a;
+    igl::pinv(t, a);
+    a = a * P * SD.col(2);
+
+    double tsqrt = sqrt(a(0) * a(0) + 1 + a(1) * a(1));
+    double E = 1 + a(0) * a(0);
+    double F1 = a(0) * a(1);
+    double G = 1 + a(1) * a(1);
+    double e = 2 * a(2) / tsqrt ;
+    double f = a(3) / tsqrt;
+    double g = 2 * a(4) / tsqrt;
+
+    Eigen::MatrixXd left(2, 2);
+    left << e, f,
+          f, g;
+    Eigen::MatrixXd right(2, 2);
+    right << E, F1, 
+          F1, G;
+
+    Eigen::MatrixXd S = -1*left*right.inverse();
+    Eigen::ComplexEigenSolver<Eigen::MatrixXd> ces;
+    ces.compute(S);
+
+    K1(i) = ces.eigenvalues()[0].real();
+    K2(i) = ces.eigenvalues()[1].real();
+    D1.row(i) = SD.leftCols(2) * ces.eigenvectors().col(0).real();
+    D2.row(i) = SD.leftCols(2) * ces.eigenvectors().col(1).real();
+  }
 }