1 // This file is part of libigl, a simple c++ geometry processing library. 2 // 3 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com> 4 // 5 // This Source Code Form is subject to the terms of the Mozilla Public License 6 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 7 // obtain one at http://mozilla.org/MPL/2.0/. 8 #ifndef IGL_ARAP_ENERGY_TYPE_DOF_H 9 #define IGL_ARAP_ENERGY_TYPE_DOF_H 10 #include "igl_inline.h" 11 12 #include <Eigen/Dense> 13 #include <Eigen/Sparse> 14 #include "ARAPEnergyType.h" 15 #include <vector> 16 17 namespace igl 18 { 19 // Caller example: 20 // 21 // Once: 22 // arap_dof_precomputation(...) 23 // 24 // Each frame: 25 // while(not satisfied) 26 // arap_dof_update(...) 27 // end 28 29 template <typename LbsMatrixType, typename SSCALAR> 30 struct ArapDOFData; 31 32 /////////////////////////////////////////////////////////////////////////// 33 // 34 // Arap DOF precomputation consists of two parts the computation. The first is 35 // that which depends solely on the mesh (V,F), the linear blend skinning 36 // weights (M) and the groups G. Then there's the part that depends on the 37 // previous precomputation and the list of free and fixed vertices. 38 // 39 /////////////////////////////////////////////////////////////////////////// 40 41 42 // The code and variables differ from the description in Section 3 of "Fast 43 // Automatic Skinning Transformations" by [Jacobson et al. 2012] 44 // 45 // Here is a useful conversion table: 46 // 47 // [article] [code] 48 // S = \tilde{K} T S = CSM * Lsep 49 // S --> R S --> R --shuffled--> Rxyz 50 // Gamma_solve RT = Pi_1 \tilde{K} RT L_part1xyz = CSolveBlock1 * Rxyz 51 // Pi_1 \tilde{K} CSolveBlock1 52 // Peq = [T_full; P_pos] 53 // T_full B_eq_fix <--- L0 54 // P_pos B_eq 55 // Pi_2 * P_eq = Lpart2and3 = Lpart2 + Lpart3 56 // Pi_2_left T_full + Lpart3 = M_fullsolve(right) * B_eq_fix 57 // Pi_2_right P_pos Lpart2 = M_fullsolve(left) * B_eq 58 // T = [Pi_1 Pi_2] [\tilde{K}TRT P_eq] L = Lpart1 + Lpart2and3 59 // 60 61 // Precomputes the system we are going to optimize. This consists of building 62 // constructor matrices (to compute covariance matrices from transformations 63 // and to build the poisson solve right hand side from rotation matrix entries) 64 // and also prefactoring the poisson system. 65 // 66 // Inputs: 67 // V #V by dim list of vertex positions 68 // F #F by {3|4} list of face indices 69 // M #V * dim by #handles * dim * (dim+1) matrix such that 70 // new_V(:) = LBS(V,W,A) = reshape(M * A,size(V)), where A is a column 71 // vectors formed by the entries in each handle's dim by dim+1 72 // transformation matrix. Specifcally, A = 73 // reshape(permute(Astack,[3 1 2]),n*dim*(dim+1),1) 74 // or A = [Lxx;Lyx;Lxy;Lyy;tx;ty], and likewise for other dim 75 // if Astack(:,:,i) is the dim by (dim+1) transformation at handle i 76 // handles are ordered according to P then BE (point handles before bone 77 // handles) 78 // G #V list of group indices (1 to k) for each vertex, such that vertex i 79 // is assigned to group G(i) 80 // Outputs: 81 // data structure containing all necessary precomputation for calling 82 // arap_dof_update 83 // Returns true on success, false on error 84 // 85 // See also: lbs_matrix_column 86 template <typename LbsMatrixType, typename SSCALAR> 87 IGL_INLINE bool arap_dof_precomputation( 88 const Eigen::MatrixXd & V, 89 const Eigen::MatrixXi & F, 90 const LbsMatrixType & M, 91 const Eigen::Matrix<int,Eigen::Dynamic,1> & G, 92 ArapDOFData<LbsMatrixType, SSCALAR> & data); 93 94 // Should always be called after arap_dof_precomputation, but may be called in 95 // between successive calls to arap_dof_update, recomputes precomputation 96 // given that there are only changes in free and fixed 97 // 98 // Inputs: 99 // fixed_dim list of transformation element indices for fixed (or partailly 100 // fixed) handles: not necessarily the complement of 'free' 101 // NOTE: the constraints for fixed transformations still need to be 102 // present in A_eq 103 // A_eq dim*#constraint_points by m*dim*(dim+1) matrix of linear equality 104 // constraint coefficients. Each row corresponds to a linear constraint, 105 // so that A_eq * L = Beq says that the linear transformation entries in 106 // the column L should produce the user supplied positional constraints 107 // for each handle in Beq. The row A_eq(i*dim+d) corresponds to the 108 // constrain on coordinate d of position i 109 // Outputs: 110 // data structure containing all necessary precomputation for calling 111 // arap_dof_update 112 // Returns true on success, false on error 113 // 114 // See also: lbs_matrix_column 115 template <typename LbsMatrixType, typename SSCALAR> 116 IGL_INLINE bool arap_dof_recomputation( 117 const Eigen::Matrix<int,Eigen::Dynamic,1> & fixed_dim, 118 const Eigen::SparseMatrix<double> & A_eq, 119 ArapDOFData<LbsMatrixType, SSCALAR> & data); 120 121 // Optimizes the transformations attached to each weight function based on 122 // precomputed system. 123 // 124 // Inputs: 125 // data precomputation data struct output from arap_dof_precomputation 126 // Beq dim*#constraint_points constraint values. 127 // L0 #handles * dim * dim+1 list of initial guess transformation entries, 128 // also holds fixed transformation entries for fixed handles 129 // max_iters maximum number of iterations 130 // tol stopping criteria parameter. If variables (linear transformation 131 // matrix entries) change by less than 'tol' the optimization terminates, 132 // 0.75 (weak tolerance) 133 // 0.0 (extreme tolerance) 134 // Outputs: 135 // L #handles * dim * dim+1 list of final optimized transformation entries, 136 // allowed to be the same as L 137 template <typename LbsMatrixType, typename SSCALAR> 138 IGL_INLINE bool arap_dof_update( 139 const ArapDOFData<LbsMatrixType,SSCALAR> & data, 140 const Eigen::Matrix<double,Eigen::Dynamic,1> & B_eq, 141 const Eigen::MatrixXd & L0, 142 const int max_iters, 143 const double tol, 144 Eigen::MatrixXd & L 145 ); 146 147 // Structure that contains fields for all precomputed data or data that needs 148 // to be remembered at update 149 template <typename LbsMatrixType, typename SSCALAR> 150 struct ArapDOFData 151 { 152 typedef Eigen::Matrix<SSCALAR, Eigen::Dynamic, Eigen::Dynamic> MatrixXS; 153 // Type of arap energy we're solving 154 igl::ARAPEnergyType energy; 155 //// LU decomposition precomptation data; note: not used by araf_dop_update 156 //// any more, replaced by M_FullSolve 157 //igl::min_quad_with_fixed_data<double> lu_data; 158 // List of indices of fixed transformation entries 159 Eigen::Matrix<int,Eigen::Dynamic,1> fixed_dim; 160 // List of precomputed covariance scatter matrices multiplied by lbs 161 // matrices 162 //std::vector<Eigen::SparseMatrix<double> > CSM_M; 163 std::vector<Eigen::MatrixXd> CSM_M; 164 LbsMatrixType M_KG; 165 // Number of mesh vertices 166 int n; 167 // Number of weight functions 168 int m; 169 // Number of dimensions 170 int dim; 171 // Effective dimensions 172 int effective_dim; 173 // List of indices into C of positional constraints 174 Eigen::Matrix<int,Eigen::Dynamic,1> interpolated; 175 std::vector<bool> free_mask; 176 // Full quadratic coefficients matrix before lagrangian (should be dense) 177 LbsMatrixType Q; 178 179 180 //// Solve matrix for the global step 181 //Eigen::MatrixXd M_Solve; // TODO: remove from here 182 183 // Full solve matrix that contains also conversion from rotations to the right hand side, 184 // i.e., solves Poisson transformations just from rotations and positional constraints 185 MatrixXS M_FullSolve; 186 187 // Precomputed condensed matrices (3x3 commutators folded to 1x1): 188 MatrixXS CSM; 189 MatrixXS CSolveBlock1; 190 191 // Print timings at each update 192 bool print_timings; 193 194 // Dynamics 195 bool with_dynamics; 196 // I'm hiding the extra dynamics stuff in this struct, which sort of defeats 197 // the purpose of this function-based coding style... 198 199 // Time step 200 double h; 201 202 // L0 #handles * dim * dim+1 list of transformation entries from 203 // previous solve 204 MatrixXS L0; 205 //// Lm1 #handles * dim * dim+1 list of transformation entries from 206 //// previous-previous solve 207 //MatrixXS Lm1; 208 // "Velocity" 209 MatrixXS Lvel0; 210 211 // #V by dim matrix of external forces 212 // fext 213 MatrixXS fext; 214 215 // Mass_tilde: MT * Mass * M 216 LbsMatrixType Mass_tilde; 217 218 // Force due to gravity (premultiplier) 219 Eigen::MatrixXd fgrav; 220 // Direction of gravity 221 Eigen::Vector3d grav_dir; 222 // Magnitude of gravity 223 double grav_mag; 224 225 // Π1 from the paper 226 MatrixXS Pi_1; 227 228 // Default values ArapDOFDataArapDOFData229 ArapDOFData(): 230 energy(igl::ARAP_ENERGY_TYPE_SPOKES), 231 with_dynamics(false), 232 h(1), 233 grav_dir(0,-1,0), 234 grav_mag(0) 235 { 236 } 237 }; 238 } 239 240 #ifndef IGL_STATIC_LIBRARY 241 # include "arap_dof.cpp" 242 #endif 243 244 #endif 245