biharmonic coordinates, impl of 'linear subspace design...', templates, overloads

Former-commit-id: 9e0a25141c4cb8575e8617293060da2665efee52

include/igl/biharmonic_coordinates.cpp

@@ -0,0 +1,117 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "biharmonic_coordinates.h"
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include "min_quad_with_fixed.h"
+#include "normal_derivative.h"
+#include "on_boundary.h"
+#include <Eigen/Sparse>
+
+template <
+  typename DerivedV,
+  typename DerivedT,
+  typename SType,
+  typename DerivedW>
+IGL_INLINE bool igl::biharmonic_coordinates(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedT> & T,
+  const std::vector<std::vector<SType> > & S,
+  Eigen::PlainObjectBase<DerivedW> & W)
+{
+  using namespace Eigen;
+  using namespace std;
+  // This is not the most efficient way to build A, but follows "Linear
+  // Subspace Design for Real-Time Shape Deformation" [Wang et al. 2015]. 
+  SparseMatrix<double> A;
+  {
+    SparseMatrix<double> N,Z,L,K,M;
+    normal_derivative(V,T,N);
+    Array<bool,Dynamic,1> I;
+    Array<bool,Dynamic,Dynamic> C;
+    on_boundary(T,I,C);
+    {
+      std::vector<Triplet<double> >ZIJV;
+      for(int t =0;t<T.rows();t++)
+      {
+        for(int f =0;f<T.cols();f++)
+        {
+          if(C(t,f))
+          {
+            const int i = t+f*T.rows();
+            for(int c = 1;c<T.cols();c++)
+            {
+              ZIJV.emplace_back(T(t,(f+c)%T.cols()),i,1);
+            }
+          }
+        }
+      }
+      Z.resize(V.rows(),N.rows());
+      Z.setFromTriplets(ZIJV.begin(),ZIJV.end());
+      N = (Z*N).eval();
+    }
+    cotmatrix(V,T,L);
+    K = N+L;
+    massmatrix(V,T,MASSMATRIX_TYPE_DEFAULT,M);
+    DiagonalMatrix<double,Dynamic> Minv = 
+      ((VectorXd)M.diagonal().array().inverse()).asDiagonal();
+    A = K.transpose() * (Minv * K);
+  }
+  // Vertices in point handles
+  const size_t mp = 
+    count_if(S.begin(),S.end(),[](const vector<int> & h){return h.size()==1;});
+  // number of region handles
+  const size_t r = S.size()-mp;
+  // Vertices in region handles
+  size_t mr = 0;
+  for(const auto & h : S)
+  {
+    if(h.size() > 1)
+    {
+      mr += h.size();
+    }
+  }
+  const size_t dim = T.cols()-1;
+  // Might as well be dense... I think...
+  MatrixXd J = MatrixXd::Zero(mp+mr,mp+r*(dim+1));
+  VectorXi b(mp+mr);
+  MatrixXd H(mp+r*(dim+1),dim);
+  {
+    int v = 0;
+    int c = 0;
+    for(int h = 0;h<S.size();h++)
+    {
+      if(S[h].size()==1)
+      {
+        H.row(c) = V.block(S[h][0],0,1,dim);
+        J(v,c++) = 1;
+        b(v) = S[h][0];
+        v++;
+      }else
+      {
+        assert(S[h].size() >= dim+1);
+        for(int p = 0;p<S[h].size();p++)
+        {
+          for(int d = 0;d<dim;d++)
+          {
+            J(v,c+d) = V(S[h][p],d);
+          }
+          J(v,c+dim) = 1;
+          b(v) = S[h][p];
+          v++;
+        }
+        H.block(c,0,dim+1,dim).setIdentity();
+        c+=dim+1;
+      }
+    }
+  }
+  // minimize    ½ W' A W' 
+  // subject to  W(b,:) = J
+  return min_quad_with_fixed(
+    A,VectorXd::Zero(A.rows()).eval(),b,J,{},VectorXd(),true,W);
+}

include/igl/biharmonic_coordinates.h

