From 1cbd20d92caa678bd11637265846eddafa952053 Mon Sep 17 00:00:00 2001 From: AdrianShe Date: Sun, 14 Oct 2018 18:50:16 -0400 Subject: [PATCH 1/4] done --- src/cotmatrix.cpp | 42 ++++++++++++++++++++++++++++++++++++++++++ src/massmatrix.cpp | 17 +++++++++++++++++ src/smooth.cpp | 29 ++++++++++++++++++++++++++++- 3 files changed, 87 insertions(+), 1 deletion(-) diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp index c888e5c..911282d 100644 --- a/src/cotmatrix.cpp +++ b/src/cotmatrix.cpp @@ -1,4 +1,7 @@ #include "cotmatrix.h" +#include "igl/doublearea.h" +#include +#include void cotmatrix( const Eigen::MatrixXd & l, @@ -6,5 +9,44 @@ void cotmatrix( Eigen::SparseMatrix & L) { // Add your code here + + // Compute areas of each face + Eigen::VectorXd A; + igl::doublearea(l, A); + + // Loop through faces to determine cotangent matrix elements + // Using Eigen's functionality that triplets at duplicate (i, j) + // are summed up when creating a sparse matrix + typedef Eigen::Triplet T; + std::vector triples; + triples.reserve(F.rows() * 12); + for (int i = 0; i < F.rows(); i++) { + // Edge e = (F(j), F(j+1)) is opposite of F(j+2) mod 3 + // cotangent of the opposite angle is (-e^2 + l_1^2 + l_2^2) / (4 A) + // where l_1, l_2 are the other sides and A is the area of current face. + double len_01 = ((pow(l(i, 0), 2) + pow(l(i,1), 2) - pow(l(i,2), 2))) / (4 * A(i)); + double len_12 = ((pow(l(i, 1), 2) + pow(l(i,2), 2) - pow(l(i,0), 2))) / (4 * A(i)); + double len_02 = ((pow(l(i, 0), 2) + pow(l(i,2), 2) - pow(l(i,1), 2))) / (4 * A(i)); + double entry_01 = len_01 * 0.5; + double entry_12 = len_12 * 0.5; + double entry_02 = len_02 * 0.5; + // Insert entries + triples.push_back(T(F(i,0), F(i,1), entry_01)); + triples.push_back(T(F(i,1), F(i,0), entry_01)); + triples.push_back(T(F(i,0), F(i,0), -entry_01)); + triples.push_back(T(F(i,1), F(i,1), -entry_01)); + + triples.push_back(T(F(i,0), F(i,2), entry_02)); + triples.push_back(T(F(i,2), F(i,0), entry_02)); + triples.push_back(T(F(i,0), F(i,0), -entry_02)); + triples.push_back(T(F(i,2), F(i,2), -entry_02)); + + triples.push_back(T(F(i,1), F(i,2), entry_12)); + triples.push_back(T(F(i,2), F(i,1), entry_12)); + triples.push_back(T(F(i,1), F(i,1), -entry_12)); + triples.push_back(T(F(i,2), F(i,2), -entry_12)); + } + L.resize(F.maxCoeff() + 1, F.maxCoeff() + 1); + L.setFromTriplets(triples.begin(), triples.end()); } diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp index 577e90f..6287b92 100644 --- a/src/massmatrix.cpp +++ b/src/massmatrix.cpp @@ -1,4 +1,5 @@ #include "massmatrix.h" +#include "igl/doublearea.h" void massmatrix( const Eigen::MatrixXd & l, @@ -6,5 +7,21 @@ void massmatrix( Eigen::DiagonalMatrix & M) { // Add your code here + + // Compute area of each face + Eigen::VectorXd A; + igl::doublearea(l, A); + M.resize(F.maxCoeff() + 1); + + // Mass matrix entry (i, i) is 1/3 * (sum of faces containing vertex i) + for (int i = 0; i <= F.maxCoeff(); i++) { + double entry = 0; + for (int j = 0; j < F.rows(); j++) { + if (F(j, 0) == i || F(j,1) == i || F(j, 2) == i) { + entry = entry + A(j); + } + } + M.diagonal()[i] = entry / 3; + } } diff --git a/src/smooth.cpp b/src/smooth.cpp index 6f72147..2967585 100644 --- a/src/smooth.cpp +++ b/src/smooth.cpp @@ -1,4 +1,10 @@ #include "smooth.h" +#include "cotmatrix.h" +#include "massmatrix.h" +#include "igl/edge_lengths.h" +#include "igl/cotmatrix.h" +#include "Eigen/SparseCholesky" +#include void smooth( const Eigen::MatrixXd & V, @@ -8,5 +14,26 @@ void smooth( Eigen::MatrixXd & U) { // Replace with your code - U = G; + + // Compute edge lengths + Eigen::MatrixXd l; + igl::edge_lengths(V, F, l); + + // Construct cot and mass matrices + Eigen::SparseMatrix L_real; + Eigen::SparseMatrix L_comp; + Eigen::DiagonalMatrix M; + cotmatrix(l, F, L_comp); + massmatrix(l, F, M); + + // Solve for U as solution of M G = (M - lambda * L) U; + Eigen::SparseMatrix I(F.maxCoeff() + 1, F.maxCoeff() + 1); + I.setIdentity(); + Eigen::SparseMatrix L_lambda = I - lambda * M.inverse() * L_comp; + L_lambda = M * L_lambda; + Eigen::SimplicialLDLT> solver; + solver.compute(L_lambda); + for (int i = 0; i < G.cols(); i++) { + U.col(i) = solver.solve(M * G.col(i)); + } } From 54f90325eadfa9bb72d780077613abc69a01ec52 Mon Sep 17 00:00:00 2001 From: AdrianShe Date: Sun, 14 Oct 2018 18:52:15 -0400 Subject: [PATCH 2/4] get rid of stuff for debugging --- src/cotmatrix.cpp | 1 - src/smooth.cpp | 2 -- 2 files changed, 3 deletions(-) diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp index 911282d..c87d4c8 100644 --- a/src/cotmatrix.cpp +++ b/src/cotmatrix.cpp @@ -1,6 +1,5 @@ #include "cotmatrix.h" #include "igl/doublearea.h" -#include #include void cotmatrix( diff --git a/src/smooth.cpp b/src/smooth.cpp index 2967585..ef4c786 100644 --- a/src/smooth.cpp +++ b/src/smooth.cpp @@ -2,9 +2,7 @@ #include "cotmatrix.h" #include "massmatrix.h" #include "igl/edge_lengths.h" -#include "igl/cotmatrix.h" #include "Eigen/SparseCholesky" -#include void smooth( const Eigen::MatrixXd & V, From fe25c82baf325afcc255af6ae6c4d1206024465c Mon Sep 17 00:00:00 2001 From: AdrianShe Date: Sun, 14 Oct 2018 18:55:03 -0400 Subject: [PATCH 3/4] get rid of stuff for debugging --- src/cotmatrix.cpp | 6 +++--- src/smooth.cpp | 5 ++--- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp index c87d4c8..00e6add 100644 --- a/src/cotmatrix.cpp +++ b/src/cotmatrix.cpp @@ -23,9 +23,9 @@ void cotmatrix( // Edge e = (F(j), F(j+1)) is opposite of F(j+2) mod 3 // cotangent of the opposite angle is (-e^2 + l_1^2 + l_2^2) / (4 A) // where l_1, l_2 are the other sides and A is the area of current face. - double len_01 = ((pow(l(i, 0), 2) + pow(l(i,1), 2) - pow(l(i,2), 2))) / (4 * A(i)); - double len_12 = ((pow(l(i, 1), 2) + pow(l(i,2), 2) - pow(l(i,0), 2))) / (4 * A(i)); - double len_02 = ((pow(l(i, 0), 2) + pow(l(i,2), 2) - pow(l(i,1), 2))) / (4 * A(i)); + double len_01 = ((pow(l(i, 0), 2) + pow(l(i,1), 2) - pow(l(i,2), 2))) / (4.0 * A(i)); + double len_12 = ((pow(l(i, 1), 2) + pow(l(i,2), 2) - pow(l(i,0), 2))) / (4.0 * A(i)); + double len_02 = ((pow(l(i, 0), 2) + pow(l(i,2), 2) - pow(l(i,1), 2))) / (4.0 * A(i)); double entry_01 = len_01 * 0.5; double entry_12 = len_12 * 0.5; double entry_02 = len_02 * 0.5; diff --git a/src/smooth.cpp b/src/smooth.cpp index ef4c786..1b3a8ab 100644 --- a/src/smooth.cpp +++ b/src/smooth.cpp @@ -18,8 +18,7 @@ void smooth( igl::edge_lengths(V, F, l); // Construct cot and mass matrices - Eigen::SparseMatrix L_real; - Eigen::SparseMatrix L_comp; + Eigen::SparseMatrix L; Eigen::DiagonalMatrix M; cotmatrix(l, F, L_comp); massmatrix(l, F, M); @@ -27,7 +26,7 @@ void smooth( // Solve for U as solution of M G = (M - lambda * L) U; Eigen::SparseMatrix I(F.maxCoeff() + 1, F.maxCoeff() + 1); I.setIdentity(); - Eigen::SparseMatrix L_lambda = I - lambda * M.inverse() * L_comp; + Eigen::SparseMatrix L_lambda = I - lambda * M.inverse() * L; L_lambda = M * L_lambda; Eigen::SimplicialLDLT> solver; solver.compute(L_lambda); From 0bc527186a82fafd041948e70a9424db74f6a7a8 Mon Sep 17 00:00:00 2001 From: AdrianShe Date: Sun, 14 Oct 2018 18:59:37 -0400 Subject: [PATCH 4/4] oops --- src/smooth.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smooth.cpp b/src/smooth.cpp index 1b3a8ab..9c606ca 100644 --- a/src/smooth.cpp +++ b/src/smooth.cpp @@ -20,7 +20,7 @@ void smooth( // Construct cot and mass matrices Eigen::SparseMatrix L; Eigen::DiagonalMatrix M; - cotmatrix(l, F, L_comp); + cotmatrix(l, F, L); massmatrix(l, F, M); // Solve for U as solution of M G = (M - lambda * L) U;