slim.cpp 30 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852
  1. // This file is part of libigl, a simple c++ geometry processing library.
  2. //
  3. // Copyright (C) 2016 Michael Rabinovich
  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. #include "slim.h"
  9. #include <igl/boundary_loop.h>
  10. #include <igl/cotmatrix.h>
  11. #include <igl/edge_lengths.h>
  12. #include <igl/grad.h>
  13. #include <igl/local_basis.h>
  14. #include <igl/readOBJ.h>
  15. #include <igl/repdiag.h>
  16. #include <igl/vector_area_matrix.h>
  17. #include <iostream>
  18. #include <Eigen/Dense>
  19. #include <Eigen/Sparse>
  20. #include <map>
  21. #include <set>
  22. #include <vector>
  23. #include "igl/arap.h"
  24. #include "igl/cat.h"
  25. #include "igl/doublearea.h"
  26. #include "igl/grad.h"
  27. #include "igl/local_basis.h"
  28. #include "igl/per_face_normals.h"
  29. #include "igl/slice_into.h"
  30. #include "igl/volume.h"
  31. #include "igl/polar_svd.h"
  32. #include <Eigen/IterativeLinearSolvers>
  33. #include <Eigen/Sparse>
  34. #include <Eigen/SparseQR>
  35. #include <Eigen/SparseCholesky>
  36. #include <Eigen/IterativeLinearSolvers>
  37. #include "flip_avoiding_line_search.h"
  38. using namespace std;
  39. using namespace Eigen;
  40. ///////// Helper functions to compute gradient matrices
  41. void compute_surface_gradient_matrix(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
  42. const Eigen::MatrixXd& F1, const Eigen::MatrixXd& F2,
  43. Eigen::SparseMatrix<double>& D1, Eigen::SparseMatrix<double>& D2) {
  44. Eigen::SparseMatrix<double> G;
  45. igl::grad(V,F,G);
  46. Eigen::SparseMatrix<double> Dx = G.block(0,0,F.rows(),V.rows());
  47. Eigen::SparseMatrix<double> Dy = G.block(F.rows(),0,F.rows(),V.rows());
  48. Eigen::SparseMatrix<double> Dz = G.block(2*F.rows(),0,F.rows(),V.rows());
  49. D1 = F1.col(0).asDiagonal()*Dx + F1.col(1).asDiagonal()*Dy + F1.col(2).asDiagonal()*Dz;
  50. D2 = F2.col(0).asDiagonal()*Dx + F2.col(1).asDiagonal()*Dy + F2.col(2).asDiagonal()*Dz;
  51. }
  52. void compute_tet_grad_matrix(const Eigen::MatrixXd& V, const Eigen::MatrixXi& T,
  53. Eigen::SparseMatrix<double>& Dx, Eigen::SparseMatrix<double>& Dy, Eigen::SparseMatrix<double>& Dz, bool uniform) {
  54. using namespace Eigen;
  55. assert(T.cols() == 4);
  56. const int n = V.rows(); int m = T.rows();
  57. /*
  58. F = [ ...
  59. T(:,1) T(:,2) T(:,3); ...
  60. T(:,1) T(:,3) T(:,4); ...
  61. T(:,1) T(:,4) T(:,2); ...
  62. T(:,2) T(:,4) T(:,3)]; */
  63. MatrixXi F(4*m,3);
  64. for (int i = 0; i < m; i++) {
  65. F.row(0*m + i) << T(i,0), T(i,1), T(i,2);
  66. F.row(1*m + i) << T(i,0), T(i,2), T(i,3);
  67. F.row(2*m + i) << T(i,0), T(i,3), T(i,1);
  68. F.row(3*m + i) << T(i,1), T(i,3), T(i,2);
  69. }
  70. // compute volume of each tet
  71. VectorXd vol; igl::volume(V,T,vol);
  72. VectorXd A(F.rows());
  73. MatrixXd N(F.rows(),3);
  74. if (!uniform) {
  75. // compute tetrahedron face normals
  76. igl::per_face_normals(V,F,N); int norm_rows = N.rows();
  77. for (int i = 0; i < norm_rows; i++)
  78. N.row(i) /= N.row(i).norm();
  79. igl::doublearea(V,F,A); A/=2.;
  80. } else {
  81. // Use a uniform tetrahedra as a reference:
  82. // V = h*[0,0,0;1,0,0;0.5,sqrt(3)/2.,0;0.5,sqrt(3)/6.,sqrt(2)/sqrt(3)]
  83. //
  84. // With normals
  85. // 0 0 1.0000
  86. // 0.8165 -0.4714 -0.3333
  87. // 0 0.9428 -0.3333
  88. // -0.8165 -0.4714 -0.3333
  89. for (int i = 0; i < m; i++) {
  90. N.row(0*m+i) << 0,0,1;
  91. double a = sqrt(2)*std::cbrt(3*vol(i));
  92. A(0*m+i) = (pow(a,2)*sqrt(3))/4.;
  93. }
  94. for (int i = 0; i < m; i++) {
  95. N.row(1*m+i) << 0.8165,-0.4714,-0.3333;
  96. double a = sqrt(2)*std::cbrt(3*vol(i));
  97. A(1*m+i) = (pow(a,2)*sqrt(3))/4.;
  98. }
  99. for (int i = 0; i < m; i++) {
  100. N.row(2*m+i) << 0,0.9428,-0.3333;
  101. double a = sqrt(2)*std::cbrt(3*vol(i));
  102. A(2*m+i) = (pow(a,2)*sqrt(3))/4.;
  103. }
  104. for (int i = 0; i < m; i++) {
  105. N.row(3*m+i) << -0.8165,-0.4714,-0.3333;
  106. double a = sqrt(2)*std::cbrt(3*vol(i));
  107. A(3*m+i) = (pow(a,2)*sqrt(3))/4.;
  108. }
  109. }
  110. /* G = sparse( ...
  111. [0*m + repmat(1:m,1,4) ...
  112. 1*m + repmat(1:m,1,4) ...
  113. 2*m + repmat(1:m,1,4)], ...
  114. repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1), ...
  115. repmat(A./(3*repmat(vol,4,1)),3,1).*N(:), ...
  116. 3*m,n);*/
  117. std::vector<Triplet<double> > Dx_t,Dy_t,Dz_t;
  118. for (int i = 0; i < 4*m; i++) {
  119. int T_j; // j indexes : repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1)
  120. switch (i/m) {
  121. case 0:
  122. T_j = 3;
  123. break;
  124. case 1:
  125. T_j = 1;
  126. break;
  127. case 2:
  128. T_j = 2;
  129. break;
  130. case 3:
  131. T_j = 0;
  132. break;
  133. }
  134. int i_idx = i%m;
  135. int j_idx = T(i_idx,T_j);
  136. double val_before_n = A(i)/(3*vol(i_idx));
  137. Dx_t.push_back(Triplet<double>(i_idx, j_idx, val_before_n * N(i,0)));
  138. Dy_t.push_back(Triplet<double>(i_idx, j_idx, val_before_n * N(i,1)));
  139. Dz_t.push_back(Triplet<double>(i_idx, j_idx, val_before_n * N(i,2)));
  140. }
  141. Dx.resize(m,n); Dy.resize(m,n); Dz.resize(m,n);
  142. Dx.setFromTriplets(Dx_t.begin(), Dx_t.end());
  143. Dy.setFromTriplets(Dy_t.begin(), Dy_t.end());
  144. Dz.setFromTriplets(Dz_t.begin(), Dz_t.end());
  145. }
  146. // Computes the weights and solve the linear system for the quadratic proxy specified in the paper
  147. // The output of this is used to generate a search direction that will be fed to the Linesearch class
  148. class WeightedGlobalLocal {
  149. public:
  150. WeightedGlobalLocal(igl::SLIMData& state);
  151. // Compute necessary information before solving the proxy quadratic
  152. void pre_calc();
  153. // Solve the weighted proxy global step
  154. // Output:
  155. // V_new #V by dim list of mesh positions (will be fed to a linesearch algorithm)
  156. void solve_weighted_proxy(Eigen::MatrixXd& V_new);
  157. // Compute the energy specified in the SLIMData structure + the soft constraint energy (in case there are soft constraints)
  158. // Input:
  159. // V_new #V by dim list of mesh positions
  160. virtual double compute_energy(Eigen::MatrixXd& V_new);
  161. //private:
  162. void compute_jacobians(const Eigen::MatrixXd& V_o);
  163. double compute_energy_with_jacobians(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
  164. const Eigen::MatrixXd& Ji, Eigen::MatrixXd& V_o, Eigen::VectorXd& areas);
  165. double compute_soft_const_energy(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
  166. Eigen::MatrixXd& V_o);
  167. void update_weights_and_closest_rotations(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, Eigen::MatrixXd& uv);
  168. void solve_weighted_arap(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, Eigen::MatrixXd& uv, Eigen::VectorXi& b,
  169. Eigen::MatrixXd& bc);
  170. void build_linear_system(Eigen::SparseMatrix<double> &L);
  171. void buildA(Eigen::SparseMatrix<double>& A);
  172. void buildRhs(const Eigen::SparseMatrix<double>& At);
  173. void add_soft_constraints(Eigen::SparseMatrix<double> &L);
  174. void add_proximal_penalty();
  175. igl::SLIMData& m_state;
  176. Eigen::VectorXd M;
  177. Eigen::VectorXd rhs;
  178. Eigen::MatrixXd Ri,Ji;
  179. Eigen::VectorXd W_11; Eigen::VectorXd W_12; Eigen::VectorXd W_13;
  180. Eigen::VectorXd W_21; Eigen::VectorXd W_22; Eigen::VectorXd W_23;
  181. Eigen::VectorXd W_31; Eigen::VectorXd W_32; Eigen::VectorXd W_33;
  182. Eigen::SparseMatrix<double> Dx,Dy,Dz;
  183. int f_n,v_n;
  184. bool first_solve;
  185. bool has_pre_calc = false;
  186. int dim;
  187. };
  188. //// Implementation
  189. WeightedGlobalLocal::WeightedGlobalLocal(igl::SLIMData& state) :
  190. m_state(state) {
  191. }
  192. void WeightedGlobalLocal::solve_weighted_proxy(Eigen::MatrixXd& V_new) {
  193. update_weights_and_closest_rotations(m_state.V,m_state.F,V_new);
  194. solve_weighted_arap(m_state.V,m_state.F,V_new,m_state.b,m_state.bc);
  195. }
  196. void WeightedGlobalLocal::compute_jacobians(const Eigen::MatrixXd& uv) {
  197. if (m_state.F.cols() == 3){
  198. // Ji=[D1*u,D2*u,D1*v,D2*v];
  199. Ji.col(0) = Dx*uv.col(0); Ji.col(1) = Dy*uv.col(0);
  200. Ji.col(2) = Dx*uv.col(1); Ji.col(3) = Dy*uv.col(1);
  201. } else /*tet mesh*/{
  202. // Ji=[D1*u,D2*u,D3*u, D1*v,D2*v, D3*v, D1*w,D2*w,D3*w];
  203. Ji.col(0) = Dx*uv.col(0); Ji.col(1) = Dy*uv.col(0); Ji.col(2) = Dz*uv.col(0);
  204. Ji.col(3) = Dx*uv.col(1); Ji.col(4) = Dy*uv.col(1); Ji.col(5) = Dz*uv.col(1);
  205. Ji.col(6) = Dx*uv.col(2); Ji.col(7) = Dy*uv.col(2); Ji.col(8) = Dz*uv.col(2);
  206. }
  207. }
  208. void WeightedGlobalLocal::update_weights_and_closest_rotations(const Eigen::MatrixXd& V,
  209. const Eigen::MatrixXi& F, Eigen::MatrixXd& uv) {
  210. compute_jacobians(uv);
  211. const double eps = 1e-8;
  212. double exp_f = m_state.exp_factor;
  213. if (dim==2) {
  214. for(int i=0; i <Ji.rows(); ++i ) {
  215. typedef Eigen::Matrix<double,2,2> Mat2;
  216. typedef Eigen::Matrix<double,2,1> Vec2;
  217. Mat2 ji,ri,ti,ui,vi; Vec2 sing; Vec2 closest_sing_vec;Mat2 mat_W;
  218. Vec2 m_sing_new;
  219. double s1,s2;
  220. ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1);
  221. ji(1,0) = Ji(i,2); ji(1,1) = Ji(i,3);
  222. igl::polar_svd(ji,ri,ti,ui,sing,vi);
  223. s1 = sing(0); s2 = sing(1);
  224. // Update Weights according to energy
  225. switch(m_state.slim_energy) {
  226. case igl::SLIMData::ARAP: {
  227. m_sing_new << 1,1;
  228. break;
  229. } case igl::SLIMData::SYMMETRIC_DIRICHLET: {
  230. double s1_g = 2* (s1-pow(s1,-3));
  231. double s2_g = 2 * (s2-pow(s2,-3));
  232. m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
  233. break;
  234. } case igl::SLIMData::LOG_ARAP: {
  235. double s1_g = 2 * (log(s1)/s1);
  236. double s2_g = 2 * (log(s2)/s2);
  237. m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
  238. break;
  239. } case igl::SLIMData::CONFORMAL: {
  240. double s1_g = 1/(2*s2) - s2/(2*pow(s1,2));
  241. double s2_g = 1/(2*s1) - s1/(2*pow(s2,2));
  242. double geo_avg = sqrt(s1*s2);
  243. double s1_min = geo_avg; double s2_min = geo_avg;
  244. m_sing_new << sqrt(s1_g/(2*(s1-s1_min))), sqrt(s2_g/(2*(s2-s2_min)));
  245. // change local step
  246. closest_sing_vec << s1_min,s2_min;
  247. ri = ui*closest_sing_vec.asDiagonal()*vi.transpose();
  248. break;
  249. } case igl::SLIMData::EXP_CONFORMAL: {
  250. double s1_g = 2* (s1-pow(s1,-3));
  251. double s2_g = 2 * (s2-pow(s2,-3));
  252. double geo_avg = sqrt(s1*s2);
  253. double s1_min = geo_avg; double s2_min = geo_avg;
  254. double in_exp = exp_f*((pow(s1,2)+pow(s2,2))/(2*s1*s2));
  255. double exp_thing = exp(in_exp);
  256. s1_g *= exp_thing*exp_f;
  257. s2_g *= exp_thing*exp_f;
  258. m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
  259. break;
  260. } case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
  261. double s1_g = 2* (s1-pow(s1,-3));
  262. double s2_g = 2 * (s2-pow(s2,-3));
  263. double in_exp = exp_f*(pow(s1,2)+pow(s1,-2)+pow(s2,2)+pow(s2,-2));
  264. double exp_thing = exp(in_exp);
  265. s1_g *= exp_thing*exp_f;
  266. s2_g *= exp_thing*exp_f;
  267. m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
  268. break;
  269. }
  270. }
  271. if (abs(s1-1) < eps) m_sing_new(0) = 1; if (abs(s2-1) < eps) m_sing_new(1) = 1;
  272. mat_W = ui*m_sing_new.asDiagonal()*ui.transpose();
  273. W_11(i) = mat_W(0,0); W_12(i) = mat_W(0,1); W_21(i) = mat_W(1,0); W_22(i) = mat_W(1,1);
  274. // 2) Update local step (doesn't have to be a rotation, for instance in case of conformal energy)
  275. Ri(i,0) = ri(0,0); Ri(i,1) = ri(1,0); Ri(i,2) = ri(0,1); Ri(i,3) = ri(1,1);
  276. }
  277. } else {
  278. typedef Eigen::Matrix<double,3,1> Vec3; typedef Eigen::Matrix<double,3,3> Mat3;
  279. Mat3 ji; Vec3 m_sing_new; Vec3 closest_sing_vec;
  280. const double sqrt_2 = sqrt(2);
  281. for(int i=0; i <Ji.rows(); ++i ) {
  282. ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1); ji(0,2) = Ji(i,2);
  283. ji(1,0) = Ji(i,3); ji(1,1) = Ji(i,4); ji(1,2) = Ji(i,5);
  284. ji(2,0) = Ji(i,6); ji(2,1) = Ji(i,7); ji(2,2) = Ji(i,8);
  285. Mat3 ri,ti,ui,vi;
  286. Vec3 sing;
  287. igl::polar_svd(ji,ri,ti,ui,sing,vi);
  288. double s1 = sing(0); double s2 = sing(1); double s3 = sing(2);
  289. // 1) Update Weights
  290. switch(m_state.slim_energy) {
  291. case igl::SLIMData::ARAP: {
  292. m_sing_new << 1,1,1;
  293. break;
  294. } case igl::SLIMData::LOG_ARAP: {
  295. double s1_g = 2 * (log(s1)/s1);
  296. double s2_g = 2 * (log(s2)/s2);
  297. double s3_g = 2 * (log(s3)/s3);
  298. m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
  299. break;
  300. } case igl::SLIMData::SYMMETRIC_DIRICHLET: {
  301. double s1_g = 2* (s1-pow(s1,-3));
  302. double s2_g = 2 * (s2-pow(s2,-3));
  303. double s3_g = 2 * (s3-pow(s3,-3));
  304. m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
  305. break;
  306. } case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
  307. double s1_g = 2* (s1-pow(s1,-3));
  308. double s2_g = 2 * (s2-pow(s2,-3));
  309. double s3_g = 2 * (s3-pow(s3,-3));
  310. m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
  311. double in_exp = exp_f*(pow(s1,2)+pow(s1,-2)+pow(s2,2)+pow(s2,-2)+pow(s3,2)+pow(s3,-2));
  312. double exp_thing = exp(in_exp);
  313. s1_g *= exp_thing*exp_f;
  314. s2_g *= exp_thing*exp_f;
  315. s3_g *= exp_thing*exp_f;
  316. m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
  317. break;
  318. }
  319. case igl::SLIMData::CONFORMAL: {
  320. double common_div = 9*(pow(s1*s2*s3,5./3.));
  321. double s1_g = (-2*s2*s3*(pow(s2,2)+pow(s3,2)-2*pow(s1,2)) ) / common_div;
  322. double s2_g = (-2*s1*s3*(pow(s1,2)+pow(s3,2)-2*pow(s2,2)) ) / common_div;
  323. double s3_g = (-2*s1*s2*(pow(s1,2)+pow(s2,2)-2*pow(s3,2)) ) / common_div;
  324. double closest_s = sqrt(pow(s1,2)+pow(s3,2)) / sqrt_2;
  325. double s1_min = closest_s; double s2_min = closest_s; double s3_min = closest_s;
  326. m_sing_new << sqrt(s1_g/(2*(s1-s1_min))), sqrt(s2_g/(2*(s2-s2_min))), sqrt(s3_g/(2*(s3-s3_min)));
  327. // change local step
  328. closest_sing_vec << s1_min,s2_min,s3_min;
  329. ri = ui*closest_sing_vec.asDiagonal()*vi.transpose();
  330. break;
  331. }
  332. case igl::SLIMData::EXP_CONFORMAL: {
  333. // E_conf = (s1^2 + s2^2 + s3^2)/(3*(s1*s2*s3)^(2/3) )
  334. // dE_conf/ds1 = (-2*(s2*s3)*(s2^2+s3^2 -2*s1^2) ) / (9*(s1*s2*s3)^(5/3))
  335. // Argmin E_conf(s1): s1 = sqrt(s1^2+s2^2)/sqrt(2)
  336. double common_div = 9*(pow(s1*s2*s3,5./3.));
  337. double s1_g = (-2*s2*s3*(pow(s2,2)+pow(s3,2)-2*pow(s1,2)) ) / common_div;
  338. double s2_g = (-2*s1*s3*(pow(s1,2)+pow(s3,2)-2*pow(s2,2)) ) / common_div;
  339. double s3_g = (-2*s1*s2*(pow(s1,2)+pow(s2,2)-2*pow(s3,2)) ) / common_div;
  340. double in_exp = exp_f*( (pow(s1,2)+pow(s2,2)+pow(s3,2))/ (3*pow((s1*s2*s3),2./3)) ); ;
  341. double exp_thing = exp(in_exp);
  342. double closest_s = sqrt(pow(s1,2)+pow(s3,2)) / sqrt_2;
  343. double s1_min = closest_s; double s2_min = closest_s; double s3_min = closest_s;
  344. s1_g *= exp_thing*exp_f;
  345. s2_g *= exp_thing*exp_f;
  346. s3_g *= exp_thing*exp_f;
  347. m_sing_new << sqrt(s1_g/(2*(s1-s1_min))), sqrt(s2_g/(2*(s2-s2_min))), sqrt(s3_g/(2*(s3-s3_min)));
  348. // change local step
  349. closest_sing_vec << s1_min,s2_min,s3_min;
  350. ri = ui*closest_sing_vec.asDiagonal()*vi.transpose();
  351. }
  352. }
  353. if (abs(s1-1) < eps) m_sing_new(0) = 1; if (abs(s2-1) < eps) m_sing_new(1) = 1; if (abs(s3-1) < eps) m_sing_new(2) = 1;
  354. Mat3 mat_W;
  355. mat_W = ui*m_sing_new.asDiagonal()*ui.transpose();
  356. W_11(i) = mat_W(0,0);
  357. W_12(i) = mat_W(0,1);
  358. W_13(i) = mat_W(0,2);
  359. W_21(i) = mat_W(1,0);
  360. W_22(i) = mat_W(1,1);
  361. W_23(i) = mat_W(1,2);
  362. W_31(i) = mat_W(2,0);
  363. W_32(i) = mat_W(2,1);
  364. W_33(i) = mat_W(2,2);
  365. // 2) Update closest rotations (not rotations in case of conformal energy)
  366. Ri(i,0) = ri(0,0); Ri(i,1) = ri(1,0); Ri(i,2) = ri(2,0);
  367. Ri(i,3) = ri(0,1); Ri(i,4) = ri(1,1); Ri(i,5) = ri(2,1);
  368. Ri(i,6) = ri(0,2); Ri(i,7) = ri(1,2); Ri(i,8) = ri(2,2);
  369. } // for loop end
  370. } // if dim end
  371. }
  372. void WeightedGlobalLocal::solve_weighted_arap(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
  373. Eigen::MatrixXd& uv, Eigen::VectorXi& soft_b_p, Eigen::MatrixXd& soft_bc_p) {
  374. using namespace Eigen;
  375. Eigen::SparseMatrix<double> L;
  376. build_linear_system(L);
  377. // solve
  378. Eigen::VectorXd Uc;
  379. if (dim == 2) {
  380. SimplicialLDLT<SparseMatrix<double> > solver;
  381. Uc = solver.compute(L).solve(rhs);
  382. } else { // seems like CG performs much worse for 2D and way better for 3D
  383. Eigen::VectorXd guess(uv.rows()*dim);
  384. for (int i = 0; i < dim; i++) for (int j = 0; j < dim; j++) guess(uv.rows()*i + j) = uv(i,j); // flatten vector
  385. ConjugateGradient<SparseMatrix<double>, Eigen::Upper> solver;
  386. solver.setTolerance(1e-8);
  387. Uc = solver.compute(L).solveWithGuess(rhs,guess);
  388. }
  389. for (int i = 0; i < dim; i++)
  390. uv.col(i) = Uc.block(i*v_n,0,v_n,1);
  391. }
  392. void WeightedGlobalLocal::pre_calc() {
  393. if (!has_pre_calc) {
  394. v_n = m_state.v_num; f_n = m_state.f_num;
  395. if (m_state.F.cols() == 3) {
  396. dim = 2;
  397. Eigen::MatrixXd F1,F2,F3;
  398. igl::local_basis(m_state.V,m_state.F,F1,F2,F3);
  399. compute_surface_gradient_matrix(m_state.V,m_state.F,F1,F2,Dx,Dy);
  400. W_11.resize(f_n); W_12.resize(f_n); W_21.resize(f_n); W_22.resize(f_n);
  401. } else {
  402. dim = 3;
  403. compute_tet_grad_matrix(m_state.V,m_state.F,Dx,Dy,Dz,
  404. m_state.mesh_improvement_3d /*use normal gradient, or one from a "regular" tet*/);
  405. W_11.resize(f_n);W_12.resize(f_n);W_13.resize(f_n);
  406. W_21.resize(f_n);W_22.resize(f_n);W_23.resize(f_n);
  407. W_31.resize(f_n);W_32.resize(f_n);W_33.resize(f_n);
  408. }
  409. Dx.makeCompressed();Dy.makeCompressed(); Dz.makeCompressed();
  410. Ri.resize(f_n, dim*dim); Ji.resize(f_n, dim*dim);
  411. rhs.resize(dim*m_state.v_num);
  412. // flattened weight matrix
  413. M.resize(dim*dim*f_n);
  414. for (int i = 0; i < dim*dim; i++)
  415. for (int j = 0; j < f_n; j++)
  416. M(i*f_n + j) = m_state.M(j);
  417. first_solve = true;
  418. has_pre_calc = true;
  419. }
  420. }
  421. void WeightedGlobalLocal::build_linear_system(Eigen::SparseMatrix<double> &L) {
  422. // formula (35) in paper
  423. Eigen::SparseMatrix<double> A(dim*dim*f_n, dim*v_n);
  424. buildA(A);
  425. Eigen::SparseMatrix<double> At = A.transpose();
  426. At.makeCompressed();
  427. Eigen::SparseMatrix<double> id_m(At.rows(),At.rows()); id_m.setIdentity();
  428. // add proximal penalty
  429. L = At*M.asDiagonal()*A + m_state.proximal_p * id_m; //add also a proximal term
  430. L.makeCompressed();
  431. buildRhs(At);
  432. Eigen::SparseMatrix<double> OldL = L;
  433. add_soft_constraints(L);
  434. L.makeCompressed();
  435. }
  436. void WeightedGlobalLocal::add_soft_constraints(Eigen::SparseMatrix<double> &L) {
  437. int v_n = m_state.v_num;
  438. for (int d = 0; d < dim; d++) {
  439. for (int i = 0; i < m_state.b.rows(); i++) {
  440. int v_idx = m_state.b(i);
  441. rhs(d*v_n + v_idx) += m_state.soft_const_p * m_state.bc(i,d); // rhs
  442. L.coeffRef(d*v_n + v_idx, d*v_n + v_idx) += m_state.soft_const_p; // diagonal of matrix
  443. }
  444. }
  445. }
  446. double WeightedGlobalLocal::compute_energy(Eigen::MatrixXd& V_new) {
  447. compute_jacobians(V_new);
  448. return compute_energy_with_jacobians(m_state.V,m_state.F, Ji, V_new,m_state.M) + compute_soft_const_energy(m_state.V,m_state.F,V_new);
  449. }
  450. double WeightedGlobalLocal::compute_soft_const_energy(const Eigen::MatrixXd& V,
  451. const Eigen::MatrixXi& F,
  452. Eigen::MatrixXd& V_o) {
  453. double e = 0;
  454. for (int i = 0; i < m_state.b.rows(); i++) {
  455. e += m_state.soft_const_p*(m_state.bc.row(i)-V_o.row(m_state.b(i))).squaredNorm();
  456. }
  457. return e;
  458. }
  459. double WeightedGlobalLocal::compute_energy_with_jacobians(const Eigen::MatrixXd& V,
  460. const Eigen::MatrixXi& F, const Eigen::MatrixXd& Ji, Eigen::MatrixXd& uv, Eigen::VectorXd& areas) {
  461. double energy = 0;
  462. if (dim == 2) {
  463. Eigen::Matrix<double,2,2> ji;
  464. for (int i = 0; i < f_n; i++) {
  465. ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1);
  466. ji(1,0) = Ji(i,2); ji(1,1) = Ji(i,3);
  467. typedef Eigen::Matrix<double,2,2> Mat2;
  468. typedef Eigen::Matrix<double,2,1> Vec2;
  469. Mat2 ri,ti,ui,vi; Vec2 sing;
  470. igl::polar_svd(ji,ri,ti,ui,sing,vi);
  471. double s1 = sing(0); double s2 = sing(1);
  472. switch(m_state.slim_energy) {
  473. case igl::SLIMData::ARAP: {
  474. energy+= areas(i) * (pow(s1-1,2) + pow(s2-1,2));
  475. break;
  476. }
  477. case igl::SLIMData::SYMMETRIC_DIRICHLET: {
  478. energy += areas(i) * (pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2));
  479. break;
  480. }
  481. case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
  482. energy += areas(i) * exp(m_state.exp_factor*(pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2)));
  483. break;
  484. }
  485. case igl::SLIMData::LOG_ARAP: {
  486. energy += areas(i) * (pow(log(s1),2) + pow(log(s2),2));
  487. break;
  488. }
  489. case igl::SLIMData::CONFORMAL: {
  490. energy += areas(i) * ( (pow(s1,2)+pow(s2,2))/(2*s1*s2) );
  491. break;
  492. }
  493. case igl::SLIMData::EXP_CONFORMAL: {
  494. energy += areas(i) * exp(m_state.exp_factor*((pow(s1,2)+pow(s2,2))/(2*s1*s2)));
  495. break;
  496. }
  497. }
  498. }
  499. } else {
  500. Eigen::Matrix<double,3,3> ji;
  501. for (int i = 0; i < f_n; i++) {
  502. ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1); ji(0,2) = Ji(i,2);
  503. ji(1,0) = Ji(i,3); ji(1,1) = Ji(i,4); ji(1,2) = Ji(i,5);
  504. ji(2,0) = Ji(i,6); ji(2,1) = Ji(i,7); ji(2,2) = Ji(i,8);
  505. typedef Eigen::Matrix<double,3,3> Mat3;
  506. typedef Eigen::Matrix<double,3,1> Vec3;
  507. Mat3 ri,ti,ui,vi; Vec3 sing;
  508. igl::polar_svd(ji,ri,ti,ui,sing,vi);
  509. double s1 = sing(0); double s2 = sing(1); double s3 = sing(2);
  510. switch(m_state.slim_energy) {
  511. case igl::SLIMData::ARAP: {
  512. energy+= areas(i) * (pow(s1-1,2) + pow(s2-1,2) + pow(s3-1,2));
  513. break;
  514. }
  515. case igl::SLIMData::SYMMETRIC_DIRICHLET: {
  516. energy += areas(i) * (pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2) + pow(s3,2) + pow(s3,-2));
  517. break;
  518. }
  519. case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
  520. energy += areas(i) * exp(m_state.exp_factor*(pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2) + pow(s3,2) + pow(s3,-2)));
  521. break;
  522. }
  523. case igl::SLIMData::LOG_ARAP: {
  524. energy += areas(i) * (pow(log(s1),2) + pow(log(abs(s2)),2) + pow(log(abs(s3)),2));
  525. break;
  526. }
  527. case igl::SLIMData::CONFORMAL: {
  528. energy += areas(i) * ( ( pow(s1,2)+pow(s2,2)+pow(s3,2) ) /(3*pow(s1*s2*s3,2./3.)) );
  529. break;
  530. }
  531. case igl::SLIMData::EXP_CONFORMAL: {
  532. energy += areas(i) * exp( ( pow(s1,2)+pow(s2,2)+pow(s3,2) ) /(3*pow(s1*s2*s3,2./3.)) );
  533. break;
  534. }
  535. }
  536. }
  537. }
  538. return energy;
  539. }
  540. void WeightedGlobalLocal::buildA(Eigen::SparseMatrix<double>& A) {
  541. // formula (35) in paper
  542. std::vector<Triplet<double> > IJV;
  543. if (dim == 2) {
  544. IJV.reserve(4*(Dx.outerSize()+ Dy.outerSize()));
  545. /*A = [W11*Dx, W12*Dx;
  546. W11*Dy, W12*Dy;
  547. W21*Dx, W22*Dx;
  548. W21*Dy, W22*Dy];*/
  549. for (int k=0; k<Dx.outerSize(); ++k) {
  550. for (SparseMatrix<double>::InnerIterator it(Dx,k); it; ++it) {
  551. int dx_r = it.row();
  552. int dx_c = it.col();
  553. double val = it.value();
  554. IJV.push_back(Triplet<double>(dx_r,dx_c, val*W_11(dx_r)));
  555. IJV.push_back(Triplet<double>(dx_r,v_n + dx_c, val*W_12(dx_r)));
  556. IJV.push_back(Triplet<double>(2*f_n+dx_r,dx_c, val*W_21(dx_r)));
  557. IJV.push_back(Triplet<double>(2*f_n+dx_r,v_n + dx_c, val*W_22(dx_r)));
  558. }
  559. }
  560. for (int k=0; k<Dy.outerSize(); ++k) {
  561. for (SparseMatrix<double>::InnerIterator it(Dy,k); it; ++it) {
  562. int dy_r = it.row();
  563. int dy_c = it.col();
  564. double val = it.value();
  565. IJV.push_back(Triplet<double>(f_n+dy_r,dy_c, val*W_11(dy_r)));
  566. IJV.push_back(Triplet<double>(f_n+dy_r,v_n + dy_c, val*W_12(dy_r)));
  567. IJV.push_back(Triplet<double>(3*f_n+dy_r,dy_c, val*W_21(dy_r)));
  568. IJV.push_back(Triplet<double>(3*f_n+dy_r,v_n + dy_c, val*W_22(dy_r)));
  569. }
  570. }
  571. } else {
  572. /*A = [W11*Dx, W12*Dx, W13*Dx;
  573. W11*Dy, W12*Dy, W13*Dy;
  574. W11*Dz, W12*Dz, W13*Dz;
  575. W21*Dx, W22*Dx, W23*Dx;
  576. W21*Dy, W22*Dy, W23*Dy;
  577. W21*Dz, W22*Dz, W23*Dz;
  578. W31*Dx, W32*Dx, W33*Dx;
  579. W31*Dy, W32*Dy, W33*Dy;
  580. W31*Dz, W32*Dz, W33*Dz;];*/
  581. IJV.reserve(9*(Dx.outerSize()+ Dy.outerSize() + Dz.outerSize()));
  582. for (int k = 0; k < Dx.outerSize(); k++) {
  583. for (SparseMatrix<double>::InnerIterator it(Dx,k); it; ++it) {
  584. int dx_r = it.row();
  585. int dx_c = it.col();
  586. double val = it.value();
  587. IJV.push_back(Triplet<double>(dx_r,dx_c, val*W_11(dx_r)));
  588. IJV.push_back(Triplet<double>(dx_r,v_n + dx_c, val*W_12(dx_r)));
  589. IJV.push_back(Triplet<double>(dx_r,2*v_n + dx_c, val*W_13(dx_r)));
  590. IJV.push_back(Triplet<double>(3*f_n+dx_r,dx_c, val*W_21(dx_r)));
  591. IJV.push_back(Triplet<double>(3*f_n+dx_r,v_n + dx_c, val*W_22(dx_r)));
  592. IJV.push_back(Triplet<double>(3*f_n+dx_r,2*v_n + dx_c, val*W_23(dx_r)));
  593. IJV.push_back(Triplet<double>(6*f_n+dx_r,dx_c, val*W_31(dx_r)));
  594. IJV.push_back(Triplet<double>(6*f_n+dx_r,v_n + dx_c, val*W_32(dx_r)));
  595. IJV.push_back(Triplet<double>(6*f_n+dx_r,2*v_n + dx_c, val*W_33(dx_r)));
  596. }
  597. }
  598. for (int k = 0; k < Dy.outerSize(); k++) {
  599. for (SparseMatrix<double>::InnerIterator it(Dy,k); it; ++it) {
  600. int dy_r = it.row();
  601. int dy_c = it.col();
  602. double val = it.value();
  603. IJV.push_back(Triplet<double>(f_n+dy_r,dy_c, val*W_11(dy_r)));
  604. IJV.push_back(Triplet<double>(f_n+dy_r,v_n + dy_c, val*W_12(dy_r)));
  605. IJV.push_back(Triplet<double>(f_n+dy_r,2*v_n + dy_c, val*W_13(dy_r)));
  606. IJV.push_back(Triplet<double>(4*f_n+dy_r,dy_c, val*W_21(dy_r)));
  607. IJV.push_back(Triplet<double>(4*f_n+dy_r,v_n + dy_c, val*W_22(dy_r)));
  608. IJV.push_back(Triplet<double>(4*f_n+dy_r,2*v_n + dy_c, val*W_23(dy_r)));
  609. IJV.push_back(Triplet<double>(7*f_n+dy_r,dy_c, val*W_31(dy_r)));
  610. IJV.push_back(Triplet<double>(7*f_n+dy_r,v_n + dy_c, val*W_32(dy_r)));
  611. IJV.push_back(Triplet<double>(7*f_n+dy_r,2*v_n + dy_c, val*W_33(dy_r)));
  612. }
  613. }
  614. for (int k = 0; k < Dz.outerSize(); k++) {
  615. for (SparseMatrix<double>::InnerIterator it(Dz,k); it; ++it) {
  616. int dz_r = it.row();
  617. int dz_c = it.col();
  618. double val = it.value();
  619. IJV.push_back(Triplet<double>(2*f_n + dz_r,dz_c, val*W_11(dz_r)));
  620. IJV.push_back(Triplet<double>(2*f_n + dz_r,v_n + dz_c, val*W_12(dz_r)));
  621. IJV.push_back(Triplet<double>(2*f_n + dz_r,2*v_n + dz_c, val*W_13(dz_r)));
  622. IJV.push_back(Triplet<double>(5*f_n+dz_r,dz_c, val*W_21(dz_r)));
  623. IJV.push_back(Triplet<double>(5*f_n+dz_r,v_n + dz_c, val*W_22(dz_r)));
  624. IJV.push_back(Triplet<double>(5*f_n+dz_r,2*v_n + dz_c, val*W_23(dz_r)));
  625. IJV.push_back(Triplet<double>(8*f_n+dz_r,dz_c, val*W_31(dz_r)));
  626. IJV.push_back(Triplet<double>(8*f_n+dz_r,v_n + dz_c, val*W_32(dz_r)));
  627. IJV.push_back(Triplet<double>(8*f_n+dz_r,2*v_n + dz_c, val*W_33(dz_r)));
  628. }
  629. }
  630. }
  631. A.setFromTriplets(IJV.begin(),IJV.end());
  632. }
  633. void WeightedGlobalLocal::buildRhs(const Eigen::SparseMatrix<double>& At) {
  634. VectorXd f_rhs(dim*dim*f_n); f_rhs.setZero();
  635. if (dim==2) {
  636. /*b = [W11*R11 + W12*R21; (formula (36))
  637. W11*R12 + W12*R22;
  638. W21*R11 + W22*R21;
  639. W21*R12 + W22*R22];*/
  640. for (int i = 0; i < f_n; i++) {
  641. f_rhs(i+0*f_n) = W_11(i) * Ri(i,0) + W_12(i)*Ri(i,1);
  642. f_rhs(i+1*f_n) = W_11(i) * Ri(i,2) + W_12(i)*Ri(i,3);
  643. f_rhs(i+2*f_n) = W_21(i) * Ri(i,0) + W_22(i)*Ri(i,1);
  644. f_rhs(i+3*f_n) = W_21(i) * Ri(i,2) + W_22(i)*Ri(i,3);
  645. }
  646. } else {
  647. /*b = [W11*R11 + W12*R21 + W13*R31;
  648. W11*R12 + W12*R22 + W13*R32;
  649. W11*R13 + W12*R23 + W13*R33;
  650. W21*R11 + W22*R21 + W23*R31;
  651. W21*R12 + W22*R22 + W23*R32;
  652. W21*R13 + W22*R23 + W23*R33;
  653. W31*R11 + W32*R21 + W33*R31;
  654. W31*R12 + W32*R22 + W33*R32;
  655. W31*R13 + W32*R23 + W33*R33;];*/
  656. for (int i = 0; i < f_n; i++) {
  657. f_rhs(i+0*f_n) = W_11(i) * Ri(i,0) + W_12(i)*Ri(i,1) + W_13(i)*Ri(i,2);
  658. f_rhs(i+1*f_n) = W_11(i) * Ri(i,3) + W_12(i)*Ri(i,4) + W_13(i)*Ri(i,5);
  659. f_rhs(i+2*f_n) = W_11(i) * Ri(i,6) + W_12(i)*Ri(i,7) + W_13(i)*Ri(i,8);
  660. f_rhs(i+3*f_n) = W_21(i) * Ri(i,0) + W_22(i)*Ri(i,1) + W_23(i)*Ri(i,2);
  661. f_rhs(i+4*f_n) = W_21(i) * Ri(i,3) + W_22(i)*Ri(i,4) + W_23(i)*Ri(i,5);
  662. f_rhs(i+5*f_n) = W_21(i) * Ri(i,6) + W_22(i)*Ri(i,7) + W_23(i)*Ri(i,8);
  663. f_rhs(i+6*f_n) = W_31(i) * Ri(i,0) + W_32(i)*Ri(i,1) + W_33(i)*Ri(i,2);
  664. f_rhs(i+7*f_n) = W_31(i) * Ri(i,3) + W_32(i)*Ri(i,4) + W_33(i)*Ri(i,5);
  665. f_rhs(i+8*f_n) = W_31(i) * Ri(i,6) + W_32(i)*Ri(i,7) + W_33(i)*Ri(i,8);
  666. }
  667. }
  668. VectorXd uv_flat(dim*v_n);
  669. for (int i = 0; i < dim; i++)
  670. for (int j = 0; j < v_n; j++)
  671. uv_flat(v_n*i+j) = m_state.V_o(j,i);
  672. rhs = (At*M.asDiagonal()*f_rhs + m_state.proximal_p * uv_flat);
  673. }
  674. // #define TwoPi 6.28318530717958648
  675. // const double eps=1e-14;
  676. /// Slim Implementation
  677. IGL_INLINE void igl::slim_precompute(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::MatrixXd& V_init, SLIMData& data,
  678. SLIMData::SLIM_ENERGY slim_energy, Eigen::VectorXi& b, Eigen::MatrixXd& bc, double soft_p) {
  679. data.V = V;
  680. data.F = F;
  681. data.V_o = V_init;
  682. data.v_num = V.rows();
  683. data.f_num = F.rows();
  684. data.slim_energy = slim_energy;
  685. data.b = b;
  686. data.bc = bc;
  687. data.soft_const_p = soft_p;
  688. data.proximal_p = 0.0001;
  689. igl::doublearea(V,F,data.M); data.M /= 2.;
  690. data.mesh_area = data.M.sum();
  691. data.mesh_improvement_3d = false; // whether to use a jacobian derived from a real mesh or an abstract regular mesh (used for mesh improvement)
  692. data.exp_factor = 1.0; // param used only for exponential energies (e.g exponential symmetric dirichlet)
  693. assert (F.cols() == 3 || F.cols() == 4);
  694. data.wGlobalLocal = new WeightedGlobalLocal(data);
  695. data.wGlobalLocal->pre_calc();
  696. data.energy = data.wGlobalLocal->compute_energy(data.V_o)/data.mesh_area;
  697. }
  698. IGL_INLINE void igl::slim_solve(SLIMData& data, int iter_num) {
  699. for (int i = 0; i < iter_num; i++) {
  700. Eigen::MatrixXd dest_res;
  701. dest_res = data.V_o;
  702. data.wGlobalLocal->solve_weighted_proxy(dest_res);
  703. double old_energy = data.energy;
  704. std::function<double(Eigen::MatrixXd&)> compute_energy = [&](Eigen::MatrixXd& aaa) { return data.wGlobalLocal->compute_energy(aaa); };
  705. data.energy = igl::flip_avoiding_line_search(data.F,data.V_o, dest_res, compute_energy,
  706. data.energy*data.mesh_area)/data.mesh_area;
  707. }
  708. }