@@ -0,0 +1,77 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_BIHARMONIC_COORDINATES_H
+#define IGL_BIHARMONIC_COORDINATES_H
+#include "igl_inline.h"
+#include <Eigen/Dense>
+#include <vector>
+namespace igl
+{
+  // Compute "discrete biharmonic generalized barycentric coordinates" as
+  // described in "Linear Subspace Design for Real-Time Shape Deformation"
+  // [Wang et al. 2015]. Not to be confused with "Bounded Biharmonic Weights
+  // for Real-Time Deformation" [Jacobson et al. 2011] or "Biharmonic
+  // Coordinates" (2D complex barycentric coordinates) [Weber et al. 2012].
+  // These weights minimize a discrete version of the squared Laplacian energy
+  // subject to positional interpolation constraints at selected vertices
+  // (point handles) and transformation interpolation constraints at regions
+  // (region handles).
+  //
+  // Templates:
+  //   HType  should be a simple index type e.g. `int`,`size_t`
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   T  #T by dim+1 list of / triangle indices into V      if dim=2
+  //                          \ tetrahedron indices into V   if dim=3
+  //   S  #point-handles+#region-handles list of lists of selected vertices for
+  //     each handle. Point handles should have singleton lists and region
+  //     handles should have lists of size at least dim+1 (and these points
+  //     should be in general position).
+  // Outputs:
+  //   W  #V by #points-handles+(#region-handles * dim+1) matrix of weights so
+  //     that columns correspond to each handles generalized barycentric
+  //     coordinates (for point-handles) or animation space weights (for region
+  //     handles).
+  // returns true only on success
+  //
+  // Example:
+  //
+  //     MatrixXd W;
+  //     igl::biharmonic_coordinates(V,F,S,W);
+  //     const size_t dim = T.cols()-1;
+  //     MatrixXd H(W.cols(),dim);
+  //     {
+  //       int c = 0;
+  //       for(int h = 0;h<S.size();h++)
+  //       {
+  //         if(S[h].size()==1)
+  //         {
+  //           H.row(c++) = V.block(S[h][0],0,1,dim);
+  //         }else
+  //         {
+  //           H.block(c,0,dim+1,dim).setIdentity();
+  //           c+=dim+1;
+  //         }
+  //       }
+  //     }
+  //     assert( (V-(W*H)).array().maxCoeff() < 1e-7 );
+  template <
+    typename DerivedV,
+    typename DerivedT,
+    typename SType,
+    typename DerivedW>
+  IGL_INLINE bool biharmonic_coordinates(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedT> & T,
+    const std::vector<std::vector<SType> > & S,
+    Eigen::PlainObjectBase<DerivedW> & W);
+};
+#  ifndef IGL_STATIC_LIBRARY
+#    include "biharmonic_coordinates.cpp"
+#  endif
+#endif

include/igl/cotmatrix.cpp

@@ -14,7 +14,6 @@
 
 // Bug in unsupported/Eigen/SparseExtra needs iostream first
 #include <iostream>
-#include <unsupported/Eigen/SparseExtra>
 
 template <typename DerivedV, typename DerivedF, typename Scalar>
 IGL_INLINE void igl::cotmatrix(

include/igl/face_occurrences.cpp

@@ -36,7 +36,7 @@ IGL_INLINE void igl::face_occurrences(
     {
       // increment count
       counts[sortedF[i]]++;
-      assert(counts[sortedF[i]] == 2);
+      assert(counts[sortedF[i]] == 2 && "Input should be manifold");
     }
   }
 

include/igl/min_quad_with_fixed.cpp

@@ -64,6 +64,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   assert((kr == 0 || known.maxCoeff() < n) && "known indices should be in [0,n)");
   assert(neq <= n && "Number of equality constraints should be less than DOFs");
 
+
   // cache known
   data.known = known;
   // get list of unknown indices
@@ -539,6 +540,31 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
   return min_quad_with_fixed_solve(data,B,Y,Beq,Z,sol);
 }
 
+template <
+  typename T,
+  typename Derivedknown,
+  typename DerivedB,
+  typename DerivedY,
+  typename DerivedBeq,
+  typename DerivedZ>
+IGL_INLINE bool igl::min_quad_with_fixed(
+  const Eigen::SparseMatrix<T>& A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
+  const Eigen::PlainObjectBase<Derivedknown> & known,
+  const Eigen::PlainObjectBase<DerivedY> & Y,
+  const Eigen::SparseMatrix<T>& Aeq,
+  const Eigen::PlainObjectBase<DerivedBeq> & Beq,
+  const bool pd,
+  Eigen::PlainObjectBase<DerivedZ> & Z)
+{
+  min_quad_with_fixed_data<T> data;
+  if(!min_quad_with_fixed_precompute(A,known,Aeq,pd,data))
+  {
+    return false;
+  }
+  return min_quad_with_fixed_solve(data,B,Y,Beq,Z);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);

