Merge branch 'Dxyk-direct-delta-mush'

cmake/LibiglDownloadExternal.cmake

@@ -184,7 +184,7 @@ function(igl_download_tutorial_data)
 	igl_download_project_aux(tutorial_data
 		"${LIBIGL_EXTERNAL}/../tutorial/data"
 		GIT_REPOSITORY https://github.com/libigl/libigl-tutorial-data
-		GIT_TAG        5c6a1ea809c043d71e5595775709c15325a7158c
+		GIT_TAG        fb5fa00bc4ede64b36002d703ce541552370b6e9
 	)
 endfunction()
 

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

@@ -0,0 +1,287 @@
+// This file is part of libigl, a simple C++ geometry processing library.
+//
+// Copyright (C) 2020 Xiangyu Kong <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "direct_delta_mush.h"
+#include "cotmatrix.h"
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedOmega,
+  typename DerivedU>
+IGL_INLINE void igl::direct_delta_mush(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const std::vector<Eigen::Affine3d, Eigen::aligned_allocator<Eigen::Affine3d> > & T,
+  const Eigen::MatrixBase<DerivedOmega> & Omega,
+  Eigen::PlainObjectBase<DerivedU> & U)
+{
+  using namespace Eigen;
+
+  // Shape checks
+  assert(V.cols() == 3 && "V should contain 3D positions.");
+  assert(F.cols() == 3 && "F should contain triangles.");
+  assert(Omega.rows() == V.rows() && "Omega contain the same number of rows as V.");
+  assert(Omega.cols() == T.size() * 10 && "Omega should have #T*10 columns.");
+
+  typedef typename DerivedV::Scalar Scalar;
+
+  int n = V.rows();
+  int m = T.size();
+
+  // V_homogeneous: #V by 4, homogeneous version of V
+  // Note:
+  // In the paper, the rest pose vertices are represented in U \in R^{4 x #V}
+  // Thus the formulae involving U would differ from the paper by a transpose.
+  Matrix<Scalar, Dynamic, 4> V_homogeneous(n, 4);
+  V_homogeneous << V, Matrix<Scalar, Dynamic, 1>::Ones(n, 1);
+  U.resize(n, 3);
+
+  for (int i = 0; i < n; ++i)
+  {
+    // Construct Q matrix using Omega and Transformations
+    Matrix<Scalar, 4, 4> Q_mat(4, 4);
+    Q_mat = Matrix<Scalar, 4, 4>::Zero(4, 4);
+    for (int j = 0; j < m; ++j)
+    {
+      Matrix<typename DerivedOmega::Scalar, 4, 4> Omega_curr(4, 4);
+      Matrix<typename DerivedOmega::Scalar, 10, 1> curr = Omega.block(i, j * 10, 1, 10).transpose();
+      Omega_curr << curr(0), curr(1), curr(2), curr(3),
+        curr(1), curr(4), curr(5), curr(6),
+        curr(2), curr(5), curr(7), curr(8),
+        curr(3), curr(6), curr(8), curr(9);
+
+      Affine3d M_curr = T[j];
+      Q_mat += M_curr.matrix() * Omega_curr;
+    }
+    // Normalize so that the last element is 1
+    Q_mat /= Q_mat(Q_mat.rows() - 1, Q_mat.cols() - 1);
+
+    Matrix<Scalar, 3, 3> Q_i = Q_mat.block(0, 0, 3, 3);
+    Matrix<Scalar, 3, 1> q_i = Q_mat.block(0, 3, 3, 1);
+    Matrix<Scalar, 3, 1> p_i = Q_mat.block(3, 0, 1, 3).transpose();
+
+    // Get rotation and translation matrices using SVD
+    Matrix<Scalar, 3, 3> SVD_i = Q_i - q_i * p_i.transpose();
+    JacobiSVD<Matrix<Scalar, 3, 3>> svd;
+    svd.compute(SVD_i, ComputeFullU | ComputeFullV);
+    Matrix<Scalar, 3, 3> R_i = svd.matrixU() * svd.matrixV().transpose();
+    Matrix<Scalar, 3, 1> t_i = q_i - R_i * p_i;
+
+    // Gamma final transformation matrix
+    Matrix<Scalar, 3, 4> Gamma_i(3, 4);
+    Gamma_i.block(0, 0, 3, 3) = R_i;
+    Gamma_i.block(0, 3, 3, 1) = t_i;
+
+    // Final deformed position
+    Matrix<Scalar, 4, 1> v_i = V_homogeneous.row(i);
+    U.row(i) = Gamma_i * v_i;
+  }
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedW,
+  typename DerivedOmega>
+IGL_INLINE void igl::direct_delta_mush_precomputation(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  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)
+{
+  using namespace Eigen;
+
+  // Shape checks
+  assert(V.cols() == 3 && "V should contain 3D positions.");
+  assert(F.cols() == 3 && "F should contain triangles.");
+  assert(W.rows() == V.rows() && "W.rows() should be equal to V.rows().");
+
+  // Parameter checks
+  assert(p > 0 && "Laplacian iteration p should be positive.");
+  assert(lambda > 0 && "lambda should be positive.");
+  assert(kappa > 0 && kappa < lambda && "kappa should be positive and less than lambda.");
+  assert(alpha >= 0 && alpha < 1 && "alpha should be non-negative and less than 1.");
+
+  typedef typename DerivedV::Scalar Scalar;
+
+  // lambda helper
+  // Given a square matrix, extract the upper triangle (including diagonal) to an array.
+  // E.g. 1   2  3  4
+  //      5   6  7  8  -> [1, 2, 3, 4, 6, 7, 8, 11, 12, 16]
+  //      9  10 11 12      0  1  2  3  4  5  6   7   8   9
+  //      13 14 15 16
+  auto extract_upper_triangle = [](
+    const Matrix<Scalar, Dynamic, Dynamic> & full) -> Matrix<Scalar, Dynamic, 1>
+  {
+    int dims = full.rows();
+    Matrix<Scalar, Dynamic, 1> upper_triangle((dims * (dims + 1)) / 2);
+    int vector_idx = 0;
+    for (int i = 0; i < dims; ++i)
+    {
+      for (int j = i; j < dims; ++j)
+      {
+        upper_triangle(vector_idx) = full(i, j);
+        vector_idx++;
+      }
+    }
+    return upper_triangle;
+  };
+
+  const int n = V.rows();
+  const int m = W.cols();
+
+  // V_homogeneous: #V by 4, homogeneous version of V
+  // Note:
+  // in the paper, the rest pose vertices are represented in U \in R^{4 \times #V}
+  // Thus the formulae involving U would differ from the paper by a transpose.
+  Matrix<Scalar, Dynamic, 4> V_homogeneous(n, 4);
+  V_homogeneous << V, Matrix<Scalar, Dynamic, 1>::Ones(n);
+
+  // Identity matrix of #V by #V
+  SparseMatrix<Scalar> I(n, n);
+  I.setIdentity();
+
+  // Laplacian matrix of #V by #V
+  // L_bar = L \times D_L^{-1}
+  SparseMatrix<Scalar> L;
+  igl::cotmatrix(V, F, L);
+  L = -L;
+  // Inverse of diagonal matrix = reciprocal elements in diagonal
+  Matrix<Scalar, Dynamic, 1> D_L = L.diagonal();
+  // D_L = D_L.array().pow(-1);  // Not using this since not sure if diagonal contains 0
+  for (int i = 0; i < D_L.size(); ++i)
+  {
+    if (D_L(i) != 0)
+    {
+      D_L(i) = 1 / D_L(i);
+    }
+  }
+  SparseMatrix<Scalar> D_L_inv = D_L.asDiagonal().toDenseMatrix().sparseView();
+  SparseMatrix<Scalar> L_bar = L * D_L_inv;
+
+  // Implicitly and iteratively solve for W'
+  // w'_{ij} = \sum_{k=1}^{n}{C_{ki} w_{kj}}      where C = (I + kappa L_bar)^{-p}:
+  // W' = C^T \times W  =>  c^T W_k = W_{k-1}     where c = (I + kappa L_bar)
+  // C positive semi-definite => ldlt solver
+  SimplicialLDLT<SparseMatrix<DerivedW>> ldlt_W_prime;
+  SparseMatrix<Scalar> c(I + kappa * L_bar);
+  Matrix<DerivedW, Dynamic, Dynamic> W_prime(W);
+  ldlt_W_prime.compute(c.transpose());
+  for (int iter = 0; iter < p; ++iter)
+  {
+    W_prime = ldlt_W_prime.solve(W_prime);
+  }
+
+  // U_precomputed: #V by 10
+  // Cache u_i^T \dot u_i \in R^{4 x 4} to reduce computation time.
+  Matrix<Scalar, Dynamic, 10> U_precomputed(n, 10);
+  for (int k = 0; k < n; ++k)
+  {
+    Matrix<Scalar, 4, 4> u_full = V_homogeneous.row(k).transpose() * V_homogeneous.row(k);
+    U_precomputed.row(k) = extract_upper_triangle(u_full);
+  }
+
+  // U_prime: #V by #T*10 of u_{jx}
+  // Each column of U_prime (u_{jx}) is the element-wise product of
+  // W_j and U_precomputed_x where j \in {1...m}, x \in {1...10}
+  Matrix<Scalar, Dynamic, Dynamic> U_prime(n, m * 10);
+  for (int j = 0; j < m; ++j)
+  {
+    Matrix<Scalar, Dynamic, 1> w_j = W.col(j);
+    for (int x = 0; x < 10; ++x)
+    {
+      Matrix<Scalar, Dynamic, 1> u_x = U_precomputed.col(x);
+      U_prime.col(10 * j + x) = w_j.array() * u_x.array();
+    }
+  }
+
+  // Implicitly and iteratively solve for Psi: #V by #T*10 of \Psi_{ij}s.
+  // Note: Using dense matrices to solve for Psi will cause the program to hang.
+  // The following won't work
+  // Matrix<Scalar, Dynamic, Dynamic> Psi(U_prime);
+  // Matrix<Scalar, Dynamic, Dynamic> b((I + lambda * L_bar).transpose());
+  // for (int iter = 0; iter < p; ++iter)
+  // {
+  //   Psi = b.ldlt().solve(Psi);  // hangs here
+  // }
+  // Convert to sparse matrices and compute
+  Matrix<Scalar, Dynamic, Dynamic> Psi = U_prime.sparseView();
+  SparseMatrix<Scalar> b = (I + lambda * L_bar).transpose();
+  SimplicialLDLT<SparseMatrix<Scalar>> ldlt_Psi;
+  ldlt_Psi.compute(b);
+  for (int iter = 0; iter < p; ++iter)
+  {
+    Psi = ldlt_Psi.solve(Psi);
+  }
+
+  // P: #V by 10 precomputed upper triangle of
+  //    p_i p_i^T , p_i
+  //    p_i^T     , 1
+  // where p_i = (\sum_{j=1}^{n} Psi_{ij})'s top right 3 by 1 column
+  Matrix<Scalar, Dynamic, 10> P(n, 10);
+  for (int i = 0; i < n; ++i)
+  {
+    Matrix<Scalar, 3, 1> p_i = Matrix<Scalar, 3, 1>::Zero(3);
+    Scalar last = 0;
+    for (int j = 0; j < m; ++j)
+    {
+      Matrix<Scalar, 3, 1> p_i_curr(3);
+      p_i_curr << Psi(i, j * 10 + 3), Psi(i, j * 10 + 6), Psi(i, j * 10 + 8);
+      p_i += p_i_curr;
+      last += Psi(i, j * 10 + 9);
+    }
+    p_i /= last;  // normalize
+    Matrix<Scalar, 4, 4> p_matrix(4, 4);
+    p_matrix.block(0, 0, 3, 3) = p_i * p_i.transpose();
+    p_matrix.block(0, 3, 3, 1) = p_i;
+    p_matrix.block(3, 0, 1, 3) = p_i.transpose();
+    p_matrix(3, 3) = 1;
+    P.row(i) = extract_upper_triangle(p_matrix);
+  }
+
+  // Omega
+  Omega.resize(n, m * 10);
+  for (int i = 0; i < n; ++i)
+  {
+    Matrix<Scalar, 10, 1> p_vector = P.row(i);
+    for (int j = 0; j < m; ++j)
+    {
+      Matrix<Scalar, 10, 1> Omega_curr(10);
+      Matrix<Scalar, 10, 1> Psi_curr = Psi.block(i, j * 10, 1, 10).transpose();
+      Omega_curr = (1. - alpha) * Psi_curr + alpha * W_prime.coeff(i, j) * p_vector;
+      Omega.block(i, j * 10, 1, 10) = Omega_curr.transpose();
+    }
+  }
+}
+
+#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<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 &,
+  std::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_precomputation<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::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> > &);
+
+#endif

