From 4246514a96053c278cc0cd9348a0339390ee2fad Mon Sep 17 00:00:00 2001
From: Feilong Song <feilong.song@mail.utoronto.ca>
Date: Thu, 25 Oct 2018 19:43:29 -0400
Subject: [PATCH] hw4

---
 src/arap_precompute.cpp       | 34 +++++++++++++++++++++++++++++++++-
 src/arap_single_iteration.cpp | 17 ++++++++++++++++-
 src/biharmonic_precompute.cpp | 15 +++++++++++++--
 src/biharmonic_solve.cpp      |  4 ++--
 4 files changed, 64 insertions(+), 6 deletions(-)

diff --git a/src/arap_precompute.cpp b/src/arap_precompute.cpp
index a580bdb..0de1672 100644
--- a/src/arap_precompute.cpp
+++ b/src/arap_precompute.cpp
@@ -1,4 +1,8 @@
 #include "arap_precompute.h"
+#include <igl/cotmatrix.h>
+#include <igl/cotmatrix_entries.h>
+#include <igl/min_quad_with_fixed.h>
+#include <vector>
 
 void arap_precompute(
   const Eigen::MatrixXd & V,
@@ -7,5 +11,33 @@ void arap_precompute(
   igl::min_quad_with_fixed_data<double> & data,
   Eigen::SparseMatrix<double> & K)
 {
-  // REPLACE WITH YOUR CODE
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix(V, F, L);
+  igl::min_quad_with_fixed_precompute(L, b, Eigen::SparseMatrix<double>(), false, data);
+
+  Eigen::MatrixXd C(F.rows(), 3);
+  igl::cotmatrix_entries(V, F, C);
+
+  std::vector<Eigen::Triplet<double>> triplets;
+
+  K.resize(V.rows(), 3 * V.rows());
+  for (int f = 0; f < F.rows(); f++) {
+    Eigen::RowVector3i face = F.row(f);
+    for (int x = 0; x < 3; x++) {
+      int i = face((x + 1) % 3);
+      int j = face((x + 2) % 3);
+      Eigen::RowVector3d edge = V.row(i) - V.row(j);
+      Eigen::Vector3d eij = C(f, x) / 3.0 * edge ;
+
+      for (int a = 0; a < 3; a++) {
+        int k = face((x + a) % 3);
+        for (int b = 0; b < 3; b++) {
+          triplets.emplace_back(Eigen::Triplet<double>(i, 3 * k + b, eij(b)));
+          triplets.emplace_back(Eigen::Triplet<double>(j, 3 * k + b, -eij(b)));
+        }
+      }
+    }
+  }
+  K.setFromTriplets(triplets.begin(), triplets.end());
+
 }
diff --git a/src/arap_single_iteration.cpp b/src/arap_single_iteration.cpp
index 4478046..d3ebec7 100644
--- a/src/arap_single_iteration.cpp
+++ b/src/arap_single_iteration.cpp
@@ -1,4 +1,6 @@
+#include <igl/min_quad_with_fixed.h>
 #include "arap_single_iteration.h"
+#include <igl/polar_svd3x3.h>
 
 void arap_single_iteration(
   const igl::min_quad_with_fixed_data<double> & data,
@@ -6,5 +8,18 @@ void arap_single_iteration(
   const Eigen::MatrixXd & bc,
   Eigen::MatrixXd & U)
 {
-  // REPLACE WITH YOUR CODE
+  Eigen::MatrixXd C = K.transpose() * U;
+  Eigen::MatrixXd R(3 * U.rows(), 3);
+
+  for (int k = 0; k < U.rows(); k++) {
+    Eigen::Matrix3d Ck, Rk;
+    for (int i = 0; i < 3; i++) {
+      Ck.row(i) = C.row(3 * k + i);
+    }
+    igl::polar_svd3x3(Ck, Rk);
+    for (int i = 0; i < 3; i++) {
+      R.row(3 * k + i) = Rk.row(i);
+    }
+  }
+  igl::min_quad_with_fixed_solve(data, K * R, bc, Eigen::MatrixXd(), U);
 }
diff --git a/src/biharmonic_precompute.cpp b/src/biharmonic_precompute.cpp
index df69d7d..b3a11b0 100644
--- a/src/biharmonic_precompute.cpp
+++ b/src/biharmonic_precompute.cpp
@@ -1,5 +1,7 @@
 #include "biharmonic_precompute.h"
 #include <igl/min_quad_with_fixed.h>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
 
 void biharmonic_precompute(
   const Eigen::MatrixXd & V,
@@ -7,7 +9,16 @@ void biharmonic_precompute(
   const Eigen::VectorXi & b,
   igl::min_quad_with_fixed_data<double> & data)
 {
-  // REPLACE WITH YOUR CODE
-  data.n = V.rows();
+  Eigen::SparseMatrix<double> L, M;
+  igl::cotmatrix(V, F, L);
+  igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+
+  Eigen::SparseMatrix<double> M_i(M.rows(), M.cols());
+  for (int i = 0; i < M.rows(); i++) {
+    M_i.coeffRef(i, i) = 1.0 / M.coeff(i, i);
+  }
+
+  Eigen::SparseMatrix<double> Q = L.transpose() * M_i * L;
+  igl::min_quad_with_fixed_precompute(Q, b, Eigen::SparseMatrix<double>(), false, data);
 }
 
diff --git a/src/biharmonic_solve.cpp b/src/biharmonic_solve.cpp
index 8351ce1..195c0bb 100644
--- a/src/biharmonic_solve.cpp
+++ b/src/biharmonic_solve.cpp
@@ -6,7 +6,7 @@ void biharmonic_solve(
   const Eigen::MatrixXd & bc,
   Eigen::MatrixXd & D)
 {
-  // REPLACE WITH YOUR CODE
-  D = Eigen::MatrixXd::Zero(data.n,3);
+  Eigen::MatrixXd B = Eigen::MatrixXd::Zero(data.n, 1);
+  igl::min_quad_with_fixed_solve(data, B, bc, Eigen::MatrixXd(), D);
 }