precomputation done

include/igl/direct_delta_mush.cpp

@@ -10,7 +10,6 @@
 #include "diag.h"
 
 
-// ===== DEBUG: START
 #include <Eigen/Geometry>
 #include <Eigen/Sparse>
 #include <Eigen/Dense>
@@ -21,10 +20,9 @@
 #include <iostream>
 #include "Timer.h"
 
-using namespace std;
-// ===== DEBUG: END
-
-// TODOs:
+// TODOs
+// 1. U_precomputed, Psi, Omega should be #V by 10 instead of #V by 16
+// 2. Vectorize Psi computation.
 
 template <
   typename DerivedV,
@@ -36,14 +34,16 @@ template <
   typename DerivedTAlloc,
   typename DerivedU>
 IGL_INLINE void igl::direct_delta_mush(
-  const Eigen::MatrixBase<DerivedV> &V,
-  const Eigen::MatrixBase<DerivedF> &F,
-  const Eigen::MatrixBase<DerivedC> &C,
-  const Eigen::MatrixBase<DerivedE> &E,
-  const Eigen::SparseMatrix<DerivedW> &W,
-  const std::vector<DerivedT, DerivedTAlloc> &T,
-  Eigen::PlainObjectBase<DerivedU> &U)
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedC> & C,
+  const Eigen::MatrixBase<DerivedE> & E,
+  const Eigen::SparseMatrix<DerivedW> & W,
+  const std::vector<DerivedT, DerivedTAlloc> & T,
+  Eigen::PlainObjectBase<DerivedU> & U)
 {
+  using namespace std;
+  using namespace Eigen;
   cout << "START DDM" << endl;
   cout << "END DDM" << endl;
 }
@@ -56,19 +56,21 @@ template <
   typename DerivedW,
   typename DerivedOmega>
 IGL_INLINE void igl::direct_delta_mush_precomputation(
-  const Eigen::MatrixBase<DerivedV> &V,
-  const Eigen::MatrixBase<DerivedF> &F,
-  const Eigen::MatrixBase<DerivedC> &C,
-  const Eigen::MatrixBase<DerivedE> &E,
-  const Eigen::SparseMatrix<DerivedW> &W,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedC> & C,
+  const Eigen::MatrixBase<DerivedE> & E,
+  const Eigen::SparseMatrix<DerivedW> & W,
   const int p,
   const typename DerivedV::Scalar lambda,
   const typename DerivedV::Scalar kappa,
   const typename DerivedV::Scalar alpha,
-  Eigen::PlainObjectBase<DerivedOmega> &Omega)
+  Eigen::PlainObjectBase<DerivedOmega> & Omega)
 {
+  using namespace std;
+  using namespace Eigen;
   assert(kappa < lambda &&
-    "kappa needs to be smaller than lambda so that optimization for R_i is well defined");
+         "kappa needs to be smaller than lambda so that optimization for R_i is well defined");
   cout << "START DDM Precomputation" << endl;
   cout << "Using params:"
        << "\np: " << p
