Finished the assignment

src/cotmatrix.cpp

@@ -1,10 +1,75 @@
 #include "cotmatrix.h"
 
+/**
+* @brief Solves with the cosine law to produce the cosine of the angle
+* opposite to the edge c.
+*/
+double cosineLaw(double a, double b, double c) {
+	double top = a * a + b * b - c * c;
+	double bottom = 2 * a * b;
+	return top / bottom;
+}
+
+/**
+* @brief Computes the cotagenent of the angle opposite to the edges with the
+* given lengths in a triangle.
+*
+* @param lengths input array of 3 edge lengths.
+* @param angles output array of 3 angles.
+*/
+void lengthsToCotangents(double* lengths, double* cots) {
+
+	// First we apply cosine law.
+	double cosines[3];
+	cosines[0] = cosineLaw(lengths[1], lengths[2], lengths[0]);
+	cosines[1] = cosineLaw(lengths[2], lengths[0], lengths[1]);
+	cosines[2] = cosineLaw(lengths[0], lengths[1], lengths[2]);
+
+	// Now we apply pythagorean and cotangent identities to produce the
+	// cotangent.
+	for (int i = 0; i < 3; i++) {
+
+		// Pythaogrean identity is valid because the angle is known to be 
+		// less than PI; so sine wouldn't have gone negative in there.
+		double cosine = cosines[i];
+		double sine = sqrt(1 - cosine * cosine);
+		double cotangent = cosine / sine;
+		cots[i] = cotangent;
+	}
+}
+
 void cotmatrix(
-  const Eigen::MatrixXd & l,
-  const Eigen::MatrixXi & F,
-  Eigen::SparseMatrix<double> & L)
-{
-  // Add your code here
+		const Eigen::MatrixXd & l,
+		const Eigen::MatrixXi & F,
+		Eigen::SparseMatrix<double> & L) {
+
+	int vertexCount = F.maxCoeff() + 1;
+	int faceCount = l.rows();
+
+	L.resize(vertexCount, vertexCount);
+	for (int i = 0; i < faceCount; i++) {
+
+		// Find cotangents opposite to each edge.
+		double lengths[] = { l(i, 0), l(i, 1), l(i, 2) };
+		double cots[3];
+		lengthsToCotangents(lengths, cots);
+
+		// For each edge, fill in the appropriate matrix stuff.
+		L.coeffRef(F(i, 0), F(i, 1)) += 0.5 * cots[2];
+		L.coeffRef(F(i, 1), F(i, 2)) += 0.5 * cots[0];
+		L.coeffRef(F(i, 2), F(i, 0)) += 0.5 * cots[1];
+
+		// And its symmetric after all...
+		L.coeffRef(F(i, 1), F(i, 0)) += 0.5 * cots[2];
+		L.coeffRef(F(i, 2), F(i, 1)) += 0.5 * cots[0];
+		L.coeffRef(F(i, 0), F(i, 2)) += 0.5 * cots[1];
+	}
+
+	// Now set the diagonal
+	for (int i = 0; i < vertexCount; i++) {
+		L.insert(i, i) = -L.row(i).sum();
+	}
+
+
 }
 

src/massmatrix.cpp

@@ -1,10 +1,35 @@
 #include "massmatrix.h"
 
 void massmatrix(
-  const Eigen::MatrixXd & l,
-  const Eigen::MatrixXi & F,
-  Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
-{
-  // Add your code here
+		const Eigen::MatrixXd & l,
+		const Eigen::MatrixXi & F,
+		Eigen::DiagonalMatrix<double, Eigen::Dynamic> & M) {
+
+	int vertexCount = F.maxCoeff() + 1;
+	int triangleCount = F.rows();
+
+	M.resize(vertexCount);
+	for (int i = 0; i < vertexCount; i++) {
+
+		double totalArea = 0;
+		for (int t = 0; t < triangleCount; t++) {
+
+			// Check if the vertex i is contained in the triangle t.
+			if (F(t, 0) != i && F(t, 1) != i && F(t, 2) != i) {
+				continue;
+			}
+
+			// Compute area of the triangle from the side lengths.
+			// Use Heron's formula.
+			double s = (l(i, 0) + l(i, 1) + l(i, 2)) / 2;
+			double area = sqrt(s * (s - l(i, 0)) * (s - l(i, 1)) * 
+					(s - l(i, 2)));
+			totalArea += area;
+		}
+		M.diagonal()[i] = totalArea;
+	}
+
+	M = M * (1 / 3.0);
+
 }
 

src/smooth.cpp

@@ -1,12 +1,39 @@
 #include "smooth.h"
+#include "igl/edge_lengths.h"
+#include "cotmatrix.h"
+#include "massmatrix.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;
+		const Eigen::MatrixXd & V,
+		const Eigen::MatrixXi & F,
+		const Eigen::MatrixXd & G,
+		double lambda,
+		Eigen::MatrixXd & U) {
+
+	// Compute matrix of edge lengths.
+	Eigen::MatrixXd edgeLengths;
+	edgeLengths.resizeLike(F);
+	igl::edge_lengths(V, F, edgeLengths);
+
+	// Compute Laplacian and Mass Matrices
+	Eigen::SparseMatrix<double> laplacian;
+	Eigen::DiagonalMatrix<double, Eigen::Dynamic> mass;
+	cotmatrix(edgeLengths, F, laplacian);
+	massmatrix(edgeLengths, F, mass);
+
+	// Treat the system given as a linear system we can solve with Ax = b.
+	// I construct A weirdly because of the data-types involved.
+	Eigen::MatrixXd b = mass * G;
+	Eigen::SparseMatrix<double> A = -lambda * laplacian;
+	for (int i = 0; i < mass.rows(); i++) {
+		A.coeffRef(i, i) += mass.diagonal()[i];
+	}
+
+	// Now solve A to obtain G for the next iteration.
+	Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> solver;
+	solver.compute(A);
+
+	U.resizeLike(G);
+	U = solver.solve(b);
+	return;
 }