Randy Hickey HW06

src/angle_defect.cpp

@@ -1,9 +1,21 @@
 #include "../include/angle_defect.h"
+#include "igl/squared_edge_lengths.h"
+#include "internal_angles.h"
+#include "igl/PI.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, X;
+    igl::squared_edge_lengths(V, F, l_sqr);
+    internal_angles(l_sqr, X);
+
+    D = 2 * igl::PI * Eigen::VectorXd::Ones(V.rows());
+    for (int i = 0; i < F.rows(); i++) {
+        D(F(i, 0)) += -X(i, 0);
+        D(F(i, 1)) += -X(i, 1);
+        D(F(i, 2)) += -X(i, 2);
+    }
 }

src/internal_angles.cpp

@@ -4,5 +4,10 @@ void internal_angles(
   const Eigen::MatrixXd & l_sqr,
   Eigen::MatrixXd & A)
 {
-  // Add with your code
+    A.resize(l_sqr.rows(), l_sqr.cols());
+    for (int i = 0; i < l_sqr.rows(); i++) {
+        A(i, 0) = acos((l_sqr(i, 1) + l_sqr(i, 2) - l_sqr(i, 0)) / (2.0 * sqrt(l_sqr(i, 1) * l_sqr(i, 2))));
+        A(i, 1) = acos((l_sqr(i, 2) + l_sqr(i, 0) - l_sqr(i, 1)) / (2.0 * sqrt(l_sqr(i, 2) * l_sqr(i, 0))));
+        A(i, 2) = acos((l_sqr(i, 0) + l_sqr(i, 1) - l_sqr(i, 2)) / (2.0 * sqrt(l_sqr(i, 0) * l_sqr(i, 1))));
+    }
 }

src/mean_curvature.cpp

@@ -1,10 +1,30 @@
 #include "../include/mean_curvature.h"
+#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, M, M_i;
+    Eigen::MatrixXd curve_normals, normals;
+    igl::cotmatrix(V, F, L);
+    igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_VORONOI, M);
+    igl::invert_diag(M, M_i);
+    curve_normals = (M_i * L) * V;
+    igl::per_vertex_normals(V, F, normals);
+
+    H.resize(curve_normals.rows());
+    for (int i = 0; i < V.rows(); i++) {
+        H(i) = 0.0;
+        if ((normals.row(i)).dot(curve_normals.row(i)) > 0) {
+            H(i) = (curve_normals.row(i)).norm();
+        }
+        else if ((normals.row(i)).dot(curve_normals.row(i)) < 0) {
+            H(i) = - ((curve_normals.row(i)).norm());
+        }
+    }
 }

src/principal_curvatures.cpp

@@ -1,4 +1,8 @@
 #include "../include/principal_curvatures.h"
+#include <set>
+#include "igl/adjacency_matrix.h"
+#include "igl/per_vertex_normals.h"
+#include "igl/pinv.h"
 
 void principal_curvatures(
   const Eigen::MatrixXd & V,
@@ -8,9 +12,62 @@ 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);
+    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);
+    Eigen::SparseMatrix<double> Ad;
+    Eigen::MatrixXd normals;
+    igl::adjacency_matrix(F, Ad);
+    igl::per_vertex_normals(V, F, normals);
+
+    for (int i = 0; i < V.rows(); i++) {
+        Eigen::MatrixXd P, X, X_i;
+        Eigen::VectorXd u, v, w, A;
+        Eigen::Matrix2d S, S_left, S_right;
+        std::set<int> neighbour, neighbour2;
+        for (Eigen::SparseMatrix<double>::InnerIterator it(Ad,i); it; ++it) {
+            if (it.row() != i) {
+                neighbour.insert(it.row());
+                neighbour2.insert(it.row());
+            }
+        }
+        for (int n: neighbour) {
+            for (Eigen::SparseMatrix<double>::InnerIterator it(Ad,n); it; ++it) {
+                if (it.row() != n) {
+                    neighbour2.insert(it.row());
+                }
+            }
+        }
+        P.resize(neighbour2.size(), 3);
+        int z = 0;
+        for (int y: neighbour2) {
+            P.row(z++) = V.row(y) - V.row(i);
+        }
+
+        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(P.transpose() * P);
+        w = es.eigenvectors().col(0);
+        v = es.eigenvectors().col(1);
+        u = es.eigenvectors().col(2);
+        if (w.dot(normals.row(i)) < 0) w *= -1;
+        X.resize(P.rows(), 5);
+        X.col(0) = P * v;
+        X.col(1) = P * u;
+        X.col(2) = (X.col(0)).array().square();
+        X.col(3) = (X.col(0)).array() * (X.col(1)).array();
+        X.col(4) = (X.col(1)).array().square();
+        igl::pinv(X, X_i);
+        A = X_i * (P * w);
+
+        double d = sqrt(A(0) * A(0) + 1 + A(1) * A(1));
+        S_right << (1 + A(0) * A(0)), A(0) * A(1), A(0) * A(1), 1 + A(1) * A(1);
+        S_left << (2 * A(2) / d), (A(3) / d), (A(3) / d), (2 * A(4) / d);
+        S = -S_left * S_right.inverse();
+
+        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es2(S);
+        K1(i) = es2.eigenvalues()(1);
+        D1.row(i) = es.eigenvectors()(1, 1) * u + es.eigenvectors()(1, 0) * v;
+        K2(i) = es2.eigenvalues()(0);
+        D2.row(i) = es.eigenvectors()(0, 1) * u + es.eigenvectors()(0, 0) * v;
+    }
 }