HW03 David Acuna

src/cotmatrix.cpp

@@ -1,10 +1,79 @@
 #include "cotmatrix.h"
+#include "iostream"
+#include <vector>
+#include <math.h>
+
+double cosineLaw(double a, double b, double c) {
+    // https://en.wikipedia.org/wiki/Law_of_cosines
+    double cosA = (-a * a + b * b + c * c) / (2 * c * b);
+    return cosA;
+}
+
+double cot(double a, double b, double c) {
+    //computing the cosA
+    double cosA = cosineLaw(a, b, c);
+    //since sin^2a + cos^2a=1.... then..
+    double sinA = sqrt(1 - cosA * cosA);
+    return cosA / sinA;
+}
+
+double cot(const Eigen::MatrixXd &l, int idx, int offset) {
+    //overloaded function for simplicity.
+    //columns correspond to edges [1,2],[2,0],[0,1] offset does a shift here.
+    assert(l.cols() == 3);
+    double b = l(idx, (offset + 0) % 3);
+    double c = l(idx, (offset + 1) % 3);
+    double a = l(idx, (offset + 2) % 3);
+    return cot(a, b, c);
+}
+
+void to_L_sym(Eigen::SparseMatrix<double> &L, int i, int j, double val) {
+    //This insert in symmetric fashion, and we need to sum Tb if Talpha already exists.
+    L.coeffRef(i, j) += val;
+    if (i != j)
+        L.coeffRef(j, i) += val;
+}
 
 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 N = F.maxCoeff() + 1;
+    L.resize(N, N);
+
+    //I am going to use coeffRef on the SparseMatrix for simplicity(I am going to use the same matrix as an accumulator)
+    //If speed were to be a  major concern instead of using coeffRef, we could fill the SparseMatrix with triplets as explained here:
+    //https://eigen.tuxfamily.org/dox/group__TutorialSparse.html#TutorialSparseFilling
+
+    //In anycase, let's reserve some memory before hand to avoid multiple reallocations.
+    L.reserve(F.rows() * 3 * 2);
+
+    //this should be already zero,since it is a newly initialized sparse matrix but ...
+    L.setZero();
+
+    for (int idx = 0; idx < F.rows(); idx++) {
+        int i = F(idx, 0);
+        int j = F(idx, 1);
+        int k = F(idx, 2);
+
+        double cotA = 0.5 * cot(l, idx, 0);
+        double cotB = 0.5 * cot(l, idx, 1);
+        double cotC = 0.5 * cot(l, idx, 2);
+
+        //to_L_sym-> add to existing value, symmetrically [eg. (i,j) (j,i)]
+        to_L_sym(L, i, j, cotA);
+        to_L_sym(L, j, k, cotB);
+        to_L_sym(L, k, i, cotC);
+
+    }
+
+    //diagonal
+    for (int idx = 0; idx < N; idx++) {
+        L.coeffRef(idx, idx) = -L.row(idx).sum();
+    }
+
+    L.makeCompressed(); //suppresses the remaining empty space(if any) and transforms the matrix into a compressed column storage.
+
 }
 

src/massmatrix.cpp

@@ -1,10 +1,33 @@
 #include "massmatrix.h"
+#include "igl/doublearea.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 N = F.maxCoeff() + 1;
+    M.resize(N);
+    M.setZero();
+
+    Eigen::MatrixXd areas;
+    igl::doublearea(l, areas);
+
+    //igl::doublearea computes the area of each triangle twice.
+    areas = areas / 2.0;
+
+    for (int idx = 0; idx < F.rows(); idx++) {
+        //let's get the Vs that are part of this triangle.
+        int v1 = F(idx, 0);
+        int v2 = F(idx, 1);
+        int v3 = F(idx, 2);
+
+        //the area of this triangle F(idx) contributes 1/3 to each of this vs(v1,v2,v3), so..
+        M.diagonal()(v1) += areas(idx) / 3.0;
+        M.diagonal()(v2) += areas(idx) / 3.0;
+        M.diagonal()(v3) += areas(idx) / 3.0;
+    }
+
 }
 

src/smooth.cpp

@@ -1,12 +1,41 @@
 #include "smooth.h"
+#include "igl/edge_lengths.h"
+#include "massmatrix.h"
+#include "cotmatrix.h"
+#include<Eigen/SparseCholesky>
+#include "iostream"
 
 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) {
+
+    //Computing edges lengths
+    //columns correspond to edges [1,2],[2,0],[0,1]
+    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);
+    // ---
+
+    //Let's create the spatial matrix (M-lambda*L)
+    //Unfortunately, Eigen doesnt support this operation (M-lambda*L) in one line. so..
+    Eigen::SparseMatrix<double> A = -lambda * L;
+    for (int idx = 0; idx < M.rows(); idx++) {
+        A.coeffRef(idx, idx) += M.diagonal()(idx);
+    }
+
+    // According to the documentation of eigen
+    // https://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html
+    // this is the recommended Cholesky solver for very sparse and not too large problems
+
+    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver(A);
+    U = solver.solve(M * G);
+
 }