arap works as parameterization (note interface changes)

Former-commit-id: 6a254e9a410c4446ef06483ecb0d81e4f810ed17

examples/arap/Makefile

@@ -13,7 +13,7 @@ LIBIGL_INC=-I$(LIBIGL)/include
 LIBIGL_LIB=-L$(LIBIGL)/lib -ligl -liglpng -liglsvd3x3
 
 EIGEN3_INC=-I/opt/local/include/eigen3 -I/opt/local/include/eigen3/unsupported
-CFLAGS+=-DIGL_HEADER_ONLY
+#CFLAGS+=-DIGL_HEADER_ONLY
 
 
 #EMBREE=$(LIBIGL)/external/embree

examples/arap/example.cpp

@@ -218,7 +218,7 @@ bool init_arap()
   //arap_data.h = 0.5;
   //arap_data.max_iter = 100;
   //partition(W,100,arap_data.G,_S,_D);
-  return arap_precomputation(V,*Ele,b,arap_data);
+  return arap_precomputation(V,*Ele,V.cols(),b,arap_data);
 }
 
 bool update_arap()

include/igl/arap_rhs.cpp

@@ -14,13 +14,14 @@
 IGL_INLINE void igl::arap_rhs(
   const Eigen::MatrixXd & V, 
   const Eigen::MatrixXi & F,
+  const int dim,
   const igl::ARAPEnergyType energy,
   Eigen::SparseMatrix<double>& K)
 {
   using namespace igl;
   using namespace Eigen;
-  // Number of dimensions
-  int dim = V.cols();
+  //// Number of dimensions
+  //int dim = V.cols();
   //// Number of mesh vertices
   //int n = V.rows();
   //// Number of mesh elements

include/igl/arap_rhs.h

@@ -18,8 +18,10 @@ namespace igl
   // ARAP_RHS build right-hand side constructor of global poisson solve for
   // various Arap energies
   // Inputs:
-  //   V  #V by dim list of initial domain positions
+  //   V  #V by Vdim list of initial domain positions
   //   F  #F by 3 list of triangle indices into V
+  //   dim  dimension being used at solve time. For deformation usually dim =
+  //     V.cols(), for surface parameterization V.cols() = 3 and dim = 2
   //   energy  igl::ARAPEnergyType enum value defining which energy is being
   //     used. See igl::ARAPEnergyType.h for valid options and explanations.
   // Outputs:
@@ -30,6 +32,7 @@ namespace igl
   IGL_INLINE void arap_rhs(
     const Eigen::MatrixXd & V, 
     const Eigen::MatrixXi & F,
+    const int dim,
     const igl::ARAPEnergyType energy,
     Eigen::SparseMatrix<double>& K);
 }

include/igl/covariance_scatter_matrix.cpp