include/igl/min_quad_with_fixed.h

@@ -100,7 +100,22 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedY> & Y,
     const Eigen::PlainObjectBase<DerivedBeq> & Beq,
     Eigen::PlainObjectBase<DerivedZ> & Z);
-
+  template <
+    typename T,
+    typename Derivedknown,
+    typename DerivedB,
+    typename DerivedY,
+    typename DerivedBeq,
+    typename DerivedZ>
+  IGL_INLINE bool min_quad_with_fixed(
+    const Eigen::SparseMatrix<T>& A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
+    const Eigen::PlainObjectBase<Derivedknown> & known,
+    const Eigen::PlainObjectBase<DerivedY> & Y,
+    const Eigen::SparseMatrix<T>& Aeq,
+    const Eigen::PlainObjectBase<DerivedBeq> & Beq,
+    const bool pd,
+    Eigen::PlainObjectBase<DerivedZ> & Z);
 }
 
 template <typename T>

include/igl/normal_derivative.cpp

@@ -0,0 +1,113 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "normal_derivative.h"
+#include "cotmatrix_entries.h"
+#include "slice.h"
+#include <cassert>
+
+template <
+  typename DerivedV, 
+  typename DerivedEle, 
+  typename Scalar>
+IGL_INLINE void igl::normal_derivative(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedEle> & Ele,
+  Eigen::SparseMatrix<Scalar>& DD)
+{
+  using namespace Eigen;
+  using namespace std;
+  // Element simplex-size 
+  const size_t ss = Ele.cols();
+  assert( ((ss==3) || (ss==4)) && "Only triangles or tets");
+  // cotangents
+  Matrix<Scalar,Dynamic,Dynamic> C;
+  cotmatrix_entries(V,Ele,C);
+  vector<Triplet<Scalar> > IJV;
+  // Number of elements
+  const size_t m = Ele.rows();
+  // Number of vertices
+  const size_t n = V.rows();
+  switch(ss)
+  {
+    default:
+      assert(false);
+      return;
+    case 4:
+    {
+      const MatrixXi DDJ = 
+        slice(
+          Ele,
+          (VectorXi(24)<<
+            1,0,2,0,3,0,2,1,3,1,0,1,3,2,0,2,1,2,0,3,1,3,2,3).finished(),
+          2);
+      MatrixXi DDI(m,24);
+      for(size_t f = 0;f<4;f++)
+      {
+        const auto & I = (VectorXi::LinSpaced(m,0,m-1).array()+f*m).eval();
+        for(size_t r = 0;r<6;r++)
+        {
+          DDI.col(f*6+r) = I;
+        }
+      }
+      const DiagonalMatrix<Scalar,24,24> S = 
+        (Matrix<Scalar,2,1>(1,-1).template replicate<12,1>()).asDiagonal();
+      Matrix<Scalar,Dynamic,Dynamic> DDV = 
+        slice(
+          C,
+          (VectorXi(24)<<
+            2,2,1,1,3,3,0,0,4,4,2,2,5,5,1,1,0,0,3,3,4,4,5,5).finished(),
+          2);
+      DDV *= S;
+
+      IJV.resize(DDV.size());
+      for(size_t f = 0;f<6*4;f++)
+      {
+        for(size_t e = 0;e<m;e++)
+        {
+          IJV.push_back(Triplet<Scalar>(DDI(e,f),DDJ(e,f),DDV(e,f)));
+        }
+      }
+      DD.resize(m*4,n);
+      DD.setFromTriplets(IJV.begin(),IJV.end());
+      break;
+    }
+    case 3:
+    {
+      const MatrixXi DDJ = 
+        slice(Ele,(VectorXi(12)<<2,0,1,0,0,1,2,1,1,2,0,2).finished(),2);
+      MatrixXi DDI(m,12);
+      for(size_t f = 0;f<3;f++)
+      {
+        const auto & I = (VectorXi::LinSpaced(m,0,m-1).array()+f*m).eval();
+        for(size_t r = 0;r<4;r++)
+        {
+          DDI.col(f*4+r) = I;
+        }
+      }
+      const DiagonalMatrix<Scalar,12,12> S = 
+        (Matrix<Scalar,2,1>(1,-1).template replicate<6,1>()).asDiagonal();
+      Matrix<Scalar,Dynamic,Dynamic> DDV = 
+        slice(C,(VectorXi(12)<<1,1,2,2,2,2,0,0,0,0,1,1).finished(),2);
+      DDV *= S;
+
+      IJV.resize(DDV.size());
+      for(size_t f = 0;f<12;f++)
+      {
+        for(size_t e = 0;e<m;e++)
+        {
+          IJV.push_back(Triplet<Scalar>(DDI(e,f),DDJ(e,f),DDV(e,f)));
+        }
+      }
+      DD.resize(m*3,n);
+      DD.setFromTriplets(IJV.begin(),IJV.end());
+      break;
+    }
+  }
+
+}
+

