Blender V4.3
slim.cpp
Go to the documentation of this file.
1/* SPDX-FileCopyrightText: 2016 Michael Rabinovich
2 * 2023 Blender Authors
3 *
4 * SPDX-License-Identifier: MPL-2.0 */
5
10#include "slim.h"
11#include "doublearea.h"
12#include "flip_avoiding_line_search.h"
13
14#include "BLI_assert.h"
15#include "BLI_math_base.h" /* M_PI */
16
17#include <vector>
18
19#include <Eigen/Geometry>
20#include <Eigen/IterativeLinearSolvers>
21#include <Eigen/SVD>
22#include <Eigen/SparseCholesky>
23
24namespace slim {
25
26/* GRAD
27 * G = grad(V,F)
28 *
29 * Compute the numerical gradient operator
30 *
31 * Inputs:
32 * V #vertices by 3 list of mesh vertex positions
33 * F #faces by 3 list of mesh face indices
34 * uniform #boolean (default false) - Use a uniform mesh instead of the vertices V
35 * Outputs:
36 * G #faces*dim by #V Gradient operator
37 *
38 *
39 * Gradient of a scalar function defined on piecewise linear elements (mesh)
40 * is constant on each triangle i,j,k:
41 * grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
42 * where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
43 * i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
44 * 90 degrees
45 */
46template<typename DerivedV, typename DerivedF>
47static inline void grad(const Eigen::PlainObjectBase<DerivedV> &V,
48 const Eigen::PlainObjectBase<DerivedF> &F,
49 Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
50 bool uniform = false)
51{
52 Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 3> eperp21(F.rows(), 3),
53 eperp13(F.rows(), 3);
54
55 for (int i = 0; i < F.rows(); ++i) {
56 /* Renaming indices of vertices of triangles for convenience. */
57 int i1 = F(i, 0);
58 int i2 = F(i, 1);
59 int i3 = F(i, 2);
60
61 /* #F x 3 matrices of triangle edge vectors, named after opposite vertices. */
62 Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v32 = V.row(i3) - V.row(i2);
63 Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v13 = V.row(i1) - V.row(i3);
64 Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v21 = V.row(i2) - V.row(i1);
65 Eigen::Matrix<typename DerivedV::Scalar, 1, 3> n = v32.cross(v13);
66 /* Area of parallelogram is twice area of triangle.
67 * Area of parallelogram is || v1 x v2 ||.
68 * This does correct l2 norm of rows, so that it contains #F list of twice.
69 * triangle areas. */
70 double dblA = std::sqrt(n.dot(n));
71 Eigen::Matrix<typename DerivedV::Scalar, 1, 3> u;
72 if (!uniform) {
73 /* Now normalize normals to get unit normals. */
74 u = n / dblA;
75 }
76 else {
77 /* Abstract equilateral triangle v1=(0,0), v2=(h,0), v3=(h/2, (sqrt(3)/2)*h) */
78
79 /* Get h (by the area of the triangle). */
80 double h = sqrt((dblA) /
81 sin(M_PI / 3.0)); /* (h^2*sin(60))/2. = Area => h = sqrt(2*Area/sin_60) */
82
83 Eigen::VectorXd v1, v2, v3;
84 v1 << 0, 0, 0;
85 v2 << h, 0, 0;
86 v3 << h / 2., (sqrt(3) / 2.) * h, 0;
87
88 /* Now fix v32,v13,v21 and the normal. */
89 v32 = v3 - v2;
90 v13 = v1 - v3;
91 v21 = v2 - v1;
92 n = v32.cross(v13);
93 }
94
95 /* Rotate each vector 90 degrees around normal. */
96 double norm21 = std::sqrt(v21.dot(v21));
97 double norm13 = std::sqrt(v13.dot(v13));
98 eperp21.row(i) = u.cross(v21);
99 eperp21.row(i) = eperp21.row(i) / std::sqrt(eperp21.row(i).dot(eperp21.row(i)));
100 eperp21.row(i) *= norm21 / dblA;
101 eperp13.row(i) = u.cross(v13);
102 eperp13.row(i) = eperp13.row(i) / std::sqrt(eperp13.row(i).dot(eperp13.row(i)));
103 eperp13.row(i) *= norm13 / dblA;
104 }
105
106 std::vector<int> rs;
107 rs.reserve(F.rows() * 4 * 3);
108 std::vector<int> cs;
109 cs.reserve(F.rows() * 4 * 3);
110 std::vector<double> vs;
111 vs.reserve(F.rows() * 4 * 3);
112
113 /* Row indices. */
114 for (int r = 0; r < 3; r++) {
115 for (int j = 0; j < 4; j++) {
116 for (int i = r * F.rows(); i < (r + 1) * F.rows(); i++)
117 rs.push_back(i);
118 }
119 }
120
121 /* Column indices. */
122 for (int r = 0; r < 3; r++) {
123 for (int i = 0; i < F.rows(); i++)
124 cs.push_back(F(i, 1));
125 for (int i = 0; i < F.rows(); i++)
126 cs.push_back(F(i, 0));
127 for (int i = 0; i < F.rows(); i++)
128 cs.push_back(F(i, 2));
129 for (int i = 0; i < F.rows(); i++)
130 cs.push_back(F(i, 0));
131 }
132
133 /* Values. */
134 for (int i = 0; i < F.rows(); i++)
135 vs.push_back(eperp13(i, 0));
136 for (int i = 0; i < F.rows(); i++)
137 vs.push_back(-eperp13(i, 0));
138 for (int i = 0; i < F.rows(); i++)
139 vs.push_back(eperp21(i, 0));
140 for (int i = 0; i < F.rows(); i++)
141 vs.push_back(-eperp21(i, 0));
142 for (int i = 0; i < F.rows(); i++)
143 vs.push_back(eperp13(i, 1));
144 for (int i = 0; i < F.rows(); i++)
145 vs.push_back(-eperp13(i, 1));
146 for (int i = 0; i < F.rows(); i++)
147 vs.push_back(eperp21(i, 1));
148 for (int i = 0; i < F.rows(); i++)
149 vs.push_back(-eperp21(i, 1));
150 for (int i = 0; i < F.rows(); i++)
151 vs.push_back(eperp13(i, 2));
152 for (int i = 0; i < F.rows(); i++)
153 vs.push_back(-eperp13(i, 2));
154 for (int i = 0; i < F.rows(); i++)
155 vs.push_back(eperp21(i, 2));
156 for (int i = 0; i < F.rows(); i++)
157 vs.push_back(-eperp21(i, 2));
158
159 /* Create sparse gradient operator matrix.. */
160 G.resize(3 * F.rows(), V.rows());
161 std::vector<Eigen::Triplet<typename DerivedV::Scalar>> triplets;
162 for (int i = 0; i < (int)vs.size(); ++i) {
163 triplets.push_back(Eigen::Triplet<typename DerivedV::Scalar>(rs[i], cs[i], vs[i]));
164 }
165 G.setFromTriplets(triplets.begin(), triplets.end());
166}
167
168/* Computes the polar decomposition (R,T) of a matrix A using SVD singular
169 * value decomposition
170 *
171 * Inputs:
172 * A 3 by 3 matrix to be decomposed
173 * Outputs:
174 * R 3 by 3 rotation matrix part of decomposition (**always rotataion**)
175 * T 3 by 3 stretch matrix part of decomposition
176 * U 3 by 3 left-singular vectors
177 * S 3 by 1 singular values
178 * V 3 by 3 right-singular vectors
179 */
180template<typename DerivedA,
181 typename DerivedR,
182 typename DerivedT,
183 typename DerivedU,
184 typename DerivedS,
185 typename DerivedV>
186static inline void polar_svd(const Eigen::PlainObjectBase<DerivedA> &A,
187 Eigen::PlainObjectBase<DerivedR> &R,
188 Eigen::PlainObjectBase<DerivedT> &T,
189 Eigen::PlainObjectBase<DerivedU> &U,
190 Eigen::PlainObjectBase<DerivedS> &S,
191 Eigen::PlainObjectBase<DerivedV> &V)
192{
193 using namespace std;
194 Eigen::JacobiSVD<DerivedA> svd;
195 svd.compute(A, Eigen::ComputeFullU | Eigen::ComputeFullV);
196 U = svd.matrixU();
197 V = svd.matrixV();
198 S = svd.singularValues();
199 R = U * V.transpose();
200 const auto &SVT = S.asDiagonal() * V.adjoint();
201 /* Check for reflection. */
202 if (R.determinant() < 0) {
203 /* Annoyingly the .eval() is necessary. */
204 auto W = V.eval();
205 W.col(V.cols() - 1) *= -1.;
206 R = U * W.transpose();
207 T = W * SVT;
208 }
209 else {
210 T = V * SVT;
211 }
212}
213
214static inline void compute_surface_gradient_matrix(const Eigen::MatrixXd &V,
215 const Eigen::MatrixXi &F,
216 const Eigen::MatrixXd &F1,
217 const Eigen::MatrixXd &F2,
218 Eigen::SparseMatrix<double> &D1,
219 Eigen::SparseMatrix<double> &D2)
220{
221 Eigen::SparseMatrix<double> G;
222 grad(V, F, G);
223 Eigen::SparseMatrix<double> Dx = G.block(0, 0, F.rows(), V.rows());
224 Eigen::SparseMatrix<double> Dy = G.block(F.rows(), 0, F.rows(), V.rows());
225 Eigen::SparseMatrix<double> Dz = G.block(2 * F.rows(), 0, F.rows(), V.rows());
226
227 D1 = F1.col(0).asDiagonal() * Dx + F1.col(1).asDiagonal() * Dy + F1.col(2).asDiagonal() * Dz;
228 D2 = F2.col(0).asDiagonal() * Dx + F2.col(1).asDiagonal() * Dy + F2.col(2).asDiagonal() * Dz;
229}
230
231static inline void compute_weighted_jacobians(SLIMData &s, const Eigen::MatrixXd &uv)
232{
233 BLI_assert(s.valid);
234
235 /* Ji=[D1*u,D2*u,D1*v,D2*v] */
236 s.Ji.col(0) = s.Dx * uv.col(0);
237 s.Ji.col(1) = s.Dy * uv.col(0);
238 s.Ji.col(2) = s.Dx * uv.col(1);
239 s.Ji.col(3) = s.Dy * uv.col(1);
240
241 /* Add weights. */
242 Eigen::VectorXd weights = s.weightPerFaceMap.cast<double>();
243 s.Ji.col(0) = weights.cwiseProduct(s.Ji.col(0));
244 s.Ji.col(1) = weights.cwiseProduct(s.Ji.col(1));
245 s.Ji.col(2) = weights.cwiseProduct(s.Ji.col(2));
246 s.Ji.col(3) = weights.cwiseProduct(s.Ji.col(3));
247}
248
249static inline void compute_unweighted_jacobians(SLIMData &s, const Eigen::MatrixXd &uv)
250{
251 BLI_assert(s.valid);
252
253 /* Ji=[D1*u,D2*u,D1*v,D2*v] */
254 s.Ji.col(0) = s.Dx * uv.col(0);
255 s.Ji.col(1) = s.Dy * uv.col(0);
256 s.Ji.col(2) = s.Dx * uv.col(1);
257 s.Ji.col(3) = s.Dy * uv.col(1);
258}
259
260static inline void compute_jacobians(SLIMData &s, const Eigen::MatrixXd &uv)
261{
262 BLI_assert(s.valid);
263
264 if (s.withWeightedParameterization) {
265 compute_weighted_jacobians(s, uv);
266 }
267 else {
268 compute_unweighted_jacobians(s, uv);
269 }
270}
271
272static inline void update_weights_and_closest_rotations(SLIMData &s, Eigen::MatrixXd &uv)
273{
274 BLI_assert(s.valid);
275 compute_jacobians(s, uv);
276
277 const double eps = 1e-8;
278 double exp_f = s.exp_factor;
279
280 for (int i = 0; i < s.Ji.rows(); ++i) {
281 typedef Eigen::Matrix<double, 2, 2> Mat2;
282 typedef Eigen::Matrix<double, 2, 1> Vec2;
283 Mat2 ji, ri, ti, ui, vi;
284 Vec2 sing;
285 Vec2 closest_sing_vec;
286 Mat2 mat_W;
287 Vec2 m_sing_new;
288 double s1, s2;
289
290 ji(0, 0) = s.Ji(i, 0);
291 ji(0, 1) = s.Ji(i, 1);
292 ji(1, 0) = s.Ji(i, 2);
293 ji(1, 1) = s.Ji(i, 3);
294
295 polar_svd(ji, ri, ti, ui, sing, vi);
296
297 s1 = sing(0);
298 s2 = sing(1);
299
300 /* Update Weights according to energy. */
301 switch (s.slim_energy) {
302 case SLIMData::ARAP: {
303 m_sing_new << 1, 1;
304 break;
305 }
306 case SLIMData::SYMMETRIC_DIRICHLET: {
307 double s1_g = 2 * (s1 - pow(s1, -3));
308 double s2_g = 2 * (s2 - pow(s2, -3));
309 m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1)));
310 break;
311 }
312 case SLIMData::LOG_ARAP: {
313 double s1_g = 2 * (log(s1) / s1);
314 double s2_g = 2 * (log(s2) / s2);
315 m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1)));
316 break;
317 }
318 case SLIMData::CONFORMAL: {
319 double s1_g = 1 / (2 * s2) - s2 / (2 * pow(s1, 2));
320 double s2_g = 1 / (2 * s1) - s1 / (2 * pow(s2, 2));
321
322 double geo_avg = sqrt(s1 * s2);
323 double s1_min = geo_avg;
324 double s2_min = geo_avg;
325
326 m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min)));
327
328 /* Change local step. */
329 closest_sing_vec << s1_min, s2_min;
330 ri = ui * closest_sing_vec.asDiagonal() * vi.transpose();
331 break;
332 }
333 case SLIMData::EXP_CONFORMAL: {
334 double s1_g = 2 * (s1 - pow(s1, -3));
335 double s2_g = 2 * (s2 - pow(s2, -3));
336
337 double in_exp = exp_f * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2));
338 double exp_thing = exp(in_exp);
339
340 s1_g *= exp_thing * exp_f;
341 s2_g *= exp_thing * exp_f;
342
343 m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1)));
344 break;
345 }
346 case SLIMData::EXP_SYMMETRIC_DIRICHLET: {
347 double s1_g = 2 * (s1 - pow(s1, -3));
348 double s2_g = 2 * (s2 - pow(s2, -3));
349
350 double in_exp = exp_f * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2));
351 double exp_thing = exp(in_exp);
352
353 s1_g *= exp_thing * exp_f;
354 s2_g *= exp_thing * exp_f;
355
356 m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1)));
357 break;
358 }
359 }
360
361 if (std::abs(s1 - 1) < eps)
362 m_sing_new(0) = 1;
363 if (std::abs(s2 - 1) < eps)
364 m_sing_new(1) = 1;
365 mat_W = ui * m_sing_new.asDiagonal() * ui.transpose();
366
367 s.W_11(i) = mat_W(0, 0);
368 s.W_12(i) = mat_W(0, 1);
369 s.W_21(i) = mat_W(1, 0);
370 s.W_22(i) = mat_W(1, 1);
371
372 /* 2) Update local step (doesn't have to be a rotation, for instance in case of conformal
373 * energy). */
374 s.Ri(i, 0) = ri(0, 0);
375 s.Ri(i, 1) = ri(1, 0);
376 s.Ri(i, 2) = ri(0, 1);
377 s.Ri(i, 3) = ri(1, 1);
378 }
379}
380
381template<typename DerivedV, typename DerivedF>
382static inline void local_basis(const Eigen::PlainObjectBase<DerivedV> &V,
383 const Eigen::PlainObjectBase<DerivedF> &F,
384 Eigen::PlainObjectBase<DerivedV> &B1,
385 Eigen::PlainObjectBase<DerivedV> &B2,
386 Eigen::PlainObjectBase<DerivedV> &B3)
387{
388 using namespace Eigen;
389 using namespace std;
390 B1.resize(F.rows(), 3);
391 B2.resize(F.rows(), 3);
392 B3.resize(F.rows(), 3);
393
394 for (unsigned i = 0; i < F.rows(); ++i) {
395 Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v1 =
396 (V.row(F(i, 1)) - V.row(F(i, 0))).normalized();
397 Eigen::Matrix<typename DerivedV::Scalar, 1, 3> t = V.row(F(i, 2)) - V.row(F(i, 0));
398 Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v3 = v1.cross(t).normalized();
399 Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v2 = v1.cross(v3).normalized();
400
401 B1.row(i) = v1;
402 B2.row(i) = -v2;
403 B3.row(i) = v3;
404 }
405}
406
407static inline void pre_calc(SLIMData &s)
408{
409 BLI_assert(s.valid);
410 if (!s.has_pre_calc) {
411 s.v_n = s.v_num;
412 s.f_n = s.f_num;
413
414 s.dim = 2;
415 Eigen::MatrixXd F1, F2, F3;
416 local_basis(s.V, s.F, F1, F2, F3);
417 compute_surface_gradient_matrix(s.V, s.F, F1, F2, s.Dx, s.Dy);
418
419 s.W_11.resize(s.f_n);
420 s.W_12.resize(s.f_n);
421 s.W_21.resize(s.f_n);
422 s.W_22.resize(s.f_n);
423
424 s.Dx.makeCompressed();
425 s.Dy.makeCompressed();
426 s.Dz.makeCompressed();
427 s.Ri.resize(s.f_n, s.dim * s.dim);
428 s.Ji.resize(s.f_n, s.dim * s.dim);
429 s.rhs.resize(s.dim * s.v_num);
430
431 /* Flattened weight matrix. */
432 s.WGL_M.resize(s.dim * s.dim * s.f_n);
433 for (int i = 0; i < s.dim * s.dim; i++)
434 for (int j = 0; j < s.f_n; j++)
435 s.WGL_M(i * s.f_n + j) = s.M(j);
436
437 s.first_solve = true;
438 s.has_pre_calc = true;
439 }
440}
441
442static inline void buildA(SLIMData &s, Eigen::SparseMatrix<double> &A)
443{
444 BLI_assert(s.valid);
445 /* Formula (35) in paper. */
446 std::vector<Eigen::Triplet<double>> IJV;
447
448 IJV.reserve(4 * (s.Dx.outerSize() + s.Dy.outerSize()));
449
450 /* A = [W11*Dx, W12*Dx;
451 * W11*Dy, W12*Dy;
452 * W21*Dx, W22*Dx;
453 * W21*Dy, W22*Dy]; */
454 for (int k = 0; k < s.Dx.outerSize(); ++k) {
455 for (Eigen::SparseMatrix<double>::InnerIterator it(s.Dx, k); it; ++it) {
456 int dx_r = it.row();
457 int dx_c = it.col();
458 double val = it.value();
459 double weight = s.weightPerFaceMap(dx_r);
460
461 IJV.push_back(Eigen::Triplet<double>(dx_r, dx_c, weight * val * s.W_11(dx_r)));
462 IJV.push_back(Eigen::Triplet<double>(dx_r, s.v_n + dx_c, weight * val * s.W_12(dx_r)));
463
464 IJV.push_back(Eigen::Triplet<double>(2 * s.f_n + dx_r, dx_c, weight * val * s.W_21(dx_r)));
465 IJV.push_back(
466 Eigen::Triplet<double>(2 * s.f_n + dx_r, s.v_n + dx_c, weight * val * s.W_22(dx_r)));
467 }
468 }
469
470 for (int k = 0; k < s.Dy.outerSize(); ++k) {
471 for (Eigen::SparseMatrix<double>::InnerIterator it(s.Dy, k); it; ++it) {
472 int dy_r = it.row();
473 int dy_c = it.col();
474 double val = it.value();
475 double weight = s.weightPerFaceMap(dy_r);
476
477 IJV.push_back(Eigen::Triplet<double>(s.f_n + dy_r, dy_c, weight * val * s.W_11(dy_r)));
478 IJV.push_back(
479 Eigen::Triplet<double>(s.f_n + dy_r, s.v_n + dy_c, weight * val * s.W_12(dy_r)));
480
481 IJV.push_back(Eigen::Triplet<double>(3 * s.f_n + dy_r, dy_c, weight * val * s.W_21(dy_r)));
482 IJV.push_back(
483 Eigen::Triplet<double>(3 * s.f_n + dy_r, s.v_n + dy_c, weight * val * s.W_22(dy_r)));
484 }
485 }
486
487 A.setFromTriplets(IJV.begin(), IJV.end());
488}
489
490static inline void buildRhs(SLIMData &s, const Eigen::SparseMatrix<double> &At)
491{
492 BLI_assert(s.valid);
493
494 Eigen::VectorXd f_rhs(s.dim * s.dim * s.f_n);
495 f_rhs.setZero();
496
497 /* b = [W11*R11 + W12*R21; (formula (36))
498 * W11*R12 + W12*R22;
499 * W21*R11 + W22*R21;
500 * W21*R12 + W22*R22]; */
501 for (int i = 0; i < s.f_n; i++) {
502 f_rhs(i + 0 * s.f_n) = s.W_11(i) * s.Ri(i, 0) + s.W_12(i) * s.Ri(i, 1);
503 f_rhs(i + 1 * s.f_n) = s.W_11(i) * s.Ri(i, 2) + s.W_12(i) * s.Ri(i, 3);
504 f_rhs(i + 2 * s.f_n) = s.W_21(i) * s.Ri(i, 0) + s.W_22(i) * s.Ri(i, 1);
505 f_rhs(i + 3 * s.f_n) = s.W_21(i) * s.Ri(i, 2) + s.W_22(i) * s.Ri(i, 3);
506 }
507
508 Eigen::VectorXd uv_flat(s.dim * s.v_n);
509 for (int i = 0; i < s.dim; i++)
510 for (int j = 0; j < s.v_n; j++)
511 uv_flat(s.v_n * i + j) = s.V_o(j, i);
512
513 s.rhs = (At * s.WGL_M.asDiagonal() * f_rhs + s.proximal_p * uv_flat);
514}
515
516static inline void add_soft_constraints(SLIMData &s, Eigen::SparseMatrix<double> &L)
517{
518 BLI_assert(s.valid);
519 int v_n = s.v_num;
520 for (int d = 0; d < s.dim; d++) {
521 for (int i = 0; i < s.b.rows(); i++) {
522 int v_idx = s.b(i);
523 s.rhs(d * v_n + v_idx) += s.soft_const_p * s.bc(i, d); /* Right hand side. */
524 L.coeffRef(d * v_n + v_idx, d * v_n + v_idx) += s.soft_const_p; /* Diagonal of matrix. */
525 }
526 }
527}
528
529static inline void build_linear_system(SLIMData &s, Eigen::SparseMatrix<double> &L)
530{
531 BLI_assert(s.valid);
532 /* Formula (35) in paper. */
533 Eigen::SparseMatrix<double> A(s.dim * s.dim * s.f_n, s.dim * s.v_n);
534 buildA(s, A);
535
536 Eigen::SparseMatrix<double> At = A.transpose();
537 At.makeCompressed();
538
539 Eigen::SparseMatrix<double> id_m(At.rows(), At.rows());
540 id_m.setIdentity();
541
542 /* Add proximal penalty. */
543 L = At * s.WGL_M.asDiagonal() * A + s.proximal_p * id_m; /* Add also a proximal term. */
544 L.makeCompressed();
545
546 buildRhs(s, At);
547 Eigen::SparseMatrix<double> OldL = L;
548 add_soft_constraints(s, L);
549 L.makeCompressed();
550}
551
552static inline double compute_energy_with_jacobians(SLIMData &s,
553 const Eigen::MatrixXd &Ji,
554 Eigen::VectorXd &areas,
555 Eigen::VectorXd &singularValues,
556 bool gatherSingularValues)
557{
558 BLI_assert(s.valid);
559 double energy = 0;
560
561 Eigen::Matrix<double, 2, 2> ji;
562 for (int i = 0; i < s.f_n; i++) {
563 ji(0, 0) = Ji(i, 0);
564 ji(0, 1) = Ji(i, 1);
565 ji(1, 0) = Ji(i, 2);
566 ji(1, 1) = Ji(i, 3);
567
568 typedef Eigen::Matrix<double, 2, 2> Mat2;
569 typedef Eigen::Matrix<double, 2, 1> Vec2;
570 Mat2 ri, ti, ui, vi;
571 Vec2 sing;
572 polar_svd(ji, ri, ti, ui, sing, vi);
573 double s1 = sing(0);
574 double s2 = sing(1);
575
576 switch (s.slim_energy) {
577 case SLIMData::ARAP: {
578 energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2));
579 break;
580 }
581 case SLIMData::SYMMETRIC_DIRICHLET: {
582 energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2));
583
584 if (gatherSingularValues) {
585 singularValues(i) = s1;
586 singularValues(i + s.F.rows()) = s2;
587 }
588 break;
589 }
590 case SLIMData::EXP_SYMMETRIC_DIRICHLET: {
591 energy += areas(i) *
592 exp(s.exp_factor * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2)));
593 break;
594 }
595 case SLIMData::LOG_ARAP: {
596 energy += areas(i) * (pow(log(s1), 2) + pow(log(s2), 2));
597 break;
598 }
599 case SLIMData::CONFORMAL: {
600 energy += areas(i) * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2));
601 break;
602 }
603 case SLIMData::EXP_CONFORMAL: {
604 energy += areas(i) * exp(s.exp_factor * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2)));
605 break;
606 }
607 }
608 }
609
610 return energy;
611}
612
613static inline double compute_soft_const_energy(SLIMData &s, Eigen::MatrixXd &V_o)
614{
615 BLI_assert(s.valid);
616 double e = 0;
617 for (int i = 0; i < s.b.rows(); i++) {
618 e += s.soft_const_p * (s.bc.row(i) - V_o.row(s.b(i))).squaredNorm();
619 }
620 return e;
621}
622
623static inline double compute_energy(SLIMData &s,
624 Eigen::MatrixXd &V_new,
625 Eigen::VectorXd &singularValues,
626 bool gatherSingularValues)
627{
628 BLI_assert(s.valid);
629 compute_jacobians(s, V_new);
630 return compute_energy_with_jacobians(s, s.Ji, s.M, singularValues, gatherSingularValues) +
631 compute_soft_const_energy(s, V_new);
632}
633
634static inline double compute_energy(SLIMData &s, Eigen::MatrixXd &V_new)
635{
636 BLI_assert(s.valid);
637 Eigen::VectorXd temp;
638 return compute_energy(s, V_new, temp, false);
639}
640
641static inline double compute_energy(SLIMData &s,
642 Eigen::MatrixXd &V_new,
643 Eigen::VectorXd &singularValues)
644{
645 BLI_assert(s.valid);
646 return compute_energy(s, V_new, singularValues, true);
647}
648
649void slim_precompute(Eigen::MatrixXd &V,
650 Eigen::MatrixXi &F,
651 Eigen::MatrixXd &V_init,
652 SLIMData &data,
653 SLIMData::SLIM_ENERGY slim_energy,
654 Eigen::VectorXi &b,
655 Eigen::MatrixXd &bc,
656 double soft_p)
657{
658 BLI_assert(data.valid);
659 data.V = V;
660 data.F = F;
661 data.V_o = V_init;
662
663 data.v_num = V.rows();
664 data.f_num = F.rows();
665
666 data.slim_energy = slim_energy;
667
668 data.b = b;
669 data.bc = bc;
670 data.soft_const_p = soft_p;
671
672 data.proximal_p = 0.0001;
673
674 doublearea(V, F, data.M);
675 data.M /= 2.;
676 data.mesh_area = data.M.sum();
677
678 data.mesh_improvement_3d = false; /* Whether to use a jacobian derived from a real mesh or an
679 * abstract regular mesh (used for mesh improvement). */
680 data.exp_factor =
681 1.0; /* Param used only for exponential energies (e.g exponential symmetric dirichlet). */
682
683 assert(F.cols() == 3);
684
685 pre_calc(data);
686
687 data.energy = compute_energy(data, data.V_o) / data.mesh_area;
688}
689
690inline double computeGlobalScaleInvarianceFactor(Eigen::VectorXd &singularValues,
691 Eigen::VectorXd &areas)
692{
693 int nFaces = singularValues.rows() / 2;
694
695 Eigen::VectorXd areasChained(2 * nFaces);
696 areasChained << areas, areas;
697
698 /* Per face energy for face i with singvals si1 and si2 and area ai when scaling geometry by x is
699 *
700 * ai*(si1*x)^2 + ai*(si2*x)^2 + ai/(si1*x)^2 + ai/(si2*x)^2)
701 *
702 * The combined Energy of all faces is therefore
703 * (s1 and s2 are the sums over all ai*(si1^2) and ai*(si2^2) respectively. t1 and t2
704 * are the sums over all ai/(si1^2) and ai/(si2^2) respectively)
705 *
706 * s1*(x^2) + s2*(x^2) + t1/(x^2) + t2/(x^2)
707 *
708 * with a = (s1 + s2) and b = (t1 + t2) we get
709 *
710 * ax^2 + b/x^2
711 *
712 * it's derivative is
713 *
714 * 2ax - 2b/(x^3)
715 *
716 * and when we set it zero we get
717 *
718 * x^4 = b/a => x = sqrt(sqrt(b/a))
719 */
720
721 Eigen::VectorXd squaredSingularValues = singularValues.cwiseProduct(singularValues);
722 Eigen::VectorXd inverseSquaredSingularValues =
723 singularValues.cwiseProduct(singularValues).cwiseInverse();
724
725 Eigen::VectorXd weightedSquaredSingularValues = squaredSingularValues.cwiseProduct(areasChained);
726 Eigen::VectorXd weightedInverseSquaredSingularValues = inverseSquaredSingularValues.cwiseProduct(
727 areasChained);
728
729 double s1 = weightedSquaredSingularValues.head(nFaces).sum();
730 double s2 = weightedSquaredSingularValues.tail(nFaces).sum();
731
732 double t1 = weightedInverseSquaredSingularValues.head(nFaces).sum();
733 double t2 = weightedInverseSquaredSingularValues.tail(nFaces).sum();
734
735 double a = s1 + s2;
736 double b = t1 + t2;
737
738 double x = sqrt(sqrt(b / a));
739
740 return 1 / x;
741}
742
743static inline void solve_weighted_arap(SLIMData &s, Eigen::MatrixXd &uv)
744{
745 BLI_assert(s.valid);
746 using namespace Eigen;
747
748 Eigen::SparseMatrix<double> L;
749 build_linear_system(s, L);
750
751 /* Solve. */
752 Eigen::VectorXd Uc;
753 SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
754 Uc = solver.compute(L).solve(s.rhs);
755
756 for (int i = 0; i < Uc.size(); i++)
757 if (!std::isfinite(Uc(i)))
758 throw SlimFailedException();
759
760 for (int i = 0; i < s.dim; i++)
761 uv.col(i) = Uc.block(i * s.v_n, 0, s.v_n, 1);
762}
763
764Eigen::MatrixXd slim_solve(SLIMData &data, int iter_num)
765{
766 BLI_assert(data.valid);
767 Eigen::VectorXd singularValues;
768 bool are_pins_present = data.b.rows() > 0;
769
770 if (are_pins_present) {
771 singularValues.resize(data.F.rows() * 2);
772 data.energy = compute_energy(data, data.V_o, singularValues) / data.mesh_area;
773 }
774
775 for (int i = 0; i < iter_num; i++) {
776 Eigen::MatrixXd dest_res;
777 dest_res = data.V_o;
778
779 /* Solve Weighted Proxy. */
780 update_weights_and_closest_rotations(data, dest_res);
781 solve_weighted_arap(data, dest_res);
782
783 std::function<double(Eigen::MatrixXd &)> compute_energy_func = [&](Eigen::MatrixXd &aaa) {
784 return are_pins_present ? compute_energy(data, aaa, singularValues) :
785 compute_energy(data, aaa);
786 };
787
788 data.energy = flip_avoiding_line_search(data.F,
789 data.V_o,
790 dest_res,
791 compute_energy_func,
792 data.energy * data.mesh_area) /
793 data.mesh_area;
794
795 if (are_pins_present) {
796 data.globalScaleInvarianceFactor = computeGlobalScaleInvarianceFactor(singularValues,
797 data.M);
798 data.Dx /= data.globalScaleInvarianceFactor;
799 data.Dy /= data.globalScaleInvarianceFactor;
800 data.energy = compute_energy(data, data.V_o, singularValues) / data.mesh_area;
801 }
802 }
803
804 return data.V_o;
805}
806
807} // namespace slim
BLI_assert.h
BLI_assert
#define BLI_assert(a)
Definition BLI_assert.h:50
sqrt
sqrt(x)+1/max(0
x
x
Definition BLI_expr_pylike_eval_test.cc:345
BLI_math_base.h
M_PI
#define M_PI
Definition BLI_math_base.h:58
double
typedef double(DMatrix)[4][4]
v2
ATTR_WARN_UNUSED_RESULT const BMVert * v2
Definition bmesh_query_inline.hh:35
e
ATTR_WARN_UNUSED_RESULT const BMVert const BMEdge * e
Definition bmesh_query_inline.hh:36
A
#define A
U
unsigned int U
Definition btGjkEpa3.h:78
normalized
SIMD_FORCE_INLINE btVector3 normalized() const
Return a normalized version of this vector.
b
local_group_size(16, 16) .push_constant(Type b
Definition compositor_morphological_distance_info.hh:16
pow
pow(value.r - subtrahend, 2.0)") .do_static_compilation(true)
int
draw_view push_constant(Type::INT, "radiance_src") .push_constant(Type capture_info_buf storage_buf(1, Qualifier::READ, "ObjectBounds", "bounds_buf[]") .push_constant(Type draw_view int
Definition eevee_lightprobe_volume_info.hh:143
flip_avoiding_line_search.h
exp
ccl_device_inline float3 exp(float3 v)
Definition math_float3.h:365
log
ccl_device_inline float3 log(float3 v)
Definition math_float3.h:370
F
#define F
Definition mball_tessellate.cc:275
R
#define R
Definition mball_tessellate.cc:271
L
#define L
Definition mball_tessellate.cc:270
G
#define G(x, y, z)
slim::compute_soft_const_energy
static double compute_soft_const_energy(SLIMData &s, Eigen::MatrixXd &V_o)
Definition slim.cpp:613
slim::slim_precompute
void slim_precompute(Eigen::MatrixXd &V, Eigen::MatrixXi &F, Eigen::MatrixXd &V_init, SLIMData &data, SLIMData::SLIM_ENERGY slim_energy, Eigen::VectorXi &b, Eigen::MatrixXd &bc, double soft_p)
Definition slim.cpp:649
slim::compute_jacobians
static void compute_jacobians(SLIMData &s, const Eigen::MatrixXd &uv)
Definition slim.cpp:260
slim::update_weights_and_closest_rotations
static void update_weights_and_closest_rotations(SLIMData &s, Eigen::MatrixXd &uv)
Definition slim.cpp:272
slim::polar_svd
static void polar_svd(const Eigen::PlainObjectBase< DerivedA > &A, Eigen::PlainObjectBase< DerivedR > &R, Eigen::PlainObjectBase< DerivedT > &T, Eigen::PlainObjectBase< DerivedU > &U, Eigen::PlainObjectBase< DerivedS > &S, Eigen::PlainObjectBase< DerivedV > &V)
Definition slim.cpp:186
slim::doublearea
void doublearea(const Eigen::PlainObjectBase< DerivedV > &V, const Eigen::PlainObjectBase< DerivedF > &F, Eigen::PlainObjectBase< DeriveddblA > &dblA)
Definition doublearea.cpp:116
slim::buildRhs
static void buildRhs(SLIMData &s, const Eigen::SparseMatrix< double > &At)
Definition slim.cpp:490
slim::compute_surface_gradient_matrix
static void compute_surface_gradient_matrix(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXd &F1, const Eigen::MatrixXd &F2, Eigen::SparseMatrix< double > &D1, Eigen::SparseMatrix< double > &D2)
Definition slim.cpp:214
slim::compute_energy_with_jacobians
static double compute_energy_with_jacobians(SLIMData &s, const Eigen::MatrixXd &Ji, Eigen::VectorXd &areas, Eigen::VectorXd &singularValues, bool gatherSingularValues)
Definition slim.cpp:552
slim::build_linear_system
static void build_linear_system(SLIMData &s, Eigen::SparseMatrix< double > &L)
Definition slim.cpp:529
slim::buildA
static void buildA(SLIMData &s, Eigen::SparseMatrix< double > &A)
Definition slim.cpp:442
slim::pre_calc
static void pre_calc(SLIMData &s)
Definition slim.cpp:407
slim::computeGlobalScaleInvarianceFactor
double computeGlobalScaleInvarianceFactor(Eigen::VectorXd &singularValues, Eigen::VectorXd &areas)
Definition slim.cpp:690
slim::solve_weighted_arap
static void solve_weighted_arap(SLIMData &s, Eigen::MatrixXd &uv)
Definition slim.cpp:743
slim::local_basis
static void local_basis(const Eigen::PlainObjectBase< DerivedV > &V, const Eigen::PlainObjectBase< DerivedF > &F, Eigen::PlainObjectBase< DerivedV > &B1, Eigen::PlainObjectBase< DerivedV > &B2, Eigen::PlainObjectBase< DerivedV > &B3)
Definition slim.cpp:382
slim::compute_weighted_jacobians
static void compute_weighted_jacobians(SLIMData &s, const Eigen::MatrixXd &uv)
Definition slim.cpp:231
slim::add_soft_constraints
static void add_soft_constraints(SLIMData &s, Eigen::SparseMatrix< double > &L)
Definition slim.cpp:516
slim::slim_solve
Eigen::MatrixXd slim_solve(SLIMData &data, int iter_num)
Definition slim.cpp:764
slim::flip_avoiding_line_search
double flip_avoiding_line_search(const Eigen::MatrixXi F, Eigen::MatrixXd &cur_v, Eigen::MatrixXd &dst_v, std::function< double(Eigen::MatrixXd &)> energy, double cur_energy)
Definition flip_avoiding_line_search.cpp:155
slim::compute_energy
static double compute_energy(SLIMData &s, Eigen::MatrixXd &V_new, Eigen::VectorXd &singularValues, bool gatherSingularValues)
Definition slim.cpp:623
slim::compute_unweighted_jacobians
static void compute_unweighted_jacobians(SLIMData &s, const Eigen::MatrixXd &uv)
Definition slim.cpp:249
slim::grad
static void grad(const Eigen::PlainObjectBase< DerivedV > &V, const Eigen::PlainObjectBase< DerivedF > &F, Eigen::SparseMatrix< typename DerivedV::Scalar > &G, bool uniform=false)
Definition slim.cpp:47
eps
const btScalar eps
Definition poly34.cpp:11
slim::SLIMData::SYMMETRIC_DIRICHLET
@ SYMMETRIC_DIRICHLET
Definition slim.h:34
slim::SLIMData::EXP_CONFORMAL
@ EXP_CONFORMAL
Definition slim.h:36
slim::SLIMData::EXP_SYMMETRIC_DIRICHLET
@ EXP_SYMMETRIC_DIRICHLET
Definition slim.h:37
Vec2
double Vec2[2]
Definition tracking_plane_tracker.cc:23
W
#define W
Definition uvedit_clipboard_graph_iso.cc:25
V
CCL_NAMESPACE_BEGIN struct Window V
Generated on Thu Feb 6 2025 07:36:39 for Blender by  doxygen 1.11.0