Sarah Kushner assignment submission

src/cotmatrix.cpp

@@ -5,6 +5,50 @@ void cotmatrix(
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double> & L)
 {
+
   // Add your code here
+	L.resize(F.maxCoeff() + 1, F.maxCoeff() + 1);
+
+	float e_i, e_j, e_ij = 0;
+	float cos_alpha, sin_alpha, cot_alpha = 0;
+
+	typedef Eigen::Triplet<double> T;
+	std::vector<T> tripletList;
+	tripletList.reserve((F.maxCoeff()+1)*12);
+
+	auto get_cot = [&](int face, int a, int b, int c)
+	{
+		e_i = l(face, a);
+		e_j = l(face, b);
+		e_ij = l(face, c);
+		cos_alpha = (pow(e_i, 2) + pow(e_j, 2) - pow(e_ij, 2)) / (2 * e_i*e_j);
+		sin_alpha = sqrt(1 - pow(cos_alpha, 2));
+		cot_alpha = cos_alpha / sin_alpha;
+		return cot_alpha / 2.;
+	};
+
+	for (int i = 0; i < F.rows(); i++)
+	{
+		auto cot = get_cot(i, 0, 1, 2);
+		tripletList.push_back(T(F(i, 0), F(i, 1), cot));
+		tripletList.push_back(T(F(i, 1), F(i, 0), cot));
+		tripletList.push_back(T(F(i, 0), F(i, 0), -cot));
+		tripletList.push_back(T(F(i, 1), F(i, 1), -cot));
+
+		cot = get_cot(i, 1, 2, 0);
+		tripletList.push_back(T(F(i, 1), F(i, 2), cot));
+		tripletList.push_back(T(F(i, 2), F(i, 1), cot));
+		tripletList.push_back(T(F(i, 1), F(i, 1), -cot));
+		tripletList.push_back(T(F(i, 2), F(i, 2), -cot));
+
+		cot = get_cot(i, 2, 0, 1);
+		tripletList.push_back(T(F(i, 2), F(i, 0), cot));
+		tripletList.push_back(T(F(i, 0), F(i, 2), cot));
+		tripletList.push_back(T(F(i, 0), F(i, 0), -cot));
+		tripletList.push_back(T(F(i, 2), F(i, 2), -cot));
+	}
+
+
+	L.setFromTriplets(tripletList.begin(), tripletList.end());
 }
 

src/massmatrix.cpp

@@ -1,10 +1,45 @@
 #include "massmatrix.h"
+#include <igl/sort.h>
 
 void massmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
   // Add your code here
+  M.resize(F.maxCoeff() + 1);
+  M.diagonal() = Eigen::VectorXd::Zero(F.maxCoeff() + 1);
+
+  float a, b, c, area = 0;
+  Eigen::ArrayXi face(3);
+  Eigen::VectorXf unsorted(3);
+  Eigen::VectorXf sorted(3);
+  Eigen::VectorXi index;
+
+  int v1, v2, v3 = 0;
+  int s;
+
+  for (int i = 0; i < F.rows(); i++)
+  {
+	  v1 = F(i, 0);
+	  v2 = F(i, 1);
+	  v3 = F(i, 2);
+
+	  unsorted << l(i, 0), l(i, 1), l(i, 2);
+	  igl::sort(unsorted, 1, 0, sorted, index);
+
+	  a = sorted[0];
+	  b = sorted[1];
+	  c = sorted[2];
+
+	  s = (a + b + c) / 2;
+	  //area = sqrt(s*(s-a)*(s-b)*(s-c));
+	  area = sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c))) / 4 ; // stable Heron's formula
+	  M.diagonal()(v1) += area;
+	  M.diagonal()(v2) += area;
+	  M.diagonal()(v3) += area;
+  }
+
+  M.diagonal() /= 3;
 }
 

src/smooth.cpp

@@ -1,4 +1,8 @@
 #include "smooth.h"
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include <igl/edge_lengths.h>
+#include <Eigen\src\Cholesky\LLT.h>
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -8,5 +12,30 @@ void smooth(
     Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = G;
+
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V, F, l);
+
+  Eigen::SparseMatrix<double> L;
+  cotmatrix(l, F, L);
+
+  Eigen::DiagonalMatrix<double, Eigen::Dynamic> M;
+  massmatrix(l, F, M);
+
+  Eigen::SparseMatrix<double> A;
+  A = -lambda * L;
+
+
+  for (int i = 0; i < A.rows(); i++)
+  {
+	  A.coeffRef(i, i) += M.diagonal()(i);
+  }
+  //Eigen::MatrixXd A = M - lambda * L;
+  Eigen::LLT<Eigen::MatrixXd> llt;
+  llt.compute(A);
+
+  Eigen::MatrixXd b = M*G;
+  U = llt.solve(b);
+
+
 }