include/igl/normal_derivative.h

@@ -0,0 +1,44 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_NORMAL_DERIVATIVE_H
+#define IGL_NORMAL_DERIVATIVE_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl 
+{
+  // NORMAL_DERIVATIVE Computes the directional derivative **normal** to
+  // **all** (half-)edges of a triangle mesh (not just boundary edges). These
+  // are integrated along the edge: they're the per-face constant gradient dot
+  // the rotated edge vector (unit rotated edge vector for direction then
+  // magnitude for integration).
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by 3|4 list of triangle|tetrahedron indices into V
+  // Outputs:
+  //   DD  #F*3|4 by #V sparse matrix representing operator to compute
+  //     directional derivative with respect to each facet of each element.
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedEle, 
+    typename Scalar>
+  IGL_INLINE void normal_derivative(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedEle> & Ele,
+    Eigen::SparseMatrix<Scalar>& DD);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "normal_derivative.cpp"
+#endif
+
+#endif
+

include/igl/on_boundary.cpp

@@ -20,43 +20,89 @@ IGL_INLINE void igl::on_boundary(
   std::vector<std::vector<bool> > & C)
 {
   using namespace std;
-
-  // Get a list of all faces
-  vector<vector<IntegerT> > F(T.size()*4,vector<IntegerT>(3));
-  // Gather faces, loop over tets
-  for(int i = 0; i< (int)T.size();i++)
+  if(T.empty())
   {
-    assert(T[i].size() == 4);
-    // get face in correct order
-    F[i*4+0][0] = T[i][1];
-    F[i*4+0][1] = T[i][3];
-    F[i*4+0][2] = T[i][2];
-    // get face in correct order
-    F[i*4+1][0] = T[i][0];
-    F[i*4+1][1] = T[i][2];
-    F[i*4+1][2] = T[i][3];
-    // get face in correct order
-    F[i*4+2][0] = T[i][0];
-    F[i*4+2][1] = T[i][3];
-    F[i*4+2][2] = T[i][1];
-    // get face in correct order
-    F[i*4+3][0] = T[i][0];
-    F[i*4+3][1] = T[i][1];
-    F[i*4+3][2] = T[i][2];
+    I.clear();
+    C.clear();
+    return;
   }
-  // Counts
-  vector<int> FC;
-  face_occurrences(F,FC);
-  C.resize(T.size(),vector<bool>(4));
-  I.resize(T.size(),false);
-  for(int i = 0; i< (int)T.size();i++)
+
+  switch(T[0].size())
   {
-    for(int j = 0;j<4;j++)
+    case 3:
+    {
+      // Get a list of all faces
+      vector<vector<IntegerT> > F(T.size()*3,vector<IntegerT>(2));
+      // Gather faces, loop over tets
+      for(int i = 0; i< (int)T.size();i++)
+      {
+        assert(T[i].size() == 3);
+        // get face in correct order
+        F[i*3+0][0] = T[i][1];
+        F[i*3+0][1] = T[i][2];
+        F[i*3+1][0] = T[i][2];
+        F[i*3+1][1] = T[i][0];
+        F[i*3+2][0] = T[i][0];
+        F[i*3+2][1] = T[i][1];
+      }
+      // Counts
+      vector<int> FC;
+      face_occurrences(F,FC);
+      C.resize(T.size(),vector<bool>(3));
+      I.resize(T.size(),false);
+      for(int i = 0; i< (int)T.size();i++)
+      {
+        for(int j = 0;j<3;j++)
+        {
+          assert(FC[i*3+j] == 2 || FC[i*3+j] == 1);
+          C[i][j] = FC[i*3+j]==1;
+          // if any are on boundary set to true
+          I[i] = I[i] || C[i][j];
+        }
+      }
+      return;
+    }
+    case 4:
     {
-      assert(FC[i*4+j] == 2 || FC[i*4+j] == 1);
-      C[i][j] = FC[i*4+j]==1;
-      // if any are on boundary set to true
-      I[i] = I[i] || C[i][j];
+      // Get a list of all faces
+      vector<vector<IntegerT> > F(T.size()*4,vector<IntegerT>(3));
+      // Gather faces, loop over tets
+      for(int i = 0; i< (int)T.size();i++)
+      {
+        assert(T[i].size() == 4);
+        // get face in correct order
+        F[i*4+0][0] = T[i][1];
+        F[i*4+0][1] = T[i][3];
+        F[i*4+0][2] = T[i][2];
+        // get face in correct order
+        F[i*4+1][0] = T[i][0];
+        F[i*4+1][1] = T[i][2];
+        F[i*4+1][2] = T[i][3];
+        // get face in correct order
+        F[i*4+2][0] = T[i][0];
+        F[i*4+2][1] = T[i][3];
+        F[i*4+2][2] = T[i][1];
+        // get face in correct order
+        F[i*4+3][0] = T[i][0];
+        F[i*4+3][1] = T[i][1];
+        F[i*4+3][2] = T[i][2];
+      }
+      // Counts
+      vector<int> FC;
+      face_occurrences(F,FC);
+      C.resize(T.size(),vector<bool>(4));
+      I.resize(T.size(),false);
+      for(int i = 0; i< (int)T.size();i++)
+      {
+        for(int j = 0;j<4;j++)
+        {
+          assert(FC[i*4+j] == 2 || FC[i*4+j] == 1);
+          C[i][j] = FC[i*4+j]==1;
+          // if any are on boundary set to true
+          I[i] = I[i] || C[i][j];
+        }
+      }
+      return;
     }
   }
 
@@ -73,7 +119,7 @@ IGL_INLINE void igl::on_boundary(
   Eigen::PlainObjectBase<DerivedI>& I,
   Eigen::PlainObjectBase<DerivedC>& C)
 {
-  assert(T.cols() == 0 || T.cols() == 4);
+  assert(T.cols() == 0 || T.cols() == 4 || T.cols() == 3);
   using namespace std;
   using namespace Eigen;
   // Cop out: use vector of vectors version

include/igl/on_boundary.h

@@ -17,16 +17,17 @@
 
 namespace igl
 {
-  // BOUNDARY_FACES Determine boundary faces of tetrahedra stored in T
+  // ON_BOUNDARY Determine boundary facets of mesh elements stored in T
   //
   // Templates:
   //   IntegerT  integer-value: i.e. int
   //   IntegerF  integer-value: i.e. int
   // Input:
-  //  T  tetrahedron index list, m by 4, where m is the number of tetrahedra
+  //  T  triangle|tetrahedron index list, m by 3|4, where m is the number of
+  //    elements
   // Output:
   //  I  m long list of bools whether tet is on boundary
-  //  C  m by 4 list of bools whether opposite face is on boundary
+  //  C  m by 3|4 list of bools whether opposite facet is on boundary
   //
   template <typename IntegerT>
   IGL_INLINE void on_boundary(

include/igl/slice.cpp

@@ -10,10 +10,6 @@
 
 #include <vector>
 
-// Bug in unsupported/Eigen/SparseExtra needs iostream first
-#include <iostream>
-#include <unsupported/Eigen/SparseExtra>
-
 template <typename T>
 IGL_INLINE void igl::slice(
   const Eigen::SparseMatrix<T>& X,

include/igl/slice_mask.h

@@ -10,6 +10,7 @@
 #include "igl_inline.h"
 
 #include <Eigen/Sparse>
+#include <Eigen/Core>
 namespace igl
 {
   // Act like the matlab X(row_mask,col_mask) operator, where