precomputation - W

cmake/LibiglFolders.cmake

@@ -82,6 +82,7 @@ igl_folder_targets("Tutorials"
     405_AsRigidAsPossible_bin
     406_FastAutomaticSkinningTransformations_bin
     407_BiharmonicCoordinates_bin
+    408_DirectDeltaMush_bin
     501_HarmonicParam_bin
     502_LSCMParam_bin
     503_ARAPParam_bin

include/igl/direct_delta_mush.cpp

@@ -1,4 +1,4 @@
-// This file is part of libigl, a simple c++ geometry processing library.
+// This file is part of libigl, a simple C++ geometry processing library.
 //
 // Copyright (C) 2020 Xiangyu Kong <[email protected]>
 //
@@ -13,14 +13,16 @@
 // ===== DEBUG: START
 #include <Eigen/Geometry>
 #include <Eigen/Sparse>
+#include <Eigen/Dense>
+#include <Eigen/Cholesky>
 #include <Eigen/LU>
 #include <iostream>
+#include "Timer.h"
 
 using namespace std;
 // ===== DEBUG: END
 
 // TODOs:
-// 1. All `Scalar`s are temporarily substituted as `DerivedV`
 
 template <
   typename DerivedV,
@@ -36,7 +38,7 @@ IGL_INLINE void igl::direct_delta_mush(
   const Eigen::MatrixBase<DerivedF> &F,
   const Eigen::MatrixBase<DerivedC> &C,
   const Eigen::MatrixBase<DerivedE> &E,
-  const Eigen::MatrixBase<DerivedW> &W,
+  const Eigen::SparseMatrix<DerivedW> &W,
   const std::vector<DerivedT, DerivedTAlloc> &T,
   Eigen::PlainObjectBase<DerivedU> &U)
 {
@@ -56,13 +58,65 @@ IGL_INLINE void igl::direct_delta_mush_precomputation(
   const Eigen::MatrixBase<DerivedF> &F,
   const Eigen::MatrixBase<DerivedC> &C,
   const Eigen::MatrixBase<DerivedE> &E,
-  const Eigen::MatrixBase<DerivedW> &W,
+  const Eigen::SparseMatrix<DerivedW> &W,
   const int p,
   const typename DerivedV::Scalar lambda,
   const typename DerivedV::Scalar kappa,
   Eigen::PlainObjectBase<DerivedOmega> &Omega)
 {
+  assert(kappa < lambda &&
+    "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
+       << "\nlambda: " << lambda
+       << "\nkappa: " << kappa
+       << endl;
+  cout << "V: " << V.rows() << " x " << V.cols() << " Sum: " << V.sum()
+       << "\nF: " << F.rows() << " x " << F.cols() << " Sum: " << F.sum()
+       << "\nC: " << C.rows() << " x " << C.cols() << " Sum: " << C.sum()
+       << "\nE: " << E.rows() << " x " << E.cols() << " Sum: " << E.sum()
+       << "\nW: " << W.rows() << " x " << W.cols() << " Sum: " << W.sum()
+       << endl;
+
+  const int n = V.rows();
+  const int m = C.rows();
+
+  // Identity of #V by #V
+  Eigen::SparseMatrix<double> I(n, n);
+  I.setIdentity();
+  cout << "I: " << I.rows() << " x " << I.cols() << " Sum: " << I.sum() << endl;
+
+  // Laplacian: L_bar = L \times D_L^{-1}
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix(V, F, L);
+  cout << "L: " << L.rows() << " x " << L.cols() << " Sum: " << L.sum() << endl;
+  Eigen::MatrixXd D_L = L.diagonal().asDiagonal();
+  cout << "D_L: " << D_L.rows() << " x " << D_L.cols() << " Sum: " << D_L.sum() << endl;
+  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);
+  Eigen::SparseMatrix<double> L_bar = ldlt.solve(L).transpose();
+  cout << "L_bar: " << L_bar.rows() << " x " << L_bar.cols() << " Sum: " << L_bar.sum() << endl;
+
+  // Implicitly and iteratively solve
+  // w'_{ij} = \sum_{k=1}^{n}{C_{ki} w_{kj}}, C = (I + kappa L_bar)^{-p}:
+  // W' = C^T \times W
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt_C;
+  Eigen::SparseMatrix<double> C_calc((I + kappa * L_bar).transpose());
+  cout << "C_calc: " << C_calc.rows() << " x " << C_calc.cols() << " Sum: " << C_calc.sum() << endl;
+  Eigen::SparseMatrix<double> W_prime(W), W_next;
+  cout << "W_prime: " << W_prime.rows() << " x " << W_prime.cols() << " Sum: " << W_prime.sum() << endl;
+  ldlt_C.compute(C_calc);
+  cout << "computing W'" << endl;
+  for (int iter = 0; iter < p; iter++) {
+    cout << "iter:" << iter << endl;
+    W_prime.makeCompressed();
+    W_prime = ldlt_C.solve(W_prime);
+  }
+  cout << "W_prime: " << W_prime.rows() << " x " << W_prime.cols() << " Sum: " << W_prime.sum() << endl;
+
   cout << "END DDM Precomputation" << endl;
 }
 
@@ -82,6 +136,7 @@ IGL_INLINE void igl::direct_delta_mush_pose_evaluation(
 
 #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>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > 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>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, 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<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::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

@@ -11,9 +11,9 @@
 #include "igl_inline.h"
 
 #include <Eigen/Core>
-// #include <Eigen/Geometry>
+#include <Eigen/Geometry>
 #include <Eigen/Sparse>
-// #include <Eigen/LU>
+#include <Eigen/LU>
 #include <vector>
 
 namespace igl {
@@ -42,7 +42,7 @@ namespace igl {
     const Eigen::MatrixBase<DerivedF> &F,
     const Eigen::MatrixBase<DerivedC> &C,
     const Eigen::MatrixBase<DerivedE> &E,
-    const Eigen::MatrixBase<DerivedW> &W,
+    const Eigen::SparseMatrix<DerivedW> &W,
     const std::vector<DerivedT, DerivedTAlloc> &T,
     Eigen::PlainObjectBase<DerivedU> &U);
 
@@ -71,7 +71,7 @@ namespace igl {
     const Eigen::MatrixBase<DerivedF> &F,
     const Eigen::MatrixBase<DerivedC> &C,
     const Eigen::MatrixBase<DerivedE> &E,
-    const Eigen::MatrixBase<DerivedW> &W,
+    const Eigen::SparseMatrix<DerivedW> &W,
     const int p,
     const typename DerivedV::Scalar lambda,
     const typename DerivedV::Scalar kappa,

tutorial/408_DirectDeltaMush/main.cpp

@@ -13,6 +13,7 @@
 
 #include <Eigen/Geometry>
 #include <Eigen/StdVector>
+#include <Eigen/Sparse>
 #include <vector>
 #include <algorithm>
 #include <iostream>
@@ -29,14 +30,15 @@ const Eigen::RowVector3d sea_green(70. / 255., 252. / 255., 167. / 255.);
 Eigen::MatrixXd V, W, C, U, M, Omega;
 Eigen::MatrixXi F, BE;
 Eigen::VectorXi P;
+Eigen::SparseMatrix<double> W_sparse;
 std::vector<RotationList> poses;
 double anim_t = 0.0;
 double anim_t_dir = 0.015;
-bool use_dqs = false;
+bool use_ddm = false;
 bool recompute = true;
-float p = 10;
-float lambda = 1.;
-float kappa = 1.;
+float p = 3.;
+float lambda = 0.5;
+float kappa = 0.4; // kappa < lambda to keep R_i well-defined
 
 bool pre_draw(igl::opengl::glfw::Viewer &viewer)
 {
@@ -68,8 +70,9 @@ 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_dqs) {
-      igl::direct_delta_mush(V, F, C, BE, W, T_list, U);
+    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 {
       U = M * T;
@@ -99,7 +102,7 @@ bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
   switch (key) {
     case 'D':
     case 'd':
-      use_dqs = !use_dqs;
+      use_ddm = !use_ddm;
       return true;
     case ' ':
       viewer.core().is_animating = !viewer.core().is_animating;
@@ -127,11 +130,12 @@ int main(int argc, char *argv[])
   poses[3][2] = rest_pose[2] * bend * rest_pose[2].conjugate();
 
   igl::readDMAT(TUTORIAL_SHARED_PATH "/arm-weights.dmat", W);
+  W_sparse = W.sparseView();
   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, p, lambda, kappa, Omega);
+  igl::direct_delta_mush_precomputation(V, F, C, BE, W_sparse, p, lambda, kappa, Omega);
 
   // Plot the mesh with pseudocolors
   igl::opengl::glfw::Viewer viewer;