include/igl/direct_delta_mush.h

@@ -0,0 +1,75 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Xiangyu Kong <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_DIRECT_DELTA_MUSH_H
+#define IGL_DIRECT_DELTA_MUSH_H
+
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+#include <Eigen/Geometry>
+#include <vector>
+
+namespace igl {
+  // Computes Direct Delta Mesh Skinning (Variant 0) from "Direct Delta Mush
+  // Skinning and Variants"
+  //
+  // Inputs:
+  //   V  #V by 3 list of rest pose vertex positions
+  //   F  #F by 3 list of triangle indices into rows of V
+  //   T  #T list of bone pose transformations
+  //   Omega #V by #T*10 list of precomputated matrix values
+  // Outputs:
+  //   U  #V by 3 list of output vertex positions
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedOmega,
+    typename DerivedU>
+  IGL_INLINE void direct_delta_mush(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const std::vector<
+      Eigen::Affine3d, Eigen::aligned_allocator<Eigen::Affine3d>
+    > & T, /* should eventually be templated more generally than double */
+    const Eigen::MatrixBase<DerivedOmega> & Omega,
+    Eigen::PlainObjectBase<DerivedU> & U);
+
+  // Precomputation
+  //
+  // Inputs:
+  //   V  #V by 3 list of rest pose vertex positions
+  //   F  #F by 3 list of triangle indices into rows of V
+  //   W  #V by #Edges list of weights
+  //   p  number of smoothing iterations
+  //   lambda  rotation smoothing step size
+  //   kappa   translation smoothness step size
+  //   alpha   translation smoothness blending weight
+  // Outputs:
+  //   Omega  #V by #T*10 list of precomputated matrix values
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedW,
+    typename DerivedOmega>
+  IGL_INLINE void direct_delta_mush_precomputation(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    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);
+} // namespace igl
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "direct_delta_mush.cpp"
+#endif
+
+#endif

