wenzheng chen hw3

main.cpp

@@ -2,7 +2,7 @@
 #include <igl/read_triangle_mesh.h>
 #include <igl/readDMAT.h>
 #include <igl/parula.h>
-#include <igl/viewer/Viewer.h>
+#include <igl/opengl/glfw/Viewer.h>
 #include <Eigen/Core>
 #include <string>
 #include <iostream>
@@ -31,7 +31,7 @@ int main(int argc, char *argv[])
     }
   }
 
-  igl::viewer::Viewer viewer;
+  igl::opengl::glfw::Viewer viewer;
   std::cout<<R"(
   D,d  smooth data
   K    decamate(?) lambda
@@ -47,19 +47,19 @@ int main(int argc, char *argv[])
     {
       std::cout<<"Too degenerate to keep smoothing. Better reset"<<std::endl;
     }
-    viewer.data.set_mesh(V,F);
-    viewer.data.compute_normals();
+    viewer.data().set_mesh(V,F);
+    viewer.data().compute_normals();
     Eigen::MatrixXd C;
     igl::parula(U,G.minCoeff(),G.maxCoeff(),C);
-    viewer.data.set_colors(C);
+    viewer.data().set_colors(C);
   };
   const auto & reset = [&]()
   {
     V = OV;
     U = G;
   };
   viewer.callback_key_pressed = 
