wenzheng chen hw6

src/angle_defect.cpp

@@ -1,9 +1,31 @@
 #include "../include/angle_defect.h"
 
+#include <igl/squared_edge_lengths.h>
+#include "internal_angles.h"
+#include <Eigen/Dense>
+
 void angle_defect(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::VectorXd & D)
 {
   D = Eigen::VectorXd::Zero(V.rows());
+
+  // calculate edge length
+  Eigen::MatrixXd Flen;
+  igl::squared_edge_lengths(V, F, Flen);
+
+  Eigen::MatrixXd VAngles;
+  internal_angles(Flen, VAngles);
+
+  D = Eigen::VectorXd::Ones(V.rows()) * 3.1416;
+  for (int i = 0; i< F.rows(); i++)
+  {
+	  for(int j = 0; j< 3; j++)
+	  {
+		  D(F(i, j)) -= VAngles(F(i, j));
+	  }
+  }
+
+
 }

src/internal_angles.cpp

@@ -1,8 +1,29 @@
 #include "../include/internal_angles.h"
 
-void internal_angles(
-  const Eigen::MatrixXd & l_sqr,
-  Eigen::MatrixXd & A)
-{
-  // Add with your code
+#include <Eigen/Dense>
+
+void internal_angles(const Eigen::MatrixXd & l_sqr, Eigen::MatrixXd & A) {
+	// Add with your code
+
+	int fnum = l_sqr.rows();
+	A.resize(fnum, 3);
+	for (int i = 0; i < fnum; i++) {
+		// c^2 = a^2 + b^2 - 2abcos
+		// length: 12, 23, 31
+		// angle order, c, a, b
+		double f12 = l_sqr(i, 0);
+		double f23 = l_sqr(i, 1);
+		double f31 = l_sqr(i, 2);
+
+		double a = (f12 * f12 + f31 * f31 - f23 * f23) / 2 / f12 / f31;
+		A(i, 0) = std::acos(a);
+
+		double b = (f12 * f12 + f23 * f23 - f31 * f31) / 2 / f12 / f23;
+		A(i, 1) = std::acos(b);
+
+		double c = (f23 * f23 + f31 * f31 - f12 * f12) / 2 / f23 / f31;
+		A(i, 2) = std::acos(c);
+	}
+
+	return;
 }

src/mean_curvature.cpp

@@ -1,10 +1,44 @@
 #include "../include/mean_curvature.h"
 
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <Eigen/Dense>
+#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());
+
+// calculate M and L
+  	// calculate L, here, w is 1
+  	Eigen::SparseMatrix<double> L;
+  	igl::cotmatrix(V, F, L);
+
+  	// M
+  	Eigen::SparseMatrix<double> M;
+  	igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_BARYCENTRIC, M);
+
+  	// MH = LV
+  	Eigen::MatrixXd b = L * V;
+
+  	// svd
+  	H = M.colPivHouseholderQr().solve(b);
+
+  	// sign
+  	Eigen::MatrixXd vnorms;
+  	igl::per_vertex_normals(V, F, igl::PER_VERTEX_NORMALS_WEIGHTING_TYPE_AREA, vnorms);
+  	int vnum = V.rows();
+  	for (int i = 0; i< vnum; i++)
+  	{
+  		double sum = H.row(i).transpose() * vnorms.row(i);
+  		if (sum > 0)
+  		{
+  			H.row(i) = -H.row(i);
+  		}
+  	}
+
 }