From eee48d7b16c1cbe710156f714876e2c6f9dfe3bb Mon Sep 17 00:00:00 2001
From: Yihuan Mao <yihuan.mao@mail.utoronto.ca>
Date: Sun, 7 Oct 2018 00:05:35 -0400
Subject: [PATCH] Yihuan Mao submission

---
 main.cpp           |  2 +-
 src/cotmatrix.cpp  | 25 ++++++++++++++++++++++++-
 src/massmatrix.cpp | 14 +++++++++++++-
 src/smooth.cpp     | 24 ++++++++++++++++++++++--
 4 files changed, 60 insertions(+), 5 deletions(-)

diff --git a/main.cpp b/main.cpp
index c1d5fed..e9401d5 100644
--- a/main.cpp
+++ b/main.cpp
@@ -14,7 +14,7 @@ int main(int argc, char *argv[])
   Eigen::MatrixXi F;
   double lambda = 1e-5;
   igl::read_triangle_mesh(
-    (argc>1?argv[1]:"../shared/data/bunny.off"),OV,F);
+    (argc>1?argv[1]:"../shared/data/decimated-knight.off"),OV,F);
   // Load data into MatrixXd rather than VectorXd for simpler `smooth` API
   // Just use y-coordinates as data to be smoothed
   Eigen::MatrixXd G = OV.col(1);
diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..ecdbb28 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,10 +1,33 @@
 #include "cotmatrix.h"
+#include "math.h"
+
+double cot(double a,double b,double c){
+	double s=(a+b+c)/2;
+	double area=sqrt(s*(s-a)*(s-b)*(s-c));
+	double sin=area*2/b/c;
+	double cos=(b*b+c*c-a*a)/2/b/c;
+	return cos/sin;
+}
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double> & L)
 {
-  // Add your code here
+  L.resize(F.maxCoeff()+1,F.maxCoeff()+1);
+  int f=F.rows();
+  double c;
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> list;
+  for (int i=0; i<f; i++){
+  	for (int j=0; j<3; j++){
+  	  c=cot(l(i,j),l(i,(j+1)%3),l(i,(j+2)%3))/2;
+  	  list.push_back(T(F(i,(j+1)%3),F(i,(j+2)%3),c));
+  	  list.push_back(T(F(i,(j+2)%3),F(i,(j+1)%3),c));
+  	  list.push_back(T(F(i,(j+1)%3),F(i,(j+1)%3),-c));
+  	  list.push_back(T(F(i,(j+2)%3),F(i,(j+2)%3),-c));
+	}
+  }
+  L.setFromTriplets(list.begin(), list.end());
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..ead81fd 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,10 +1,22 @@
 #include "massmatrix.h"
+#include "math.h"
 
 void massmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
-  // Add your code here
+  int f=F.rows();
+  M.resize(F.maxCoeff()+1);
+  Eigen::VectorXd area=Eigen::VectorXd::Zero(F.maxCoeff()+1);
+  double s,a;
+  for (int i=0; i<f; i++){
+  	s=(l(i,0)+l(i,1)+l(i,2))/2;
+  	a=sqrt(s*(s-l(i,0))*(s-l(i,1))*(s-l(i,2))); 
+  	for (int j=0; j<3; j++){
+  	  area(F(i,j))+=a;
+	}
+  }
+  M=1.0/3*area.asDiagonal();
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..b3bad9a 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,4 +1,8 @@
 #include "smooth.h"
+#include <Eigen/Cholesky>
+#include <igl/edge_lengths.h>
+#include "massmatrix.h"
+#include "cotmatrix.h"
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -7,6 +11,22 @@ void smooth(
     double lambda,
     Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = G;
+  Eigen::MatrixXd l;
+  Eigen::SparseMatrix<double> L,A;
+  Eigen::DiagonalMatrix<double,Eigen::Dynamic> M;
+  igl::edge_lengths(V, F, l);
+  cotmatrix(l,F,L);
+  massmatrix(l,F,M);
+  Eigen::MatrixXd b=M*G;
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> list;
+  Eigen::VectorXd M1=M.diagonal();
+  for (int i=0; i<M1.size(); i++)
+    list.push_back(T(i,i,M1(i)));
+  A.resize(M.rows(),M.rows());
+  A.setFromTriplets(list.begin(), list.end()); 
+  A-=lambda*L;
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
+  solver.compute(A);
+  U=solver.solve(b);
 }