slim.cpp 34 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961
  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 "boundary_loop.h"
  10. #include "cotmatrix.h"
  11. #include "edge_lengths.h"
  12. #include "grad.h"
  13. #include "local_basis.h"
  14. #include "repdiag.h"
  15. #include "vector_area_matrix.h"
  16. #include "arap.h"
  17. #include "cat.h"
  18. #include "doublearea.h"
  19. #include "grad.h"
  20. #include "local_basis.h"
  21. #include "per_face_normals.h"
  22. #include "slice_into.h"
  23. #include "volume.h"
  24. #include "polar_svd.h"
  25. #include "flip_avoiding_line_search.h"
  26. #include <iostream>
  27. #include <map>
  28. #include <set>
  29. #include <vector>
  30. #include <Eigen/IterativeLinearSolvers>
  31. #include <Eigen/SparseCholesky>
  32. #include <Eigen/IterativeLinearSolvers>
  33. #include "Timer.h"
  34. #include "sparse_cached.h"
  35. #include "AtA_cached.h"
  36. #ifdef CHOLMOD
  37. #include <Eigen/CholmodSupport>
  38. #endif
  39. namespace igl
  40. {
  41. namespace slim
  42. {
  43. // Definitions of internal functions
  44. IGL_INLINE void compute_surface_gradient_matrix(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
  45. const Eigen::MatrixXd &F1, const Eigen::MatrixXd &F2,
  46. Eigen::SparseMatrix<double> &D1, Eigen::SparseMatrix<double> &D2);
  47. IGL_INLINE void buildA(igl::SLIMData& s, std::vector<Eigen::Triplet<double> > & IJV);
  48. IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &A);
  49. IGL_INLINE void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L);
  50. IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new);
  51. IGL_INLINE double compute_soft_const_energy(igl::SLIMData& s,
  52. const Eigen::MatrixXd &V,
  53. const Eigen::MatrixXi &F,
  54. Eigen::MatrixXd &V_o);
  55. IGL_INLINE double compute_energy_with_jacobians(igl::SLIMData& s,
  56. const Eigen::MatrixXd &V,
  57. const Eigen::MatrixXi &F, const Eigen::MatrixXd &Ji,
  58. Eigen::MatrixXd &uv, Eigen::VectorXd &areas);
  59. IGL_INLINE void solve_weighted_arap(igl::SLIMData& s,
  60. const Eigen::MatrixXd &V,
  61. const Eigen::MatrixXi &F,
  62. Eigen::MatrixXd &uv,
  63. Eigen::VectorXi &soft_b_p,
  64. Eigen::MatrixXd &soft_bc_p);
  65. IGL_INLINE void update_weights_and_closest_rotations( igl::SLIMData& s,
  66. const Eigen::MatrixXd &V,
  67. const Eigen::MatrixXi &F,
  68. Eigen::MatrixXd &uv);
  69. IGL_INLINE void compute_jacobians(igl::SLIMData& s, const Eigen::MatrixXd &uv);
  70. IGL_INLINE void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix<double> &L);
  71. IGL_INLINE void pre_calc(igl::SLIMData& s);
  72. // Implementation
  73. IGL_INLINE void compute_surface_gradient_matrix(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
  74. const Eigen::MatrixXd &F1, const Eigen::MatrixXd &F2,
  75. Eigen::SparseMatrix<double> &D1, Eigen::SparseMatrix<double> &D2)
  76. {
  77. Eigen::SparseMatrix<double> G;
  78. igl::grad(V, F, G);
  79. Eigen::SparseMatrix<double> Dx = G.block(0, 0, F.rows(), V.rows());
  80. Eigen::SparseMatrix<double> Dy = G.block(F.rows(), 0, F.rows(), V.rows());
  81. Eigen::SparseMatrix<double> Dz = G.block(2 * F.rows(), 0, F.rows(), V.rows());
  82. D1 = F1.col(0).asDiagonal() * Dx + F1.col(1).asDiagonal() * Dy + F1.col(2).asDiagonal() * Dz;
  83. D2 = F2.col(0).asDiagonal() * Dx + F2.col(1).asDiagonal() * Dy + F2.col(2).asDiagonal() * Dz;
  84. }
  85. IGL_INLINE void compute_jacobians(igl::SLIMData& s, const Eigen::MatrixXd &uv)
  86. {
  87. if (s.F.cols() == 3)
  88. {
  89. // Ji=[D1*u,D2*u,D1*v,D2*v];
  90. s.Ji.col(0) = s.Dx * uv.col(0);
  91. s.Ji.col(1) = s.Dy * uv.col(0);
  92. s.Ji.col(2) = s.Dx * uv.col(1);
  93. s.Ji.col(3) = s.Dy * uv.col(1);
  94. }
  95. else /*tet mesh*/{
  96. // Ji=[D1*u,D2*u,D3*u, D1*v,D2*v, D3*v, D1*w,D2*w,D3*w];
  97. s.Ji.col(0) = s.Dx * uv.col(0);
  98. s.Ji.col(1) = s.Dy * uv.col(0);
  99. s.Ji.col(2) = s.Dz * uv.col(0);
  100. s.Ji.col(3) = s.Dx * uv.col(1);
  101. s.Ji.col(4) = s.Dy * uv.col(1);
  102. s.Ji.col(5) = s.Dz * uv.col(1);
  103. s.Ji.col(6) = s.Dx * uv.col(2);
  104. s.Ji.col(7) = s.Dy * uv.col(2);
  105. s.Ji.col(8) = s.Dz * uv.col(2);
  106. }
  107. }
  108. IGL_INLINE void update_weights_and_closest_rotations(igl::SLIMData& s,
  109. const Eigen::MatrixXd &V,
  110. const Eigen::MatrixXi &F,
  111. Eigen::MatrixXd &uv)
  112. {
  113. compute_jacobians(s, uv);
  114. const double eps = 1e-8;
  115. double exp_f = s.exp_factor;
  116. if (s.dim == 2)
  117. {
  118. for (int i = 0; i < s.Ji.rows(); ++i)
  119. {
  120. typedef Eigen::Matrix<double, 2, 2> Mat2;
  121. typedef Eigen::Matrix<double, 2, 1> Vec2;
  122. Mat2 ji, ri, ti, ui, vi;
  123. Vec2 sing;
  124. Vec2 closest_sing_vec;
  125. Mat2 mat_W;
  126. Vec2 m_sing_new;
  127. double s1, s2;
  128. ji(0, 0) = s.Ji(i, 0);
  129. ji(0, 1) = s.Ji(i, 1);
  130. ji(1, 0) = s.Ji(i, 2);
  131. ji(1, 1) = s.Ji(i, 3);
  132. igl::polar_svd(ji, ri, ti, ui, sing, vi);
  133. s1 = sing(0);
  134. s2 = sing(1);
  135. // Update Weights according to energy
  136. switch (s.slim_energy)
  137. {
  138. case igl::SLIMData::ARAP:
  139. {
  140. m_sing_new << 1, 1;
  141. break;
  142. }
  143. case igl::SLIMData::SYMMETRIC_DIRICHLET:
  144. {
  145. double s1_g = 2 * (s1 - pow(s1, -3));
  146. double s2_g = 2 * (s2 - pow(s2, -3));
  147. m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1)));
  148. break;
  149. }
  150. case igl::SLIMData::LOG_ARAP:
  151. {
  152. double s1_g = 2 * (log(s1) / s1);
  153. double s2_g = 2 * (log(s2) / s2);
  154. m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1)));
  155. break;
  156. }
  157. case igl::SLIMData::CONFORMAL:
  158. {
  159. double s1_g = 1 / (2 * s2) - s2 / (2 * pow(s1, 2));
  160. double s2_g = 1 / (2 * s1) - s1 / (2 * pow(s2, 2));
  161. double geo_avg = sqrt(s1 * s2);
  162. double s1_min = geo_avg;
  163. double s2_min = geo_avg;
  164. m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min)));
  165. // change local step
  166. closest_sing_vec << s1_min, s2_min;
  167. ri = ui * closest_sing_vec.asDiagonal() * vi.transpose();
  168. break;
  169. }
  170. case igl::SLIMData::EXP_CONFORMAL:
  171. {
  172. double s1_g = 2 * (s1 - pow(s1, -3));
  173. double s2_g = 2 * (s2 - pow(s2, -3));
  174. double geo_avg = sqrt(s1 * s2);
  175. double s1_min = geo_avg;
  176. double s2_min = geo_avg;
  177. double in_exp = exp_f * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2));
  178. double exp_thing = exp(in_exp);
  179. s1_g *= exp_thing * exp_f;
  180. s2_g *= exp_thing * exp_f;
  181. m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1)));
  182. break;
  183. }
  184. case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET:
  185. {
  186. double s1_g = 2 * (s1 - pow(s1, -3));
  187. double s2_g = 2 * (s2 - pow(s2, -3));
  188. double in_exp = exp_f * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2));
  189. double exp_thing = exp(in_exp);
  190. s1_g *= exp_thing * exp_f;
  191. s2_g *= exp_thing * exp_f;
  192. m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1)));
  193. break;
  194. }
  195. }
  196. if (std::abs(s1 - 1) < eps) m_sing_new(0) = 1;
  197. if (std::abs(s2 - 1) < eps) m_sing_new(1) = 1;
  198. mat_W = ui * m_sing_new.asDiagonal() * ui.transpose();
  199. s.W_11(i) = mat_W(0, 0);
  200. s.W_12(i) = mat_W(0, 1);
  201. s.W_21(i) = mat_W(1, 0);
  202. s.W_22(i) = mat_W(1, 1);
  203. // 2) Update local step (doesn't have to be a rotation, for instance in case of conformal energy)
  204. s.Ri(i, 0) = ri(0, 0);
  205. s.Ri(i, 1) = ri(1, 0);
  206. s.Ri(i, 2) = ri(0, 1);
  207. s.Ri(i, 3) = ri(1, 1);
  208. }
  209. }
  210. else
  211. {
  212. typedef Eigen::Matrix<double, 3, 1> Vec3;
  213. typedef Eigen::Matrix<double, 3, 3> Mat3;
  214. Mat3 ji;
  215. Vec3 m_sing_new;
  216. Vec3 closest_sing_vec;
  217. const double sqrt_2 = sqrt(2);
  218. for (int i = 0; i < s.Ji.rows(); ++i)
  219. {
  220. ji(0, 0) = s.Ji(i, 0);
  221. ji(0, 1) = s.Ji(i, 1);
  222. ji(0, 2) = s.Ji(i, 2);
  223. ji(1, 0) = s.Ji(i, 3);
  224. ji(1, 1) = s.Ji(i, 4);
  225. ji(1, 2) = s.Ji(i, 5);
  226. ji(2, 0) = s.Ji(i, 6);
  227. ji(2, 1) = s.Ji(i, 7);
  228. ji(2, 2) = s.Ji(i, 8);
  229. Mat3 ri, ti, ui, vi;
  230. Vec3 sing;
  231. igl::polar_svd(ji, ri, ti, ui, sing, vi);
  232. double s1 = sing(0);
  233. double s2 = sing(1);
  234. double s3 = sing(2);
  235. // 1) Update Weights
  236. switch (s.slim_energy)
  237. {
  238. case igl::SLIMData::ARAP:
  239. {
  240. m_sing_new << 1, 1, 1;
  241. break;
  242. }
  243. case igl::SLIMData::LOG_ARAP:
  244. {
  245. double s1_g = 2 * (log(s1) / s1);
  246. double s2_g = 2 * (log(s2) / s2);
  247. double s3_g = 2 * (log(s3) / s3);
  248. m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1)));
  249. break;
  250. }
  251. case igl::SLIMData::SYMMETRIC_DIRICHLET:
  252. {
  253. double s1_g = 2 * (s1 - pow(s1, -3));
  254. double s2_g = 2 * (s2 - pow(s2, -3));
  255. double s3_g = 2 * (s3 - pow(s3, -3));
  256. m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1)));
  257. break;
  258. }
  259. case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET:
  260. {
  261. double s1_g = 2 * (s1 - pow(s1, -3));
  262. double s2_g = 2 * (s2 - pow(s2, -3));
  263. double s3_g = 2 * (s3 - pow(s3, -3));
  264. m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1)));
  265. double in_exp = exp_f * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2));
  266. double exp_thing = exp(in_exp);
  267. s1_g *= exp_thing * exp_f;
  268. s2_g *= exp_thing * exp_f;
  269. s3_g *= exp_thing * exp_f;
  270. m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1)));
  271. break;
  272. }
  273. case igl::SLIMData::CONFORMAL:
  274. {
  275. double common_div = 9 * (pow(s1 * s2 * s3, 5. / 3.));
  276. double s1_g = (-2 * s2 * s3 * (pow(s2, 2) + pow(s3, 2) - 2 * pow(s1, 2))) / common_div;
  277. double s2_g = (-2 * s1 * s3 * (pow(s1, 2) + pow(s3, 2) - 2 * pow(s2, 2))) / common_div;
  278. double s3_g = (-2 * s1 * s2 * (pow(s1, 2) + pow(s2, 2) - 2 * pow(s3, 2))) / common_div;
  279. double closest_s = sqrt(pow(s1, 2) + pow(s3, 2)) / sqrt_2;
  280. double s1_min = closest_s;
  281. double s2_min = closest_s;
  282. double s3_min = closest_s;
  283. m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min))), sqrt(
  284. s3_g / (2 * (s3 - s3_min)));
  285. // change local step
  286. closest_sing_vec << s1_min, s2_min, s3_min;
  287. ri = ui * closest_sing_vec.asDiagonal() * vi.transpose();
  288. break;
  289. }
  290. case igl::SLIMData::EXP_CONFORMAL:
  291. {
  292. // E_conf = (s1^2 + s2^2 + s3^2)/(3*(s1*s2*s3)^(2/3) )
  293. // dE_conf/ds1 = (-2*(s2*s3)*(s2^2+s3^2 -2*s1^2) ) / (9*(s1*s2*s3)^(5/3))
  294. // Argmin E_conf(s1): s1 = sqrt(s1^2+s2^2)/sqrt(2)
  295. double common_div = 9 * (pow(s1 * s2 * s3, 5. / 3.));
  296. double s1_g = (-2 * s2 * s3 * (pow(s2, 2) + pow(s3, 2) - 2 * pow(s1, 2))) / common_div;
  297. double s2_g = (-2 * s1 * s3 * (pow(s1, 2) + pow(s3, 2) - 2 * pow(s2, 2))) / common_div;
  298. double s3_g = (-2 * s1 * s2 * (pow(s1, 2) + pow(s2, 2) - 2 * pow(s3, 2))) / common_div;
  299. double in_exp = exp_f * ((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow((s1 * s2 * s3), 2. / 3)));;
  300. double exp_thing = exp(in_exp);
  301. double closest_s = sqrt(pow(s1, 2) + pow(s3, 2)) / sqrt_2;
  302. double s1_min = closest_s;
  303. double s2_min = closest_s;
  304. double s3_min = closest_s;
  305. s1_g *= exp_thing * exp_f;
  306. s2_g *= exp_thing * exp_f;
  307. s3_g *= exp_thing * exp_f;
  308. m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min))), sqrt(
  309. s3_g / (2 * (s3 - s3_min)));
  310. // change local step
  311. closest_sing_vec << s1_min, s2_min, s3_min;
  312. ri = ui * closest_sing_vec.asDiagonal() * vi.transpose();
  313. }
  314. }
  315. if (std::abs(s1 - 1) < eps) m_sing_new(0) = 1;
  316. if (std::abs(s2 - 1) < eps) m_sing_new(1) = 1;
  317. if (std::abs(s3 - 1) < eps) m_sing_new(2) = 1;
  318. Mat3 mat_W;
  319. mat_W = ui * m_sing_new.asDiagonal() * ui.transpose();
  320. s.W_11(i) = mat_W(0, 0);
  321. s.W_12(i) = mat_W(0, 1);
  322. s.W_13(i) = mat_W(0, 2);
  323. s.W_21(i) = mat_W(1, 0);
  324. s.W_22(i) = mat_W(1, 1);
  325. s.W_23(i) = mat_W(1, 2);
  326. s.W_31(i) = mat_W(2, 0);
  327. s.W_32(i) = mat_W(2, 1);
  328. s.W_33(i) = mat_W(2, 2);
  329. // 2) Update closest rotations (not rotations in case of conformal energy)
  330. s.Ri(i, 0) = ri(0, 0);
  331. s.Ri(i, 1) = ri(1, 0);
  332. s.Ri(i, 2) = ri(2, 0);
  333. s.Ri(i, 3) = ri(0, 1);
  334. s.Ri(i, 4) = ri(1, 1);
  335. s.Ri(i, 5) = ri(2, 1);
  336. s.Ri(i, 6) = ri(0, 2);
  337. s.Ri(i, 7) = ri(1, 2);
  338. s.Ri(i, 8) = ri(2, 2);
  339. } // for loop end
  340. } // if dim end
  341. }
  342. IGL_INLINE void solve_weighted_arap(igl::SLIMData& s,
  343. const Eigen::MatrixXd &V,
  344. const Eigen::MatrixXi &F,
  345. Eigen::MatrixXd &uv,
  346. Eigen::VectorXi &soft_b_p,
  347. Eigen::MatrixXd &soft_bc_p)
  348. {
  349. using namespace Eigen;
  350. Eigen::SparseMatrix<double> L;
  351. build_linear_system(s,L);
  352. igl::Timer t;
  353. //t.start();
  354. // solve
  355. Eigen::VectorXd Uc;
  356. #ifndef CHOLMOD
  357. if (s.dim == 2)
  358. {
  359. SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
  360. Uc = solver.compute(L).solve(s.rhs);
  361. }
  362. else
  363. { // seems like CG performs much worse for 2D and way better for 3D
  364. Eigen::VectorXd guess(uv.rows() * s.dim);
  365. for (int i = 0; i < s.v_num; i++) for (int j = 0; j < s.dim; j++) guess(uv.rows() * j + i) = uv(i, j); // flatten vector
  366. ConjugateGradient<Eigen::SparseMatrix<double>, Lower | Upper> cg;
  367. cg.setTolerance(1e-8);
  368. cg.compute(L);
  369. Uc = cg.solveWithGuess(s.rhs, guess);
  370. }
  371. #else
  372. CholmodSimplicialLDLT<Eigen::SparseMatrix<double> > solver;
  373. Uc = solver.compute(L).solve(s.rhs);
  374. #endif
  375. for (int i = 0; i < s.dim; i++)
  376. uv.col(i) = Uc.block(i * s.v_n, 0, s.v_n, 1);
  377. // t.stop();
  378. // std::cerr << "solve: " << t.getElapsedTime() << std::endl;
  379. }
  380. IGL_INLINE void pre_calc(igl::SLIMData& s)
  381. {
  382. if (!s.has_pre_calc)
  383. {
  384. s.v_n = s.v_num;
  385. s.f_n = s.f_num;
  386. if (s.F.cols() == 3)
  387. {
  388. s.dim = 2;
  389. Eigen::MatrixXd F1, F2, F3;
  390. igl::local_basis(s.V, s.F, F1, F2, F3);
  391. compute_surface_gradient_matrix(s.V, s.F, F1, F2, s.Dx, s.Dy);
  392. s.W_11.resize(s.f_n);
  393. s.W_12.resize(s.f_n);
  394. s.W_21.resize(s.f_n);
  395. s.W_22.resize(s.f_n);
  396. }
  397. else
  398. {
  399. s.dim = 3;
  400. Eigen::SparseMatrix<double> G;
  401. igl::grad(s.V, s.F, G,
  402. s.mesh_improvement_3d /*use normal gradient, or one from a "regular" tet*/);
  403. s.Dx = G.block(0, 0, s.F.rows(), s.V.rows());
  404. s.Dy = G.block(s.F.rows(), 0, s.F.rows(), s.V.rows());
  405. s.Dz = G.block(2 * s.F.rows(), 0, s.F.rows(), s.V.rows());
  406. s.W_11.resize(s.f_n);
  407. s.W_12.resize(s.f_n);
  408. s.W_13.resize(s.f_n);
  409. s.W_21.resize(s.f_n);
  410. s.W_22.resize(s.f_n);
  411. s.W_23.resize(s.f_n);
  412. s.W_31.resize(s.f_n);
  413. s.W_32.resize(s.f_n);
  414. s.W_33.resize(s.f_n);
  415. }
  416. s.Dx.makeCompressed();
  417. s.Dy.makeCompressed();
  418. s.Dz.makeCompressed();
  419. s.Ri.resize(s.f_n, s.dim * s.dim);
  420. s.Ji.resize(s.f_n, s.dim * s.dim);
  421. s.rhs.resize(s.dim * s.v_num);
  422. // flattened weight matrix
  423. s.WGL_M.resize(s.dim * s.dim * s.f_n);
  424. for (int i = 0; i < s.dim * s.dim; i++)
  425. for (int j = 0; j < s.f_n; j++)
  426. s.WGL_M(i * s.f_n + j) = s.M(j);
  427. s.first_solve = true;
  428. s.has_pre_calc = true;
  429. }
  430. }
  431. IGL_INLINE void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
  432. {
  433. // formula (35) in paper
  434. std::vector<Eigen::Triplet<double> > IJV;
  435. #ifdef SLIM_CACHED
  436. buildA(s,IJV);
  437. if (s.A.rows() == 0)
  438. {
  439. s.A = Eigen::SparseMatrix<double>(s.dim * s.dim * s.f_n, s.dim * s.v_n);
  440. igl::sparse_cached_precompute(IJV,s.A_data,s.A);
  441. }
  442. else
  443. igl::sparse_cached(IJV,s.A_data,s.A);
  444. #else
  445. Eigen::SparseMatrix<double> A(s.dim * s.dim * s.f_n, s.dim * s.v_n);
  446. buildA(s,IJV);
  447. A.setFromTriplets(IJV.begin(),IJV.end());
  448. A.makeCompressed();
  449. #endif
  450. #ifdef SLIM_CACHED
  451. #else
  452. Eigen::SparseMatrix<double> At = A.transpose();
  453. At.makeCompressed();
  454. #endif
  455. #ifdef SLIM_CACHED
  456. Eigen::SparseMatrix<double> id_m(s.A.cols(), s.A.cols());
  457. #else
  458. Eigen::SparseMatrix<double> id_m(A.cols(), A.cols());
  459. #endif
  460. id_m.setIdentity();
  461. // add proximal penalty
  462. #ifdef SLIM_CACHED
  463. s.AtA_data.W = s.WGL_M;
  464. if (s.AtA.rows() == 0)
  465. igl::AtA_cached_precompute(s.A,s.AtA_data,s.AtA);
  466. else
  467. igl::AtA_cached(s.A,s.AtA_data,s.AtA);
  468. L = s.AtA + s.proximal_p * id_m; //add also a proximal
  469. L.makeCompressed();
  470. #else
  471. L = At * s.WGL_M.asDiagonal() * A + s.proximal_p * id_m; //add also a proximal term
  472. L.makeCompressed();
  473. #endif
  474. #ifdef SLIM_CACHED
  475. buildRhs(s, s.A);
  476. #else
  477. buildRhs(s, A);
  478. #endif
  479. Eigen::SparseMatrix<double> OldL = L;
  480. add_soft_constraints(s,L);
  481. L.makeCompressed();
  482. }
  483. IGL_INLINE void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
  484. {
  485. int v_n = s.v_num;
  486. for (int d = 0; d < s.dim; d++)
  487. {
  488. for (int i = 0; i < s.b.rows(); i++)
  489. {
  490. int v_idx = s.b(i);
  491. s.rhs(d * v_n + v_idx) += s.soft_const_p * s.bc(i, d); // rhs
  492. L.coeffRef(d * v_n + v_idx, d * v_n + v_idx) += s.soft_const_p; // diagonal of matrix
  493. }
  494. }
  495. }
  496. IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new)
  497. {
  498. compute_jacobians(s,V_new);
  499. return compute_energy_with_jacobians(s, s.V, s.F, s.Ji, V_new, s.M) +
  500. compute_soft_const_energy(s, s.V, s.F, V_new);
  501. }
  502. IGL_INLINE double compute_soft_const_energy(igl::SLIMData& s,
  503. const Eigen::MatrixXd &V,
  504. const Eigen::MatrixXi &F,
  505. Eigen::MatrixXd &V_o)
  506. {
  507. double e = 0;
  508. for (int i = 0; i < s.b.rows(); i++)
  509. {
  510. e += s.soft_const_p * (s.bc.row(i) - V_o.row(s.b(i))).squaredNorm();
  511. }
  512. return e;
  513. }
  514. IGL_INLINE double compute_energy_with_jacobians(igl::SLIMData& s,
  515. const Eigen::MatrixXd &V,
  516. const Eigen::MatrixXi &F, const Eigen::MatrixXd &Ji,
  517. Eigen::MatrixXd &uv, Eigen::VectorXd &areas)
  518. {
  519. double energy = 0;
  520. if (s.dim == 2)
  521. {
  522. Eigen::Matrix<double, 2, 2> ji;
  523. for (int i = 0; i < s.f_n; i++)
  524. {
  525. ji(0, 0) = Ji(i, 0);
  526. ji(0, 1) = Ji(i, 1);
  527. ji(1, 0) = Ji(i, 2);
  528. ji(1, 1) = Ji(i, 3);
  529. typedef Eigen::Matrix<double, 2, 2> Mat2;
  530. typedef Eigen::Matrix<double, 2, 1> Vec2;
  531. Mat2 ri, ti, ui, vi;
  532. Vec2 sing;
  533. igl::polar_svd(ji, ri, ti, ui, sing, vi);
  534. double s1 = sing(0);
  535. double s2 = sing(1);
  536. switch (s.slim_energy)
  537. {
  538. case igl::SLIMData::ARAP:
  539. {
  540. energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2));
  541. break;
  542. }
  543. case igl::SLIMData::SYMMETRIC_DIRICHLET:
  544. {
  545. energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2));
  546. break;
  547. }
  548. case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET:
  549. {
  550. energy += areas(i) * exp(s.exp_factor * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2)));
  551. break;
  552. }
  553. case igl::SLIMData::LOG_ARAP:
  554. {
  555. energy += areas(i) * (pow(log(s1), 2) + pow(log(s2), 2));
  556. break;
  557. }
  558. case igl::SLIMData::CONFORMAL:
  559. {
  560. energy += areas(i) * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2));
  561. break;
  562. }
  563. case igl::SLIMData::EXP_CONFORMAL:
  564. {
  565. energy += areas(i) * exp(s.exp_factor * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2)));
  566. break;
  567. }
  568. }
  569. }
  570. }
  571. else
  572. {
  573. Eigen::Matrix<double, 3, 3> ji;
  574. for (int i = 0; i < s.f_n; i++)
  575. {
  576. ji(0, 0) = Ji(i, 0);
  577. ji(0, 1) = Ji(i, 1);
  578. ji(0, 2) = Ji(i, 2);
  579. ji(1, 0) = Ji(i, 3);
  580. ji(1, 1) = Ji(i, 4);
  581. ji(1, 2) = Ji(i, 5);
  582. ji(2, 0) = Ji(i, 6);
  583. ji(2, 1) = Ji(i, 7);
  584. ji(2, 2) = Ji(i, 8);
  585. typedef Eigen::Matrix<double, 3, 3> Mat3;
  586. typedef Eigen::Matrix<double, 3, 1> Vec3;
  587. Mat3 ri, ti, ui, vi;
  588. Vec3 sing;
  589. igl::polar_svd(ji, ri, ti, ui, sing, vi);
  590. double s1 = sing(0);
  591. double s2 = sing(1);
  592. double s3 = sing(2);
  593. switch (s.slim_energy)
  594. {
  595. case igl::SLIMData::ARAP:
  596. {
  597. energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2) + pow(s3 - 1, 2));
  598. break;
  599. }
  600. case igl::SLIMData::SYMMETRIC_DIRICHLET:
  601. {
  602. energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2));
  603. break;
  604. }
  605. case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET:
  606. {
  607. energy += areas(i) * exp(s.exp_factor *
  608. (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2)));
  609. break;
  610. }
  611. case igl::SLIMData::LOG_ARAP:
  612. {
  613. energy += areas(i) * (pow(log(s1), 2) + pow(log(std::abs(s2)), 2) + pow(log(std::abs(s3)), 2));
  614. break;
  615. }
  616. case igl::SLIMData::CONFORMAL:
  617. {
  618. energy += areas(i) * ((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.)));
  619. break;
  620. }
  621. case igl::SLIMData::EXP_CONFORMAL:
  622. {
  623. energy += areas(i) * exp((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.)));
  624. break;
  625. }
  626. }
  627. }
  628. }
  629. return energy;
  630. }
  631. IGL_INLINE void buildA(igl::SLIMData& s, std::vector<Eigen::Triplet<double> > & IJV)
  632. {
  633. // formula (35) in paper
  634. if (s.dim == 2)
  635. {
  636. IJV.reserve(4 * (s.Dx.outerSize() + s.Dy.outerSize()));
  637. /*A = [W11*Dx, W12*Dx;
  638. W11*Dy, W12*Dy;
  639. W21*Dx, W22*Dx;
  640. W21*Dy, W22*Dy];*/
  641. for (int k = 0; k < s.Dx.outerSize(); ++k)
  642. {
  643. for (Eigen::SparseMatrix<double>::InnerIterator it(s.Dx, k); it; ++it)
  644. {
  645. int dx_r = it.row();
  646. int dx_c = it.col();
  647. double val = it.value();
  648. IJV.push_back(Eigen::Triplet<double>(dx_r, dx_c, val * s.W_11(dx_r)));
  649. IJV.push_back(Eigen::Triplet<double>(dx_r, s.v_n + dx_c, val * s.W_12(dx_r)));
  650. IJV.push_back(Eigen::Triplet<double>(2 * s.f_n + dx_r, dx_c, val * s.W_21(dx_r)));
  651. IJV.push_back(Eigen::Triplet<double>(2 * s.f_n + dx_r, s.v_n + dx_c, val * s.W_22(dx_r)));
  652. }
  653. }
  654. for (int k = 0; k < s.Dy.outerSize(); ++k)
  655. {
  656. for (Eigen::SparseMatrix<double>::InnerIterator it(s.Dy, k); it; ++it)
  657. {
  658. int dy_r = it.row();
  659. int dy_c = it.col();
  660. double val = it.value();
  661. IJV.push_back(Eigen::Triplet<double>(s.f_n + dy_r, dy_c, val * s.W_11(dy_r)));
  662. IJV.push_back(Eigen::Triplet<double>(s.f_n + dy_r, s.v_n + dy_c, val * s.W_12(dy_r)));
  663. IJV.push_back(Eigen::Triplet<double>(3 * s.f_n + dy_r, dy_c, val * s.W_21(dy_r)));
  664. IJV.push_back(Eigen::Triplet<double>(3 * s.f_n + dy_r, s.v_n + dy_c, val * s.W_22(dy_r)));
  665. }
  666. }
  667. }
  668. else
  669. {
  670. /*A = [W11*Dx, W12*Dx, W13*Dx;
  671. W11*Dy, W12*Dy, W13*Dy;
  672. W11*Dz, W12*Dz, W13*Dz;
  673. W21*Dx, W22*Dx, W23*Dx;
  674. W21*Dy, W22*Dy, W23*Dy;
  675. W21*Dz, W22*Dz, W23*Dz;
  676. W31*Dx, W32*Dx, W33*Dx;
  677. W31*Dy, W32*Dy, W33*Dy;
  678. W31*Dz, W32*Dz, W33*Dz;];*/
  679. IJV.reserve(9 * (s.Dx.outerSize() + s.Dy.outerSize() + s.Dz.outerSize()));
  680. for (int k = 0; k < s.Dx.outerSize(); k++)
  681. {
  682. for (Eigen::SparseMatrix<double>::InnerIterator it(s.Dx, k); it; ++it)
  683. {
  684. int dx_r = it.row();
  685. int dx_c = it.col();
  686. double val = it.value();
  687. IJV.push_back(Eigen::Triplet<double>(dx_r, dx_c, val * s.W_11(dx_r)));
  688. IJV.push_back(Eigen::Triplet<double>(dx_r, s.v_n + dx_c, val * s.W_12(dx_r)));
  689. IJV.push_back(Eigen::Triplet<double>(dx_r, 2 * s.v_n + dx_c, val * s.W_13(dx_r)));
  690. IJV.push_back(Eigen::Triplet<double>(3 * s.f_n + dx_r, dx_c, val * s.W_21(dx_r)));
  691. IJV.push_back(Eigen::Triplet<double>(3 * s.f_n + dx_r, s.v_n + dx_c, val * s.W_22(dx_r)));
  692. IJV.push_back(Eigen::Triplet<double>(3 * s.f_n + dx_r, 2 * s.v_n + dx_c, val * s.W_23(dx_r)));
  693. IJV.push_back(Eigen::Triplet<double>(6 * s.f_n + dx_r, dx_c, val * s.W_31(dx_r)));
  694. IJV.push_back(Eigen::Triplet<double>(6 * s.f_n + dx_r, s.v_n + dx_c, val * s.W_32(dx_r)));
  695. IJV.push_back(Eigen::Triplet<double>(6 * s.f_n + dx_r, 2 * s.v_n + dx_c, val * s.W_33(dx_r)));
  696. }
  697. }
  698. for (int k = 0; k < s.Dy.outerSize(); k++)
  699. {
  700. for (Eigen::SparseMatrix<double>::InnerIterator it(s.Dy, k); it; ++it)
  701. {
  702. int dy_r = it.row();
  703. int dy_c = it.col();
  704. double val = it.value();
  705. IJV.push_back(Eigen::Triplet<double>(s.f_n + dy_r, dy_c, val * s.W_11(dy_r)));
  706. IJV.push_back(Eigen::Triplet<double>(s.f_n + dy_r, s.v_n + dy_c, val * s.W_12(dy_r)));
  707. IJV.push_back(Eigen::Triplet<double>(s.f_n + dy_r, 2 * s.v_n + dy_c, val * s.W_13(dy_r)));
  708. IJV.push_back(Eigen::Triplet<double>(4 * s.f_n + dy_r, dy_c, val * s.W_21(dy_r)));
  709. IJV.push_back(Eigen::Triplet<double>(4 * s.f_n + dy_r, s.v_n + dy_c, val * s.W_22(dy_r)));
  710. IJV.push_back(Eigen::Triplet<double>(4 * s.f_n + dy_r, 2 * s.v_n + dy_c, val * s.W_23(dy_r)));
  711. IJV.push_back(Eigen::Triplet<double>(7 * s.f_n + dy_r, dy_c, val * s.W_31(dy_r)));
  712. IJV.push_back(Eigen::Triplet<double>(7 * s.f_n + dy_r, s.v_n + dy_c, val * s.W_32(dy_r)));
  713. IJV.push_back(Eigen::Triplet<double>(7 * s.f_n + dy_r, 2 * s.v_n + dy_c, val * s.W_33(dy_r)));
  714. }
  715. }
  716. for (int k = 0; k < s.Dz.outerSize(); k++)
  717. {
  718. for (Eigen::SparseMatrix<double>::InnerIterator it(s.Dz, k); it; ++it)
  719. {
  720. int dz_r = it.row();
  721. int dz_c = it.col();
  722. double val = it.value();
  723. IJV.push_back(Eigen::Triplet<double>(2 * s.f_n + dz_r, dz_c, val * s.W_11(dz_r)));
  724. IJV.push_back(Eigen::Triplet<double>(2 * s.f_n + dz_r, s.v_n + dz_c, val * s.W_12(dz_r)));
  725. IJV.push_back(Eigen::Triplet<double>(2 * s.f_n + dz_r, 2 * s.v_n + dz_c, val * s.W_13(dz_r)));
  726. IJV.push_back(Eigen::Triplet<double>(5 * s.f_n + dz_r, dz_c, val * s.W_21(dz_r)));
  727. IJV.push_back(Eigen::Triplet<double>(5 * s.f_n + dz_r, s.v_n + dz_c, val * s.W_22(dz_r)));
  728. IJV.push_back(Eigen::Triplet<double>(5 * s.f_n + dz_r, 2 * s.v_n + dz_c, val * s.W_23(dz_r)));
  729. IJV.push_back(Eigen::Triplet<double>(8 * s.f_n + dz_r, dz_c, val * s.W_31(dz_r)));
  730. IJV.push_back(Eigen::Triplet<double>(8 * s.f_n + dz_r, s.v_n + dz_c, val * s.W_32(dz_r)));
  731. IJV.push_back(Eigen::Triplet<double>(8 * s.f_n + dz_r, 2 * s.v_n + dz_c, val * s.W_33(dz_r)));
  732. }
  733. }
  734. }
  735. }
  736. IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &A)
  737. {
  738. Eigen::VectorXd f_rhs(s.dim * s.dim * s.f_n);
  739. f_rhs.setZero();
  740. if (s.dim == 2)
  741. {
  742. /*b = [W11*R11 + W12*R21; (formula (36))
  743. W11*R12 + W12*R22;
  744. W21*R11 + W22*R21;
  745. W21*R12 + W22*R22];*/
  746. for (int i = 0; i < s.f_n; i++)
  747. {
  748. f_rhs(i + 0 * s.f_n) = s.W_11(i) * s.Ri(i, 0) + s.W_12(i) * s.Ri(i, 1);
  749. f_rhs(i + 1 * s.f_n) = s.W_11(i) * s.Ri(i, 2) + s.W_12(i) * s.Ri(i, 3);
  750. f_rhs(i + 2 * s.f_n) = s.W_21(i) * s.Ri(i, 0) + s.W_22(i) * s.Ri(i, 1);
  751. f_rhs(i + 3 * s.f_n) = s.W_21(i) * s.Ri(i, 2) + s.W_22(i) * s.Ri(i, 3);
  752. }
  753. }
  754. else
  755. {
  756. /*b = [W11*R11 + W12*R21 + W13*R31;
  757. W11*R12 + W12*R22 + W13*R32;
  758. W11*R13 + W12*R23 + W13*R33;
  759. W21*R11 + W22*R21 + W23*R31;
  760. W21*R12 + W22*R22 + W23*R32;
  761. W21*R13 + W22*R23 + W23*R33;
  762. W31*R11 + W32*R21 + W33*R31;
  763. W31*R12 + W32*R22 + W33*R32;
  764. W31*R13 + W32*R23 + W33*R33;];*/
  765. for (int i = 0; i < s.f_n; i++)
  766. {
  767. f_rhs(i + 0 * s.f_n) = s.W_11(i) * s.Ri(i, 0) + s.W_12(i) * s.Ri(i, 1) + s.W_13(i) * s.Ri(i, 2);
  768. f_rhs(i + 1 * s.f_n) = s.W_11(i) * s.Ri(i, 3) + s.W_12(i) * s.Ri(i, 4) + s.W_13(i) * s.Ri(i, 5);
  769. f_rhs(i + 2 * s.f_n) = s.W_11(i) * s.Ri(i, 6) + s.W_12(i) * s.Ri(i, 7) + s.W_13(i) * s.Ri(i, 8);
  770. f_rhs(i + 3 * s.f_n) = s.W_21(i) * s.Ri(i, 0) + s.W_22(i) * s.Ri(i, 1) + s.W_23(i) * s.Ri(i, 2);
  771. f_rhs(i + 4 * s.f_n) = s.W_21(i) * s.Ri(i, 3) + s.W_22(i) * s.Ri(i, 4) + s.W_23(i) * s.Ri(i, 5);
  772. f_rhs(i + 5 * s.f_n) = s.W_21(i) * s.Ri(i, 6) + s.W_22(i) * s.Ri(i, 7) + s.W_23(i) * s.Ri(i, 8);
  773. f_rhs(i + 6 * s.f_n) = s.W_31(i) * s.Ri(i, 0) + s.W_32(i) * s.Ri(i, 1) + s.W_33(i) * s.Ri(i, 2);
  774. f_rhs(i + 7 * s.f_n) = s.W_31(i) * s.Ri(i, 3) + s.W_32(i) * s.Ri(i, 4) + s.W_33(i) * s.Ri(i, 5);
  775. f_rhs(i + 8 * s.f_n) = s.W_31(i) * s.Ri(i, 6) + s.W_32(i) * s.Ri(i, 7) + s.W_33(i) * s.Ri(i, 8);
  776. }
  777. }
  778. Eigen::VectorXd uv_flat(s.dim *s.v_n);
  779. for (int i = 0; i < s.dim; i++)
  780. for (int j = 0; j < s.v_n; j++)
  781. uv_flat(s.v_n * i + j) = s.V_o(j, i);
  782. s.rhs = (f_rhs.transpose() * s.WGL_M.asDiagonal() * A).transpose() + s.proximal_p * uv_flat;
  783. }
  784. }
  785. }
  786. /// Slim Implementation
  787. IGL_INLINE void igl::slim_precompute(
  788. const Eigen::MatrixXd &V,
  789. const Eigen::MatrixXi &F,
  790. const Eigen::MatrixXd &V_init,
  791. SLIMData &data,
  792. SLIMData::SLIM_ENERGY slim_energy,
  793. Eigen::VectorXi &b,
  794. Eigen::MatrixXd &bc,
  795. double soft_p)
  796. {
  797. data.V = V;
  798. data.F = F;
  799. data.V_o = V_init;
  800. data.v_num = V.rows();
  801. data.f_num = F.rows();
  802. data.slim_energy = slim_energy;
  803. data.b = b;
  804. data.bc = bc;
  805. data.soft_const_p = soft_p;
  806. data.proximal_p = 0.0001;
  807. igl::doublearea(V, F, data.M);
  808. data.M /= 2.;
  809. data.mesh_area = data.M.sum();
  810. data.mesh_improvement_3d = false; // whether to use a jacobian derived from a real mesh or an abstract regular mesh (used for mesh improvement)
  811. data.exp_factor = 1.0; // param used only for exponential energies (e.g exponential symmetric dirichlet)
  812. assert (F.cols() == 3 || F.cols() == 4);
  813. igl::slim::pre_calc(data);
  814. data.energy = igl::slim::compute_energy(data,data.V_o) / data.mesh_area;
  815. }
  816. IGL_INLINE Eigen::MatrixXd igl::slim_solve(SLIMData &data, int iter_num)
  817. {
  818. for (int i = 0; i < iter_num; i++)
  819. {
  820. Eigen::MatrixXd dest_res;
  821. dest_res = data.V_o;
  822. // Solve Weighted Proxy
  823. igl::slim::update_weights_and_closest_rotations(data,data.V, data.F, dest_res);
  824. igl::slim::solve_weighted_arap(data,data.V, data.F, dest_res, data.b, data.bc);
  825. double old_energy = data.energy;
  826. std::function<double(Eigen::MatrixXd &)> compute_energy = [&](
  827. Eigen::MatrixXd &aaa) { return igl::slim::compute_energy(data,aaa); };
  828. data.energy = igl::flip_avoiding_line_search(data.F, data.V_o, dest_res, compute_energy,
  829. data.energy * data.mesh_area) / data.mesh_area;
  830. }
  831. return data.V_o;
  832. }
  833. #ifdef IGL_STATIC_LIBRARY
  834. // Explicit template instantiation
  835. #endif