-    [&](igl::viewer::Viewer &, unsigned int key, int)
+    [&](igl::opengl::glfw::Viewer &, unsigned int key, int)
   {
     switch(key)
     {

src/cotmatrix.cpp

@@ -1,10 +1,104 @@
 #include "cotmatrix.h"
+typedef Eigen::Triplet<double> T;
 
-void cotmatrix(
-  const Eigen::MatrixXd & l,
-  const Eigen::MatrixXi & F,
-  Eigen::SparseMatrix<double> & L)
-{
-  // Add your code here
+void setV2F(Eigen::MatrixXi & V2F1, Eigen::MatrixXi & V2F2, int i, int j,
+		int k) {
+
+	if (V2F1(i, j) == -1) {
+		V2F1(i, j) = k;
+		V2F1(j, i) = k;
+	} else {
+		V2F2(i, j) = k;
+		V2F2(j, i) = k;
+	}
+	return;
+}
+
+void cotmatrix(const Eigen::MatrixXd & l, const Eigen::MatrixXi & F,
+		Eigen::SparseMatrix<double> & L) {
+	// Add your code here
+
+	int vnum = F.maxCoeff() + 1;
+	int fnum = F.rows();
+
+	// calculate area of faces
+	Eigen::VectorXd FA(fnum);
+	for (int i = 0; i < fnum; i++) {
+		int v1 = F(i, 0);
+		int v2 = F(i, 1);
+		int v3 = F(i, 2);
+		double e1len = l(v1, v2);
+		double e2len = l(v1, v3);
+		double e3len = l(v2, v3);
+		double s = e1len + e2len + e3len;
+		s = s / 2;
+		double A2 = s * (s - e1len) * (s - e2len) * (s - e3len);
+		double A = sqrt(A2);
+		FA[i] = A;
+	}
+
+	// build point2face connection
+	// 2 point i, j will link to 2 faces
+	Eigen::MatrixXi V2F1, V2F2;
+	V2F1 = -Eigen::MatrixXi::Ones(vnum, vnum);
+	V2F2 = -Eigen::MatrixXi::Ones(vnum, vnum);
+	for (int i = 0; i < fnum; i++) {
+		int v1 = F(i, 0);
+		int v2 = F(i, 1);
+		int v3 = F(i, 2);
+		setV2F(V2F1, V2F2, v1, v2, i);
+	}
+
+	std::vector<T> coef;
+	for (int i = 0; i < vnum; i++) {
+		// diag
+		double ii = 0;
+		for (int j = 0; j < vnum; j++) {
+
+			if (i == j)
+				continue;
+
+			// is i and j connected?
+			int k1 = V2F1(i, j);
+			int k2 = V2F2(i, j);
+			// cot c = (a^2 + b^2 - c^2) / 4A
+
+			int k = k1;
+			if (k != -1) {
+				int v1 = F(k, 0);
+				int v2 = F(k, 1);
+				int v3 = F(k, 2);
+				double e1len = l(v1, v2);
+				double e2len = l(v1, v3);
+				double e3len = l(v2, v3);
+				double A = FA[k];
+				double cot = e1len * e1len + e2len * e2len - e3len * e3len;
+				cot = cot / (4 * A);
+				T tmp(i, j, cot);
+				coef.push_back(tmp);
+				ii += cot;
+			}
+
+			k = k2;
+			if (k != -1) {
+				int v1 = F(k, 0);
+				int v2 = F(k, 1);
+				int v3 = F(k, 2);
+				double e1len = l(v1, v2);
+				double e2len = l(v1, v3);
+				double e3len = l(v2, v3);
+				double A = FA[k];
+				double cot = e1len * e1len + e2len * e2len - e3len * e3len;
+				cot = cot / (4 * A);
+				T tmp(i, j, cot);
+				coef.push_back(tmp);
+				ii += cot;
+			}
+		}
+		T tmp(i, i, ii);
+		coef.push_back(tmp);
+	}
+
+	L.setFromTriplets(coef.begin(), coef.end());
 }
 

src/massmatrix.cpp

@@ -1,10 +1,48 @@
 #include "massmatrix.h"
 
-void massmatrix(
-  const Eigen::MatrixXd & l,
-  const Eigen::MatrixXi & F,
-  Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
-{
-  // Add your code here
+void massmatrix(const Eigen::MatrixXd & l, const Eigen::MatrixXi & F,
+		Eigen::DiagonalMatrix<double, Eigen::Dynamic> & M) {
+	// Add your code here
+	int vnum = F.maxCoeff() + 1;
+	int fnum = F.rows();
+
+	// calculate area of faces
+	Eigen::VectorXd FA(fnum);
+	for (int i = 0; i < fnum; i++) {
+		int v1 = F(i, 0);
+		int v2 = F(i, 1);
+		int v3 = F(i, 2);
+		double e1len = l(v1, v2);
+		double e2len = l(v1, v3);
+		double e3len = l(v2, v3);
+		double s = e1len + e2len + e3len;
+		s = s / 2;
+		double A2 = s * (s - e1len) * (s - e2len) * (s - e3len);
+		double A = sqrt(A2);
+		FA[i] = A;
+	}
+
+	Eigen::MatrixXi V2F = Eigen::MatrixXi::Zero(vnum, fnum);
+	for (int i = 0; i < fnum; i++) {
+		int v1 = F(i, 0);
+		int v2 = F(i, 1);
+		int v3 = F(i, 2);
+		V2F(v1, i) = 1;
+		V2F(v2, i) = 1;
+		V2F(v3, i) = 1;
+	}
+	Eigen::VectorXd diag(vnum);
+	for (int i = 0; i < vnum; i++) {
+		double ii = 0;
+		for (int j = 0; j < fnum; j++) {
+			if (V2F(i, j) == 1) {
+				ii += FA[j];
+			}
+		}
+
+		diag[i] = ii / 3;
+	}
+
+	M = diag.asDiagonal();
 }
 

src/smooth.cpp

@@ -1,12 +1,35 @@
 #include "smooth.h"
+#include "massmatrix.h"
+#include "cotmatrix.h"
 
-void smooth(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
-    const Eigen::MatrixXd & G,
-    double lambda,
-    Eigen::MatrixXd & U)
-{
-  // Replace with your code
-  U = G;
+#include "igl/edge_lengths.h"
+#include <Eigen/IterativeLinearSolvers>
+
+void smooth(const Eigen::MatrixXd & V, const Eigen::MatrixXi & F,
+		const Eigen::MatrixXd & G, double lambda, Eigen::MatrixXd & U) {
+	// Replace with your code
+// first, calculate face edges;
+	Eigen::MatrixXd l;
+	igl::edge_lengths(V, F, l);
+
+	// second, calculate L
+	Eigen::SparseMatrix<double> L;
+	cotmatrix(l, F, L);
+
+	// third , M
+	Eigen::DiagonalMatrix<double, Eigen::Dynamic> M;
+	massmatrix(l, F, M);
+
+	// solve equation
+	// next, we combine G with M
+	Eigen::SparseMatrix<double> A = -lambda * L;
+	for (int i = 0; i < V.rows(); i++) {
+		A.coeffRef(i, i) += M.diagonal()(i);
+	}
+	// last, we solve Ax = B
+	Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > solver;
+	solver.compute(A);
+
+	Eigen::MatrixXd B = M * G;
+	U = solver.solve(B);
 }