@@ -83,7 +85,7 @@ IGL_INLINE void igl::direct_delta_mush_precomputation(
        << endl;
 
   const int n = V.rows();
-  const int m = C.rows();
+  const int m = E.rows();
 
   // U: 4 by #V homogeneous transposed version of V
   Eigen::MatrixXd U(4, n);
@@ -91,7 +93,7 @@ IGL_INLINE void igl::direct_delta_mush_precomputation(
   Eigen::VectorXd ones(n);
   for (int i = 0; i < n; i++)
   {
-      ones(i) = 1;
+    ones(i) = 1;
   }
   U.row(U.rows() - 1) = ones;
   cout << "U: " << U.rows() << " x " << U.cols() << " Sum: " << U.sum() << endl;
@@ -113,7 +115,7 @@ IGL_INLINE void igl::direct_delta_mush_precomputation(
   Eigen::SparseMatrix<double> D_L_sparse = D_L.sparseView();
   cout << "D_L_sparse: " << D_L_sparse.rows() << " x " << D_L_sparse.cols() << " Sum: " << D_L_sparse.sum() << endl;
   Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
-  ldlt.compute(D_L_sparse);
+  ldlt.compute(D_L_sparse); // this will take quite a while to compute
   Eigen::SparseMatrix<double> L_bar = ldlt.solve(L).transpose();
   cout << "L_bar: " << L_bar.rows() << " x " << L_bar.cols() << " Sum: " << L_bar.sum() << endl;
 
@@ -127,14 +129,15 @@ IGL_INLINE void igl::direct_delta_mush_precomputation(
   cout << "W_prime: " << W_prime.rows() << " x " << W_prime.cols() << " Sum: " << W_prime.sum() << endl;
   ldlt_W_prime.compute(C_calc);
   cout << "computing W'" << endl;
-  for (int iter = 0; iter < p; iter++) {
+  for (int iter = 0; iter < p; iter++)
+  {
     cout << "iter:" << iter << endl;
     W_prime.makeCompressed();
     W_prime = ldlt_W_prime.solve(W_prime);
   }
   cout << "W_prime: " << W_prime.rows() << " x " << W_prime.cols() << " Sum: " << W_prime.sum() << endl;
 
-  // Psi was hard to solve iteratively since i couldnt express u_k \times u_k^T as matrix form
+  // Psi was hard to solve iteratively since i couldn't express u_k \times u_k^T as matrix form
   Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt_U_tilde;
   Eigen::SparseMatrix<double> B_calc((I + lambda * L_bar).transpose());
   cout << "B_calc: " << B_calc.rows() << " x " << B_calc.cols() << " Sum: " << B_calc.sum() << endl;
@@ -143,7 +146,8 @@ IGL_INLINE void igl::direct_delta_mush_precomputation(
   Eigen::SparseMatrix<double> U_tilde(U_sparse.transpose());
   cout << "U_tilde: " << U_tilde.rows() << " x " << U_tilde.cols() << " Sum: " << U_tilde.sum() << endl;
   cout << "computing U_tilde'" << endl;
-  for (int i = 0; i < p; i++) {
+  for (int i = 0; i < p; i++)
+  {
     cout << "i = " << i << endl;
     U_tilde.makeCompressed();
     U_tilde = ldlt_U_tilde.solve(U_tilde);
@@ -151,45 +155,60 @@ IGL_INLINE void igl::direct_delta_mush_precomputation(
   U_tilde = U_tilde.transpose();
   cout << "U_tilde: " << U_tilde.rows() << " x " << U_tilde.cols() << " Sum: " << U_tilde.sum() << endl;
 
-  // TODO: sparse didn't work here
-  Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> qr_B;
+  // TODO: sparse didn't work here - malloc error
   Eigen::MatrixXd U_dense(U_sparse);
   Eigen::MatrixXd U_tilde_dense(U_tilde);
   Eigen::MatrixXd B_inv_dense = U_dense.householderQr().solve(U_tilde_dense);
   Eigen::SparseMatrix<double> B_inv = B_inv_dense.sparseView();
+  // Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> qr_B;
   // U_sparse.makeCompressed();
   // Eigen::SparseMatrix<double> V_sparse_transpose(U_sparse), V_tilde_transpose(U_tilde);
   // qr_B.compute(V_sparse_transpose);
   // Eigen::SparseMatrix<double> B_inv = qr_B.solve(V_tilde_transpose);
   cout << "B_inv: " << B_inv.rows() << " x " << B_inv.cols() << " Sum: " << B_inv.sum() << endl;
 
+  // U_precomputed: #V by 10
+  // U_precomputed.row(i) = u_i \dot u_i^T \in R^{4 x 4}
+  Eigen::MatrixXd U_precomputed(n, 16);
+  for (int k = 0; k < n; k++)
+  {
+    Eigen::Vector4d u = U.col(k);
+    Eigen::MatrixXd u_full = U.col(k) * U.col(k).transpose();
+    Eigen::VectorXd u_full_vector(Eigen::Map<Eigen::VectorXd>(u_full.data(), u_full.cols() * u_full.rows()));
+    U_precomputed.row(k) = u_full_vector;
+  }
+  cout << "U_precomputed: " << U_precomputed.rows() << " x " << U_precomputed.cols() << " Sum: " << U_precomputed.sum()
+       << endl;
+
   // Psi: #V * #C by 16 (10) of \Psi_{ij}s.
   // To access \Psi_{ij}, look at row (i * m + j)
-  Eigen::MatrixXd Psi(n * m, 16); // TODO: this could be 10 to reduce storage
-
-  // TODO: FIXME: VECTORIZE THIS
+  // this takes even longer since it is not vectorized
+  Eigen::MatrixXd Psi(n * m, 16);
   for (int i = 0; i < n; i++)
   {
-    cout << "i = " << i << " / n = " << n << endl;
     for (int j = 0; j < m; j++)
     {
-      Eigen::MatrixXd Psi_curr = Eigen::MatrixXd::Zero(4, 4);
+      Eigen::VectorXd Psi_curr = Eigen::VectorXd::Zero(16);
       for (int k = 0; k < n; k++)
       {
-        Psi_curr += B_inv.coeff(k, i) * W.coeff(k, j) * U.col(k) * U.col(k).transpose();
+        // FAILED at k = 0, j = 4 since
+        // - in the paper, W.shape = #V by #C (num bones)
+        // - in our codebase, W.shape = #V by #E (num joints)
+        Psi_curr += B_inv.coeff(k, i) * W.coeff(k, j) * U_precomputed.row(k);
+//        Psi_curr += B_inv.coeff(k, i) * U_precomputed.row(k);
       }
-      Eigen::VectorXd Psi_vector(Eigen::Map<Eigen::VectorXd>(Psi_curr.data(), Psi.cols()*Psi.rows()));
-      Psi.row(i * m + j) = Psi_vector;
+      Psi.row(i * m + j) = Psi_curr;
     }
   }
   cout << "Psi: " << Psi.rows() << " x " << Psi.cols() << " Sum: " << Psi.sum() << endl;
 
-  // vector ps
+  // vector p_i:
+  // the sum of 3th, 7th and 11th of the stored 16-float vector representing Psis
+  // from j = 1 to m.
   Eigen::MatrixXd P_vectors(n, 3);
   for (int i = 0; i < n; i++)
   {
     Eigen::Vector3d p_i = Eigen::Vector3d::Zero(3);
-    // Eigen::seq not usable???
     for (int j = 0; j < m; j++)
     {
       Eigen::Vector3d p_i_curr(3);
@@ -198,26 +217,29 @@ IGL_INLINE void igl::direct_delta_mush_precomputation(
     }
     P_vectors.row(i) = p_i;
   }
+  cout << "P_vectors: " << P_vectors.rows() << " x " << P_vectors.cols() << " Sum: " << P_vectors.sum() << endl;
 
   // Omega
-  Omega.resize(n, 16);
+  Omega.resize(n * m, 16);
   for (int i = 0; i < n; i++)
   {
-    Eigen::Matrix3d p_matrix(4, 4);
-    Eigen::Vector3d curr_p = P_vectors.row(i);
+    Eigen::MatrixXd p_matrix(4, 4);
+    Eigen::VectorXd curr_p = P_vectors.row(i);
     p_matrix.block(0, 0, 3, 3) = curr_p * curr_p.transpose();
-    p_matrix.block(3, 0, 3, 1) = curr_p.transpose();
-    p_matrix.block(0, 3, 3, 3) = curr_p;
+    p_matrix.block(3, 0, 1, 3) = curr_p.transpose();
+    p_matrix.block(0, 3, 3, 1) = curr_p;
     p_matrix(3, 3) = 1;
     for (int j = 0; j < m; j++)
     {
       Eigen::MatrixXd Omega_curr(4, 4);
-      Omega_curr = (1 - alpha) * Psi.block(i * m + j, 0, 1, 16).eval() + alpha * W_prime.coeff(i, j) * p_matrix;
-      Eigen::VectorXd Omega_vector(Eigen::Map<Eigen::VectorXd>(Omega_curr.data(), Omega.cols()*Omega.rows()));
+      Eigen::MatrixXd Psi_curr = Psi.block(i * m + j, 0, 1, 16);
+      Psi_curr.resize(4, 4);
+      Omega_curr = (1 - alpha) * Psi_curr + alpha * W_prime.coeff(i, j) * p_matrix;
+      Eigen::VectorXd Omega_vector(Eigen::Map<Eigen::VectorXd>(Omega_curr.data(), Omega_curr.cols() * Omega_curr.rows()));
       Omega.row(i * m + j) = Omega_vector;
     }
   }
-
+  cout << "Omega: " << Omega.rows() << " x " << Omega.cols() << " Sum: " << Omega.sum() << endl;
 
   cout << "END DDM Precomputation" << endl;
 }
@@ -228,17 +250,42 @@ template <
   typename DerivedOmega,
   typename DerivedU>
 IGL_INLINE void igl::direct_delta_mush_pose_evaluation(
-  const std::vector<DerivedT, DerivedTAlloc> &T,
-  const Eigen::MatrixBase<DerivedOmega> &Omega,
-  Eigen::PlainObjectBase<DerivedU> &U)
+  const std::vector<DerivedT, DerivedTAlloc> & T,
+  const Eigen::MatrixBase<DerivedOmega> & Omega,
+  Eigen::PlainObjectBase<DerivedU> & U)
 {
+  using namespace std;
+  using namespace Eigen;
   cout << "START DDM Pose Eval" << endl;
   cout << "END DDM Pose Eval" << endl;
 }
 
 #ifdef IGL_STATIC_LIBRARY
+
 // Explicit template instantiation
-template void igl::direct_delta_mush<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Transform<double, 3, 2, 0>, Eigen::aligned_allocator<Eigen::Transform<double, 3, 2, 0> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int> const&, std::__1::vector<Eigen::Transform<double, 3, 2, 0>, Eigen::aligned_allocator<Eigen::Transform<double, 3, 2, 0> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::direct_delta_mush_precomputation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int> const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::direct_delta_mush_pose_evaluation<Eigen::Transform<double, 3, 2, 0>, Eigen::aligned_allocator<Eigen::Transform<double, 3, 2, 0> >, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::__1::vector<Eigen::Transform<double, 3, 2, 0>, Eigen::aligned_allocator<Eigen::Transform<double, 3, 2, 0> > > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void
+igl::direct_delta_mush<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Transform<double, 3, 2, 0>, Eigen::aligned_allocator<Eigen::Transform<double, 3, 2, 0> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(
+  Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const &,
+  Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const &,
+  Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const &,
+  Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const &, Eigen::SparseMatrix<double, 0, int> const &,
+  std::__1::vector<Eigen::Transform<double, 3, 2, 0>, Eigen::aligned_allocator<Eigen::Transform<double, 3, 2, 0> > > const &,
+  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > &);
+
+template void
+igl::direct_delta_mush_precomputation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(
+  Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const &,
+  Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const &,
+  Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const &,
+  Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const &, Eigen::SparseMatrix<double, 0, int> const &, int,
+  Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar,
+  Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar,
+  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > &);
+
+template void
+igl::direct_delta_mush_pose_evaluation<Eigen::Transform<double, 3, 2, 0>, Eigen::aligned_allocator<Eigen::Transform<double, 3, 2, 0> >, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(
+  std::__1::vector<Eigen::Transform<double, 3, 2, 0>, Eigen::aligned_allocator<Eigen::Transform<double, 3, 2, 0> > > const &,
+  Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const &,
+  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > &);
+
 #endif

include/igl/direct_delta_mush.h

@@ -38,13 +38,13 @@ namespace igl {
     typename DerivedTAlloc,
     typename DerivedU>
   IGL_INLINE void direct_delta_mush(
-    const Eigen::MatrixBase<DerivedV> &V,
-    const Eigen::MatrixBase<DerivedF> &F,
-    const Eigen::MatrixBase<DerivedC> &C,
-    const Eigen::MatrixBase<DerivedE> &E,
-    const Eigen::SparseMatrix<DerivedW> &W,
-    const std::vector<DerivedT, DerivedTAlloc> &T,
-    Eigen::PlainObjectBase<DerivedU> &U);
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedC> & C,
+    const Eigen::MatrixBase<DerivedE> & E,
+    const Eigen::SparseMatrix<DerivedW> & W,
+    const std::vector<DerivedT, DerivedTAlloc> & T,
+    Eigen::PlainObjectBase<DerivedU> & U);
 
   // Precomputation
   //
@@ -52,8 +52,8 @@ namespace igl {
   //   V  #V by 3 list of rest pose vertex positions
   //   F  #F by 3 list of triangle indices into rows of V
   //   C  #C by 3 list of rest pose bone endpoint positions
-  //   E  #T by 2 list of bone edge indices into rows of C
-  //   W  #V by #C list of weights // TODO: this could be a sparse matrix?
+  //   E  #E by 2 list of bone edge indices into rows of C
+  //   W  #V by #E list of weights
   //   p  number of smoothing iterations
   //   lambda  smoothing step size
   //   kappa  smoothness parameter (section 3.3)
@@ -67,16 +67,16 @@ namespace igl {
     typename DerivedW,
     typename DerivedOmega>
   IGL_INLINE void direct_delta_mush_precomputation(
-    const Eigen::MatrixBase<DerivedV> &V,
-    const Eigen::MatrixBase<DerivedF> &F,
-    const Eigen::MatrixBase<DerivedC> &C,
-    const Eigen::MatrixBase<DerivedE> &E,
-    const Eigen::SparseMatrix<DerivedW> &W,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedC> & C,
+    const Eigen::MatrixBase<DerivedE> & E,
+    const Eigen::SparseMatrix<DerivedW> & W,
     const int p,
     const typename DerivedV::Scalar lambda,
     const typename DerivedV::Scalar kappa,
     const typename DerivedV::Scalar alpha,
-    Eigen::PlainObjectBase<DerivedOmega> &Omega);
+    Eigen::PlainObjectBase<DerivedOmega> & Omega);
 
   // Pose evaluation
   //   Omega  #V by #T*10 list of precomputated matrix values
@@ -89,9 +89,9 @@ namespace igl {
     typename DerivedOmega,
     typename DerivedU>
   IGL_INLINE void direct_delta_mush_pose_evaluation(
-    const std::vector<DerivedT, DerivedTAlloc> &T,
-    const Eigen::MatrixBase<DerivedOmega> &Omega,
-    Eigen::PlainObjectBase<DerivedU> &U);
+    const std::vector<DerivedT, DerivedTAlloc> & T,
+    const Eigen::MatrixBase<DerivedOmega> & Omega,
+    Eigen::PlainObjectBase<DerivedU> & U);
 } // namespace igl
 
 #ifndef IGL_STATIC_LIBRARY

tutorial/408_DirectDeltaMush/main.cpp

@@ -10,6 +10,7 @@
 #include <igl/readOBJ.h>
 #include <igl/readTGF.h>
 #include <igl/opengl/glfw/Viewer.h>
+#include <igl/embree/bone_heat.h>
 
 #include <Eigen/Geometry>
 #include <Eigen/StdVector>
@@ -41,19 +42,21 @@ float lambda = 0.5;
 float kappa = 0.4; // kappa < lambda to keep R_i well-defined
 float alpha = 0.5;
 
-bool pre_draw(igl::opengl::glfw::Viewer &viewer)
+bool pre_draw(igl::opengl::glfw::Viewer & viewer)
 {
   using namespace Eigen;
   using namespace std;
-  if (recompute) {
+  if (recompute)
+  {
     // Find pose interval
-    const int begin = (int)floor(anim_t) % poses.size();
-    const int end = (int)(floor(anim_t) + 1) % poses.size();
+    const int begin = (int) floor(anim_t) % poses.size();
+    const int end = (int) (floor(anim_t) + 1) % poses.size();
     const double t = anim_t - floor(anim_t);
 
     // Interpolate pose and identity
     RotationList anim_pose(poses[begin].size());
-    for (int e = 0; e < poses[begin].size(); e++) {
+    for (int e = 0; e < poses[begin].size(); e++)
+    {
       anim_pose[e] = poses[begin][e].slerp(t, poses[end][e]);
     }
     // Propagate relative rotations via FK to retrieve absolute transformations
@@ -63,7 +66,8 @@ bool pre_draw(igl::opengl::glfw::Viewer &viewer)
     const int dim = C.cols();
     MatrixXd T(BE.rows() * (dim + 1), dim);
     TransformationList T_list(BE.rows());
-    for (int e = 0; e < BE.rows(); e++) {
+    for (int e = 0; e < BE.rows(); e++)
+    {
       Affine3d a = Affine3d::Identity();
       a.translate(vT[e]);
       a.rotate(vQ[e]);
@@ -71,11 +75,12 @@ bool pre_draw(igl::opengl::glfw::Viewer &viewer)
         a.matrix().transpose().block(0, 0, dim + 1, dim);
     }
     // Compute deformation via LBS as matrix multiplication
-    if (use_ddm) {
+    if (use_ddm)
+    {
       igl::direct_delta_mush_pose_evaluation(T_list, Omega, U);
       igl::direct_delta_mush(V, F, C, BE, W_sparse, T_list, U);
-    }
-    else {
+    } else
+    {
       U = M * T;
     }
 
@@ -87,20 +92,22 @@ bool pre_draw(igl::opengl::glfw::Viewer &viewer)
     viewer.data().set_vertices(U);
     viewer.data().set_edges(CT, BET, sea_green);
     viewer.data().compute_normals();
-    if (viewer.core().is_animating) {
+    if (viewer.core().is_animating)
+    {
       anim_t += anim_t_dir;
-    }
-    else {
+    } else
+    {
       recompute = false;
     }
   }
   return false;
 }
 
-bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
+bool key_down(igl::opengl::glfw::Viewer & viewer, unsigned char key, int mods)
 {
   recompute = true;
-  switch (key) {
+  switch (key)
+  {
     case 'D':
     case 'd':
       use_ddm = !use_ddm;
@@ -135,7 +142,6 @@ int main(int argc, char *argv[])
   igl::lbs_matrix(V, W, M);
 
   // Precomputation for Direct Delta Mush
-  cout<<"Initializing Direct Delta Mush..."<<endl;
   igl::direct_delta_mush_precomputation(V, F, C, BE, W_sparse, p, lambda, kappa, alpha, Omega);
 
   // Plot the mesh with pseudocolors