include/igl/extract_non_manifold_edge_curves.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2016 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "extract_non_manifold_edge_curves.h"
 #include <algorithm>
@@ -17,8 +17,8 @@ typename DerivedF,
 typename DerivedEMAP,
 typename uE2EType >
 IGL_INLINE void igl::extract_non_manifold_edge_curves(
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedEMAP>& /*EMAP*/,
+        const Eigen::MatrixBase<DerivedF>& F,
+        const Eigen::MatrixBase<DerivedEMAP>& /*EMAP*/,
         const std::vector<std::vector<uE2EType> >& uE2E,
         std::vector<std::vector<size_t> >& curves) {
     const size_t num_faces = F.rows();

include/igl/extract_non_manifold_edge_curves.h

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2016 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_NON_MANIFOLD_EDGE_CURVES
 #define IGL_NON_MANIFOLD_EDGE_CURVES
@@ -34,8 +34,8 @@ namespace igl {
         typename DerivedEMAP,
         typename uE2EType>
     IGL_INLINE void extract_non_manifold_edge_curves(
-            const Eigen::PlainObjectBase<DerivedF>& F,
-            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+            const Eigen::MatrixBase<DerivedF>& F,
+            const Eigen::MatrixBase<DerivedEMAP>& EMAP,
             const std::vector<std::vector<uE2EType> >& uE2E,
             std::vector<std::vector<size_t> >& curves);
 }

include/igl/outer_element.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2015 Qingnan Zhou <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "outer_element.h"
 #include <iostream>
@@ -17,13 +17,13 @@ template <
      typename DerivedA
      >
 IGL_INLINE void igl::outer_vertex(
-        const Eigen::PlainObjectBase<DerivedV> & V,
-        const Eigen::PlainObjectBase<DerivedF> & F,
-        const Eigen::PlainObjectBase<DerivedI> & I,
+        const Eigen::MatrixBase<DerivedV> & V,
+        const Eigen::MatrixBase<DerivedF> & F,
+        const Eigen::MatrixBase<DerivedI> & I,
         IndexType & v_index,
         Eigen::PlainObjectBase<DerivedA> & A)
 {
-    // Algorithm: 
+    // Algorithm:
     //    Find an outer vertex (i.e. vertex reachable from infinity)
     //    Return the vertex with the largest X value.
     //    If there is a tie, pick the one with largest Y value.
@@ -38,7 +38,7 @@ IGL_INLINE void igl::outer_vertex(
     typename DerivedV::Scalar outer_val = 0;
     for (size_t i=0; i<num_selected_faces; i++)
     {
-        size_t f = I(i);
+        size_t f = I(i, 0);
         for (size_t j=0; j<3; j++)
         {
             Index v = F(f, j);
@@ -74,7 +74,7 @@ IGL_INLINE void igl::outer_vertex(
     assert(outer_vid != INVALID);
     assert(candidate_faces.size() > 0);
     v_index = outer_vid;
-    A.resize(candidate_faces.size());
+    A.resize(candidate_faces.size(),1);
     std::copy(candidate_faces.begin(), candidate_faces.end(), A.data());
 }
 
@@ -86,9 +86,9 @@ template<
     typename DerivedA
     >
 IGL_INLINE void igl::outer_edge(
-        const Eigen::PlainObjectBase<DerivedV> & V,
-        const Eigen::PlainObjectBase<DerivedF> & F,
-        const Eigen::PlainObjectBase<DerivedI> & I,
+        const Eigen::MatrixBase<DerivedV> & V,
+        const Eigen::MatrixBase<DerivedF> & F,
+        const Eigen::MatrixBase<DerivedI> & I,
         IndexType & v1,
         IndexType & v2,
         Eigen::PlainObjectBase<DerivedA> & A) {
@@ -112,7 +112,7 @@ IGL_INLINE void igl::outer_edge(
     const ScalarArray3& outer_v = V.row(outer_vid);
     assert(candidate_faces.size() > 0);
 
-    auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index 
+    auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index
     {
         if (f[0] == vid) return 0;
         if (f[1] == vid) return 1;
@@ -200,7 +200,7 @@ IGL_INLINE void igl::outer_edge(
 
     v1 = outer_vid;
     v2 = outer_opp_vid;
-    A.resize(incident_faces.size());
+    A.resize(incident_faces.size(),1);
     std::copy(incident_faces.begin(), incident_faces.end(), A.data());
 }
 
@@ -212,10 +212,10 @@ template<
     typename IndexType
     >
 IGL_INLINE void igl::outer_facet(
-        const Eigen::PlainObjectBase<DerivedV> & V,
-        const Eigen::PlainObjectBase<DerivedF> & F,
-        const Eigen::PlainObjectBase<DerivedN> & N,
-        const Eigen::PlainObjectBase<DerivedI> & I,
+        const Eigen::MatrixBase<DerivedV> & V,
+        const Eigen::MatrixBase<DerivedF> & F,
+        const Eigen::MatrixBase<DerivedN> & N,
+        const Eigen::MatrixBase<DerivedI> & I,
         IndexType & f,
         bool & flipped) {
     // Algorithm:
@@ -241,7 +241,7 @@ IGL_INLINE void igl::outer_facet(
     Scalar max_nx = 0;
     size_t outer_fid = INVALID;
     const size_t num_incident_faces = incident_faces.size();
-    for (size_t i=0; i<num_incident_faces; i++) 
+    for (size_t i=0; i<num_incident_faces; i++)
     {
         const auto& fid = incident_faces(i);
         const Scalar nx = N(fid, 0);
@@ -273,8 +273,8 @@ IGL_INLINE void igl::outer_facet(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::outer_facet<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
-template void igl::outer_facet<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, unsigned long&, bool&);
-template void igl::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
-template void igl::outer_facet<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>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::outer_facet<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::outer_facet<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, unsigned long>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, unsigned long&, bool&);
+template void igl::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(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, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::outer_facet<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>, int>(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&, int&, bool&);
 #endif

include/igl/outer_element.h

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2015 Qingan Zhou <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_OUTER_ELEMENT_H
 #define IGL_OUTER_ELEMENT_H
@@ -33,9 +33,9 @@ namespace igl
       typename DerivedA
       >
   IGL_INLINE void outer_vertex(
-          const Eigen::PlainObjectBase<DerivedV> & V,
-          const Eigen::PlainObjectBase<DerivedF> & F,
-          const Eigen::PlainObjectBase<DerivedI> & I,
+          const Eigen::MatrixBase<DerivedV> & V,
+          const Eigen::MatrixBase<DerivedF> & F,
+          const Eigen::MatrixBase<DerivedI> & I,
           IndexType & v_index,
           Eigen::PlainObjectBase<DerivedA> & A);
 
@@ -64,9 +64,9 @@ namespace igl
       typename DerivedA
       >
   IGL_INLINE void outer_edge(
-          const Eigen::PlainObjectBase<DerivedV> & V,
-          const Eigen::PlainObjectBase<DerivedF> & F,
-          const Eigen::PlainObjectBase<DerivedI> & I,
+          const Eigen::MatrixBase<DerivedV> & V,
+          const Eigen::MatrixBase<DerivedF> & F,
+          const Eigen::MatrixBase<DerivedI> & I,
           IndexType & v1,
           IndexType & v2,
           Eigen::PlainObjectBase<DerivedA> & A);
@@ -96,10 +96,10 @@ namespace igl
       typename IndexType
       >
   IGL_INLINE void outer_facet(
-          const Eigen::PlainObjectBase<DerivedV> & V,
-          const Eigen::PlainObjectBase<DerivedF> & F,
-          const Eigen::PlainObjectBase<DerivedN> & N,
-          const Eigen::PlainObjectBase<DerivedI> & I,
+          const Eigen::MatrixBase<DerivedV> & V,
+          const Eigen::MatrixBase<DerivedF> & F,
+          const Eigen::MatrixBase<DerivedN> & N,
+          const Eigen::MatrixBase<DerivedI> & I,
           IndexType & f,
           bool & flipped);
 }

include/igl/project_isometrically_to_plane.cpp

@@ -1,30 +1,30 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2015 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "project_isometrically_to_plane.h"
 #include "edge_lengths.h"
 
 template <
-  typename DerivedV, 
+  typename DerivedV,
   typename DerivedF,
   typename DerivedU,
   typename DerivedUF,
   typename Scalar>
 IGL_INLINE void igl::project_isometrically_to_plane(
-  const Eigen::PlainObjectBase<DerivedV> & V, 
-  const Eigen::PlainObjectBase<DerivedF> & F, 
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
   Eigen::PlainObjectBase<DerivedU> & U,
-  Eigen::PlainObjectBase<DerivedUF> & UF, 
+  Eigen::PlainObjectBase<DerivedUF> & UF,
   Eigen::SparseMatrix<Scalar>& I)
 {
   using namespace std;
   using namespace Eigen;
   assert(F.cols() == 3 && "F should contain triangles");
-  typedef DerivedV MatrixX;
+  typedef Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> MatrixX;
   MatrixX l;
   edge_lengths(V,F,l);
   // Number of faces
@@ -35,9 +35,9 @@ IGL_INLINE void igl::project_isometrically_to_plane(
   // Second corner along x-axis
   U.block(m,0,m,1) = l.col(2);
   // Third corner rotated onto plane
-  U.block(m*2,0,m,1) = 
-    (-l.col(0).array().square() + 
-     l.col(1).array().square() + 
+  U.block(m*2,0,m,1) =
+    (-l.col(0).array().square() +
+     l.col(1).array().square() +
      l.col(2).array().square())/(2.*l.col(2).array());
   U.block(m*2,1,m,1) =
     (l.col(1).array().square()-U.block(m*2,0,m,1).array().square()).sqrt();
@@ -59,5 +59,5 @@ IGL_INLINE void igl::project_isometrically_to_plane(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-template void igl::project_isometrically_to_plane<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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::project_isometrically_to_plane<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/project_isometrically_to_plane.h

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2014 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_PROJECT_ISOMETRICALLY_TO_PLANE_H
 #define IGL_PROJECT_ISOMETRICALLY_TO_PLANE_H
@@ -12,7 +12,8 @@
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
-namespace igl 
+
+namespace igl
 {
   // Project each triangle to the plane
   //
@@ -24,19 +25,19 @@ namespace igl
   // Outputs:
   //   U  #F*3 by 2 list of triangle positions
   //   UF  #F by 3 list of mesh indices into U
-  //   I  #V by #F such that I(i,j) = 1 implies U(j,:) corresponds to V(i,:)
+  //   I  #V by #F*3 such that I(i,j) = 1 implies U(j,:) corresponds to V(i,:)
   //
   template <
-    typename DerivedV, 
+    typename DerivedV,
     typename DerivedF,
     typename DerivedU,
     typename DerivedUF,
     typename Scalar>
   IGL_INLINE void project_isometrically_to_plane(
-    const Eigen::PlainObjectBase<DerivedV> & V, 
-    const Eigen::PlainObjectBase<DerivedF> & F, 
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
     Eigen::PlainObjectBase<DerivedU> & U,
-    Eigen::PlainObjectBase<DerivedUF> & UF, 
+    Eigen::PlainObjectBase<DerivedUF> & UF,
     Eigen::SparseMatrix<Scalar>& I);
 }
 

tests/include/igl/direct_delta_mush.cpp

@@ -0,0 +1,25 @@
+#include <test_common.h>
+#include <igl/direct_delta_mush.h>
+#include <igl/adjacency_list.h>
+#include <iostream>
+
+TEST_CASE("direct_delta_mush: cube", "[igl]")
+{
+  Eigen::MatrixXd V, W, U, Omega;
+  Eigen::MatrixXi F;
+
+  // Test that direct delta mush with identity transform reproduces the rest state of the mesh.
+  // For simplicity, testing on a cube of dimensions 1.0 x 1.0 x 1.0,
+  // with all vertices bound strictly to one and only one imaginary bone (weight is a column of 1s)
+  igl::read_triangle_mesh(test_common::data_path("cube.off"), V, F);
+  Eigen::SparseMatrix<double> W_sparse = Eigen::MatrixXd::Ones(V.rows(), 1).sparseView();
+
+  // Parameters such as p, lambda, kappa, alpha do not matter given identity transform
+  igl::direct_delta_mush_precomputation(V, F, W_sparse, 1, 1., 0.5, 0.5, Omega);
+
+  std::vector<Eigen::Affine3d, Eigen::aligned_allocator<Eigen::Affine3d>> T_list;
+  T_list.push_back(Eigen::Affine3d::Identity());
+  igl::direct_delta_mush(V, F, T_list, Omega, U);
+
+  test_common::assert_near(U, V, 1e-4);
+}

tests/test_common.h

@@ -16,6 +16,7 @@
 #include <algorithm>
 #include <tuple>
 
+
 // Disable lengthy tests in debug mode
 #ifdef NDEBUG
 #define IGL_DEBUG_OFF ""

tutorial/408_DirectDeltaMush/CMakeLists.txt

@@ -0,0 +1,6 @@
+get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+project(${PROJECT_NAME})
+
+add_executable(${PROJECT_NAME}_bin main.cpp)
+target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw tutorials)
+

tutorial/408_DirectDeltaMush/main.cpp

@@ -0,0 +1,125 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/readTGF.h>
+#include <igl/readDMAT.h>
+#include <igl/lbs_matrix.h>
+#include <igl/deform_skeleton.h>
+#include <igl/direct_delta_mush.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include "tutorial_shared_path.h"
+#include <Eigen/Geometry>
+#include <vector>
+
+int main(int argc, char * argv[]) 
+{
+
+  Eigen::MatrixXd V,U,C,W,T,M,Omega;
+  Eigen::MatrixXi F,BE;
+  igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/elephant.obj",V,F);
+  igl::readTGF(           TUTORIAL_SHARED_PATH "/elephant.tgf",C,BE);
+  igl::readDMAT(          TUTORIAL_SHARED_PATH "/elephant-weights.dmat",W);
+  igl::readDMAT(          TUTORIAL_SHARED_PATH "/elephant-anim.dmat",T);
+  // convert weight to piecewise-rigid weights to stress test DDM
+  for (int i = 0; i < W.rows(); ++i)
+  {
+    int maxj;
+    W.row(i).maxCoeff(&maxj);
+    for (int j = 0; j < W.cols(); j++)
+    {
+      W(i, j) = double(maxj == j);
+    }
+  }
+
+  igl::lbs_matrix(V,W,M);
+
+  int p = 20;
+  float lambda = 3; // 0 < lambda
+  float kappa = 1; // 0 < kappa < lambda
+  float alpha = 0.8; // 0 <= alpha < 1
+  Eigen::SparseMatrix<double> Wsparse =  W.sparseView();
+  igl::direct_delta_mush_precomputation(V, F,Wsparse, p, lambda, kappa, alpha, Omega);
+
+  igl::opengl::glfw::Viewer viewer;
+  int frame = 0;
+  const int pr_id = viewer.selected_data_index;
+  viewer.append_mesh();
+  const int ddm_id = viewer.selected_data_index;
+  Eigen::RowVector3d offset(1.1*(V.col(0).maxCoeff()-V.col(0).minCoeff()),0,0);
+
+  viewer.callback_pre_draw = [&](igl::opengl::glfw::Viewer &) -> bool
+  {
+    if(viewer.core().is_animating)
+    {
+      const Eigen::Map<Eigen::MatrixXd> Tf(
+        T.data()+frame*T.rows(),4*W.cols(),3);
+      U = (M*Tf).rowwise()-offset;
+
+      Eigen::MatrixXd Cf;
+      Eigen::MatrixXi BEf;
+      igl::deform_skeleton(C,BE,Tf,Cf,BEf);
+      viewer.data(pr_id).set_edges(Cf,BEf, Eigen::RowVector3d(1,1,1));
+
+      viewer.data(pr_id).set_vertices(U);
+      viewer.data(pr_id).compute_normals();
+
+      {
+        std::vector<Eigen::Affine3d, Eigen::aligned_allocator<Eigen::Affine3d>> 
+          T_list(BE.rows());
+        for (int e = 0; e < BE.rows(); e++)
+        {
+          T_list[e] = Eigen::Affine3d::Identity();
+          T_list[e].matrix().block(0, 0, 3, 4) = Tf.block(e*4,0,4,3).transpose();
+        }
+        igl::direct_delta_mush(V, F, T_list, Omega, U);
+      }
+      U.rowwise() += offset;
+      viewer.data(ddm_id).set_vertices(U);
+      viewer.data(ddm_id).compute_normals();
+
+      frame++;
+      if(frame == T.cols())
+      {
+        frame = 0;
+        viewer.core().is_animating = false;
+      }
+    }
+    return false;
+  };
+  viewer.callback_key_pressed = [&](igl::opengl::glfw::Viewer &, unsigned int key, int mod)
+  {
+    switch(key)
+    {
+    case ' ':
+      viewer.core().is_animating = !viewer.core().is_animating;
+      return true;
+    }
+    return false;
+  };
+
+
+
+  for(auto & id : {pr_id,ddm_id})
+  {
+    if(id == pr_id)
+    {
+      viewer.data(id).set_mesh( (V.rowwise()-offset*1.0).eval() ,F);
+      viewer.data(id).set_colors(Eigen::RowVector3d(214./255.,170./255.,148./255.));
+      viewer.data(id).set_edges(C,BE, Eigen::RowVector3d(1,1,1));
+    }else if(id == ddm_id){
+      viewer.data(id).set_mesh( (V.rowwise()+offset*1.0).eval() ,F);
+      viewer.data(id).set_colors(Eigen::RowVector3d(132./255.,214./255.,105./255.));
+    }
+    viewer.data(id).show_lines = false;
+    viewer.data(id).set_face_based(true);
+    viewer.data(id).show_overlay_depth = false;
+  }
+  viewer.core().is_animating = false;
+  viewer.core().animation_max_fps = 24.;
+  //viewer.data().set_colors(V,F);
+
+
+  viewer.launch_init();
+  viewer.core().align_camera_center(V);
+
+  viewer.launch_rendering(true);
+  viewer.launch_shut();
+}

tutorial/CMakeLists.txt

@@ -92,6 +92,7 @@ if(TUTORIALS_CHAPTER4)
   add_subdirectory("405_AsRigidAsPossible")
   add_subdirectory("406_FastAutomaticSkinningTransformations")
   add_subdirectory("407_BiharmonicCoordinates")
+  add_subdirectory("408_DirectDeltaMush")
 endif()
 
 # Chapter 5