@@ -17,6 +17,7 @@
 IGL_INLINE void igl::covariance_scatter_matrix(
   const Eigen::MatrixXd & V, 
   const Eigen::MatrixXi & F,
+  const int dim,
   const ARAPEnergyType energy,
   Eigen::SparseMatrix<double>& CSM)
 {
@@ -25,8 +26,8 @@ IGL_INLINE void igl::covariance_scatter_matrix(
   // number of mesh vertices
   int n = V.rows();
   assert(n > F.maxCoeff());
-  // dimension of mesh
-  int dim = V.cols();
+  //// dimension of mesh
+  //int dim = V.cols();
   // Number of mesh elements
   int m = F.rows();
 

include/igl/covariance_scatter_matrix.h

@@ -19,6 +19,8 @@ namespace igl
   // Inputs:
   //   V  #V by dim list of initial domain positions
   //   F  #F by 3 list of triangle indices into V
+  //   dim  dimension being used at solve time. For deformation usually dim =
+  //     V.cols(), for surface parameterization V.cols() = 3 and dim = 2
   //   energy  ARAPEnergyType enum value defining which energy is being used.
   //     See ARAPEnergyType.h for valid options and explanations.
   // Outputs:
@@ -28,6 +30,7 @@ namespace igl
   IGL_INLINE void covariance_scatter_matrix(
     const Eigen::MatrixXd & V, 
     const Eigen::MatrixXi & F,
+    const int dim,
     const ARAPEnergyType energy,
     Eigen::SparseMatrix<double>& CSM);
 }

include/igl/svd3x3/arap.cpp

@@ -29,6 +29,7 @@ template <
 IGL_INLINE bool igl::arap_precomputation(
   const Eigen::PlainObjectBase<DerivedV> & V,
   const Eigen::PlainObjectBase<DerivedF> & F,
+  const int dim,
   const Eigen::PlainObjectBase<Derivedb> & b,
   ARAPData & data)
 {
@@ -45,10 +46,12 @@ IGL_INLINE bool igl::arap_precomputation(
   data.b = b;
   //assert(F.cols() == 3 && "For now only triangles");
   // dimension
-  const int dim = V.cols();
+  //const int dim = V.cols();
+  assert((dim == 3 || dim ==2) && "dim should be 2 or 3");
+  data.dim = dim;
   //assert(dim == 3 && "Only 3d supported");
   // Defaults
-  data.f_ext = Eigen::MatrixXd::Zero(n,dim);
+  data.f_ext = Eigen::MatrixXd::Zero(n,data.dim);
 
   SparseMatrix<Scalar> L;
   cotmatrix(V,F,L);
@@ -59,7 +62,7 @@ IGL_INLINE bool igl::arap_precomputation(
     switch(F.cols())
     {
       case 3:
-        if(dim == 3)
+        if(data.dim == 3)
         {
           eff_energy = ARAP_ENERGY_TYPE_SPOKES_AND_RIMS;
         }else
@@ -78,7 +81,7 @@ IGL_INLINE bool igl::arap_precomputation(
 
   // Get covariance scatter matrix, when applied collects the covariance
   // matrices used to fit rotations to during optimization
-  covariance_scatter_matrix(V,F,eff_energy,data.CSM);
+  covariance_scatter_matrix(V,F,data.dim,eff_energy,data.CSM);
 
   // Get group sum scatter matrix, when applied sums all entries of the same
   // group according to G
@@ -112,12 +115,12 @@ IGL_INLINE bool igl::arap_precomputation(
     group_sum_matrix(data.G,G_sum);
   }
   SparseMatrix<double> G_sum_dim;
-  repdiag(G_sum,dim,G_sum_dim);
+  repdiag(G_sum,data.dim,G_sum_dim);
   assert(G_sum_dim.cols() == data.CSM.rows());
   data.CSM = (G_sum_dim * data.CSM).eval();
 
 
-  arap_rhs(V,F,eff_energy,data.K);
+  arap_rhs(V,F,data.dim,eff_energy,data.K);
 
   SparseMatrix<double> Q = (-0.5*L).eval();
 
@@ -130,8 +133,8 @@ IGL_INLINE bool igl::arap_precomputation(
     SparseMatrix<double> DQ = 0.5/(h*h)*data.M;
     Q += DQ;
     // Dummy external forces
-    data.f_ext = MatrixXd::Zero(n,dim);
-    data.vel = MatrixXd::Zero(n,dim);
+    data.f_ext = MatrixXd::Zero(n,data.dim);
+    data.vel = MatrixXd::Zero(n,data.dim);
   }
 
   return min_quad_with_fixed_precompute(
@@ -150,7 +153,10 @@ IGL_INLINE bool igl::arap_solve(
   using namespace Eigen;
   using namespace std;
   assert(data.b.size() == bc.rows());
-  const int dim = bc.cols();
+  if(bc.size() > 0)
+  {
+    assert(bc.cols() == data.dim && "bc.cols() match data.dim");
+  }
   const int n = data.n;
   int iter = 0;
   if(U.size() == 0)
@@ -158,8 +164,11 @@ IGL_INLINE bool igl::arap_solve(
     // terrible initial guess.. should at least copy input mesh
 #ifndef NDEBUG
     cerr<<"arap_solve: Using terrible initial guess for U. Try U = V."<<endl;
-    U = MatrixXd::Zero(data.n,dim);
 #endif
+    U = MatrixXd::Zero(data.n,data.dim);
+  }else
+  {
+    assert(U.cols() == data.dim && "U.cols() match data.dim");
   }
   // changes each arap iteration
   MatrixXd U_prev = U;
@@ -178,11 +187,11 @@ IGL_INLINE bool igl::arap_solve(
       U.row(data.b(bi)) = bc.row(bi);
     }
 
-    const auto & Udim = U.replicate(dim,1);
+    const auto & Udim = U.replicate(data.dim,1);
     MatrixXd S = data.CSM * Udim;
 
-    MatrixXd R(dim,data.CSM.rows());
-    if(dim == 2)
+    MatrixXd R(data.dim,data.CSM.rows());
+    if(data.dim == 2)
     {
       fit_rotations_planar(S,R);
     }else
@@ -200,7 +209,7 @@ IGL_INLINE bool igl::arap_solve(
 
 
     // Number of rotations: #vertices or #elements
-    int num_rots = data.K.cols()/dim/dim;
+    int num_rots = data.K.cols()/data.dim/data.dim;
     // distribute group rotations to vertices in each group
     MatrixXd eff_R;
     if(data.G.size() == 0)
@@ -209,11 +218,11 @@ IGL_INLINE bool igl::arap_solve(
       eff_R = R;
     }else
     {
-      eff_R.resize(dim,num_rots*dim);
+      eff_R.resize(data.dim,num_rots*data.dim);
       for(int r = 0;r<num_rots;r++)
       {
-        eff_R.block(0,dim*r,dim,dim) = 
-          R.block(0,dim*data.G(r),dim,dim);
+        eff_R.block(0,data.dim*r,data.dim,data.dim) = 
+          R.block(0,data.dim*data.G(r),data.dim,data.dim);
       }
     }
 
@@ -235,7 +244,7 @@ IGL_INLINE bool igl::arap_solve(
     VectorXd Rcol;
     columnize(eff_R,num_rots,2,Rcol);
     VectorXd Bcol = -data.K * Rcol;
-    for(int c = 0;c<dim;c++)
+    for(int c = 0;c<data.dim;c++)
     {
       VectorXd Uc,Bc,bcc,Beq;
       Bc = Bcol.block(c*n,0,n,1);
@@ -262,6 +271,6 @@ IGL_INLINE bool igl::arap_solve(
 }
 
 #ifndef IGL_HEADER_ONLY
-template bool igl::arap_precomputation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(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<int, -1, 1, 0, -1, 1> > const&, igl::ARAPData&);
 template bool igl::arap_solve<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::ARAPData&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::arap_precomputation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, igl::ARAPData&);
 #endif

include/igl/svd3x3/arap.h

@@ -31,6 +31,7 @@ namespace igl
     // M  mass matrix
     // solver_data  quadratic solver data
     // b  list of boundary indices into V
+    // dim  dimension being used for solving
     int n;
     Eigen::VectorXi G;
     ARAPEnergyType energy;
@@ -42,6 +43,7 @@ namespace igl
     Eigen::SparseMatrix<double> CSM;
     min_quad_with_fixed_data<double> solver_data;
     Eigen::VectorXi b;
+    int dim;
       ARAPData():
         n(0),
         G(),
@@ -53,7 +55,8 @@ namespace igl
         K(),
         CSM(),
         solver_data(),
-        b()
+        b(),
+        dim(-1) // force this to be set by _precomputation
     {
     };
   };
@@ -63,6 +66,8 @@ namespace igl
   // Inputs:
   //   V  #V by dim list of mesh positions
   //   F  #F by simplex-size list of triangle|tet indices into V
+  //   dim  dimension being used at solve time. For deformation usually dim =
+  //     V.cols(), for surface parameterization V.cols() = 3 and dim = 2
   //   b  #b list of "boundary" fixed vertex indices into V
   // Outputs:
   //   data  struct containing necessary precomputation
@@ -73,6 +78,7 @@ namespace igl
   IGL_INLINE bool arap_precomputation(
     const Eigen::PlainObjectBase<DerivedV> & V,
     const Eigen::PlainObjectBase<DerivedF> & F,
+    const int dim,
     const Eigen::PlainObjectBase<Derivedb> & b,
     ARAPData & data);
   // Inputs:

include/igl/svd3x3/arap_dof.cpp

@@ -123,7 +123,7 @@ IGL_INLINE bool igl::arap_dof_precomputation(
   // used to fit rotations to during optimization
   SparseMatrix<double> CSM;
   //printf("covariance_scatter_matrix()\n");
-  covariance_scatter_matrix(V,F,data.energy,CSM);
+  covariance_scatter_matrix(V,F,V.cols(),data.energy,CSM);
 #ifdef EXTREME_VERBOSE
   cout<<"CSMIJV=["<<endl;print_ijv(CSM,1);cout<<endl<<"];"<<
     endl<<"CSM=sparse(CSMIJV(:,1),CSMIJV(:,2),CSMIJV(:,3),"<<
@@ -226,7 +226,7 @@ IGL_INLINE bool igl::arap_dof_precomputation(
   // precompute arap_rhs matrix
   //printf("arap_rhs()\n");
   SparseMatrix<double> K;
-  arap_rhs(V,F,data.energy,K);
+  arap_rhs(V,F,V.cols(),data.energy,K);
 //#ifdef EXTREME_VERBOSE
 //  cout<<"KIJV=["<<endl;print_ijv(K,1);cout<<endl<<"];"<<
 //    endl<<"K=sparse(KIJV(:,1),KIJV(:,2),KIJV(:,3),"<<