GCC Code Coverage Report

Directory: ./
File: lib/geogram/parameterization/mesh_ABF.cpp
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 0 365 0.0%
Functions: 0 25 0.0%
Branches: 0 300 0.0%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/parameterization/mesh_ABF.h>
41 #include <geogram/parameterization/mesh_LSCM.h>
42 #include <geogram/mesh/mesh.h>
43 #include <geogram/mesh/mesh_geometry.h>
44 #include <geogram/bibliography/bibliography.h>
45 #include <geogram/basic/memory.h>
46 #include <geogram/NL/nl.h>
47 #include <geogram/NL/nl_matrix.h>
48
49
50 namespace {
51 using namespace GEO;
52
53 class ABFPlusPlus {
54 public:
55 ✗ ABFPlusPlus(Mesh& mesh) :
56 ✗ mesh_(mesh),
57 ✗ epsilon_(1e-5),
58 ✗ newton_tolf_(1.0),
59 ✗ newton_tolx_(1.0),
60 ✗ max_newton_iter_(10),
61 ✗ positive_angle_ro_(1.2),
62 ✗ step_length_factor_(1.0) {
63 ✗ v_on_border_.assign(mesh.vertices.nb(),false);
64 ✗ v_to_c_.assign(mesh.vertices.nb(),NO_VERTEX);
65 ✗ next_c_around_v_.assign(mesh.facet_corners.nb(), NO_CORNER);
66 ✗ for(index_t c: mesh.facet_corners) {
67 index_t v = mesh.facet_corners.vertex(c);
68 ✗ if(mesh.facet_corners.adjacent_facet(c) == NO_FACET) {
69 ✗ v_on_border_[v] = true;
70 }
71 ✗ next_c_around_v_[c] = v_to_c_[v];
72 ✗ v_to_c_[v] = c;
73 }
74 ✗ if(!mesh_.facets.are_simplices()) {
75 ✗ c_to_f_.resize(mesh_.facet_corners.nb());
76 ✗ for(index_t f: mesh_.facets) {
77 ✗ for(index_t c: mesh_.facets.corners(f)) {
78 ✗ c_to_f_[c] = f;
79 }
80 }
81 }
82 ✗ verbose_ = false;
83 ✗ }
84
85 ✗ ~ABFPlusPlus() {
86 ✗ deallocate_variables();
87 ✗ }
88
89 void set_verbose(bool x) {
90 ✗ verbose_ = x;
91 }
92
93 ✗ bool parameterize() {
94 ✗ geo_cite("DBLP:journals/tog/ShefferLMB05");
95 ✗ allocate_variables();
96 ✗ compute_beta();
97 ✗ angle_.bind(mesh_.facet_corners.attributes(),"angle");
98 ✗ if(!solve_angles()) {
99 ✗ if(verbose_) {
100 ✗ Logger::err("ABF++") << "Did not converge." << std::endl ;
101 ✗ Logger::err("ABF++") << "Switching to LSCM" << std::endl ;
102 }
103 // Note: AnglesToUV with angles measured on the mesh
104 // (i.e. beta's) = LSCM !!!
105 ✗ for(index_t c: mesh_.facet_corners) {
106 ✗ angle_[c] = beta_[c];
107 }
108 }
109 ✗ deallocate_variables() ;
110 ✗ angle_.unbind();
111 ✗ return true ;
112 }
113
114 protected:
115 index_t c_to_f(index_t c) {
116 geo_debug_assert(c < mesh_.facet_corners.nb());
117 ✗ return mesh_.facets.are_simplices() ? (c/3) : c_to_f_[c];
118 }
119
120 index_t nb_interior_vertices(const Mesh& M) const {
121 index_t result=0;
122 ✗ for(index_t v: M.vertices) {
123 ✗ if(!v_on_border_[v]) {
124 ✗ ++result;
125 }
126 }
127 return result;
128 }
129
130 ✗ void allocate_variables() {
131 // ------- sizes & indexes ------
132 ✗ nf_ = mesh_.facets.nb();
133 ✗ nalpha_ = mesh_.facet_corners.nb();
134 ✗ nint_ = nb_interior_vertices(mesh_);
135 ✗ nlambda_ = nf_ + 2*nint_;
136 ✗ ntot_ = nalpha_ + nlambda_;
137
138 // ------- ABF variables --------
139 ✗ alpha_.resize(nalpha_) ;
140 ✗ lambda_.resize(nlambda_) ;
141 ✗ beta_.resize(nalpha_) ;
142 ✗ w_.resize(nalpha_) ;
143
144 // ------- Step vectors ---------
145 ✗ dalpha_.resize(nalpha_) ;
146 ✗ dlambda1_.resize(nf_) ;
147 ✗ dlambda2_.resize(2*nint_) ;
148
149 // ------- Gradients ------------
150 ✗ b1_.resize(nalpha_) ;
151 ✗ b2_.resize(nlambda_) ;
152
153 // ------- Jacobian -------------
154 ✗ nlSparseMatrixConstruct(
155 ✗ &J2_, NLuint(2*nint_), NLuint(nalpha_), NL_MATRIX_STORE_COLUMNS
156 );
157
158 // ------- ABF++ ----------------
159 ✗ nlSparseMatrixConstruct(
160 ✗ &J_star_, NLuint(2*nint_), NLuint(nf_), NL_MATRIX_STORE_COLUMNS
161 );
162 ✗ nlSparseMatrixConstruct(
163 ✗ &M_, NLuint(2*nint_), NLuint(2*nint_), NL_MATRIX_STORE_ROWS
164 );
165 ✗ }
166
167 ✗ void deallocate_variables() {
168 // ------- ABF variables --------
169 alpha_.clear() ;
170 lambda_.clear() ;
171 beta_.clear() ;
172 w_.clear() ;
173
174 // ------- Step vectors ---------
175 dalpha_.clear() ;
176 dlambda1_.clear() ;
177 dlambda2_.clear() ;
178
179 // ------- Gradients ------------
180 b1_.clear() ;
181 b2_.clear() ;
182 ✗ nlSparseMatrixDestroy(&J2_);
183
184 // ------- ABF++ ----------------
185 ✗ nlSparseMatrixDestroy(&J_star_);
186 ✗ nlSparseMatrixDestroy(&M_);
187 ✗ }
188
189
190 ✗ void compute_beta() {
191
192 ✗ for(index_t v: mesh_.vertices) {
193 // Compute sum_angles
194 double sum_angle = 0.0 ;
195 {
196 ✗ index_t c = v_to_c_[v];
197 do {
198 ✗ double angle = corner_angle(c) ;
199 ✗ sum_angle += angle ;
200 ✗ c = next_c_around_v_[c];
201 ✗ } while(c != NO_CORNER) ;
202 }
203
204 double ratio = 1.0 ;
205
206 ✗ if(!v_on_border_[v]) {
207 ✗ ratio = 2.0 * M_PI / sum_angle ;
208 }
209
210 {
211 index_t c = v_to_c_[v];
212 do {
213 ✗ beta_[c] = corner_angle(c) * ratio ;
214 ✗ beta_[c] = std::max(beta_[c], 3.0 * M_PI / 180.0);
215 ✗ beta_[c] = std::min(beta_[c], 175.0 * M_PI / 180.0);
216 ✗ c = next_c_around_v_[c];
217 ✗ } while(c != NO_CORNER) ;
218 }
219 }
220 ✗ }
221
222 ✗ double corner_angle(index_t c) {
223 index_t f = c_to_f(c);
224 ✗ index_t c_prev = mesh_.facets.prev_corner_around_facet(f,c);
225 index_t c_next = mesh_.facets.next_corner_around_facet(f,c);
226
227 const vec3& p1 = mesh_.facet_corners.point(c);
228 const vec3& p2 = mesh_.facet_corners.point(c_next);
229 const vec3& p3 = mesh_.facet_corners.point(c_prev);
230 ✗ double result = Geom::angle(p2-p1,p3-p1);
231 result = std::max(result, 2.0 * M_PI / 360.0) ;
232 ✗ return result ;
233 }
234
235
236 ✗ bool solve_angles() {
237
238 ✗ if(!nlInitExtension("SUPERLU")) {
239 ✗ Logger::warn("ABF++")
240 << "Could not initialize SuperLU extension"
241 << std::endl;
242 ✗ return false;
243 }
244
245 // Initial values
246 ✗ lambda_.assign(lambda_.size(),0.0);
247 ✗ for(index_t i=0; i<nalpha_; i++) {
248 ✗ alpha_[i] = beta_[i] ;
249 ✗ w_[i] = 1.0 / (beta_[i] * beta_[i]) ;
250 }
251
252 ✗ for(index_t k=0; k<max_newton_iter_; k++) {
253
254 // Compute Jacobian
255 ✗ nlSparseMatrixZero(&J2_);
256 ✗ add_JC2();
257 ✗ add_JC3();
258
259 // Compute rhs
260 ✗ b1_.assign(b1_.size(), 0.0);
261 ✗ b2_.assign(b2_.size(), 0.0);
262 ✗ sub_grad_F();
263 ✗ sub_grad_C1();
264 ✗ sub_grad_C2();
265 ✗ sub_grad_C3();
266
267 ✗ double errf_k = errf() ;
268
269
270 ✗ if(verbose_) {
271 ✗ Logger::out("ABF++")
272 << "iter= " << k << " errf= " << errf_k
273 << std::endl ;
274 }
275
276
277 ✗ if(Numeric::is_nan(errf_k) || errf_k > 1e18) {
278 ✗ if(verbose_) {
279 ✗ Logger::err("ABF++") << "errf=" << errf_k << std::endl ;
280 }
281 ✗ return false ;
282 }
283
284
285 ✗ if(errf_k <= newton_tolf_) {
286 ✗ if(verbose_) {
287 ✗ Logger::out("ABF++") << "converged" << std::endl ;
288 }
289 ✗ return test_and_commit_solution() ;
290 }
291
292 ✗ solve_current_iteration() ;
293
294 double errx ;
295 ✗ double s = compute_step_length_and_update_weights() ;
296 ✗ s *= step_length_factor_ ;
297 ✗ errx = compute_errx_and_update_x(s) ;
298
299 // TODO: since errx is dependent on the size of the mesh,
300 // weight the threshold by the size of the mesh.
301
302 ✗ if(Numeric::is_nan(errx) || errx > 1e15) {
303 ✗ if(verbose_) {
304 ✗ Logger::err("ABF++") << "errx: " << errx << std::endl ;
305 }
306 ✗ return false ;
307 }
308
309 ✗ if(verbose_) {
310 ✗ Logger::out("ABF++") << "iter= " << k << " errx= " << errx
311 << std::endl ;
312 }
313
314 ✗ if(errx <= newton_tolx_) {
315 ✗ if(verbose_) {
316 ✗ Logger::out("ABF++") << "converged" << std::endl ;
317 }
318 ✗ return test_and_commit_solution();
319 }
320
321 }
322
323 ✗ if(verbose_) {
324 ✗ Logger::out("ABF++") << "ran out of Newton iters" << std::endl ;
325 }
326
327 ✗ return test_and_commit_solution() ;
328 }
329
330 ✗ bool test_and_commit_solution() {
331 ✗ for(index_t c: mesh_.facet_corners) {
332 ✗ if(Numeric::is_nan(alpha_[c])) {
333 ✗ if(verbose_) {
334 ✗ Logger::err("ABF++") << "solution has nan"
335 << std::endl;
336 }
337 ✗ return false;
338 }
339 ✗ if(alpha_[c] == 0.0) {
340 ✗ if(verbose_) {
341 ✗ Logger::err("ABF++") << "solution has null angle"
342 << std::endl;
343 }
344 ✗ return false;
345 }
346 }
347 ✗ if(verbose_) {
348 ✗ Logger::out("ABF++") << "solution OK"
349 << std::endl;
350 }
351 if(angle_.is_bound()) {
352 ✗ for(index_t c: mesh_.facet_corners) {
353 ✗ angle_[c] = alpha_[c];
354 }
355 }
356 return true;
357 }
358
359 ✗ void solve_current_iteration() {
360
361 ✗ Delta_inv_.resize(nalpha_) ;
362 ✗ for(index_t i=0; i<nalpha_; i++) {
363 ✗ Delta_inv_[i] = 1.0 / (2.0 * w_[i]) ;
364 }
365
366
367 // 1) Create the pieces of J.Delta^-1.Jt
368 // 1.1) Diagonal part: Delta*^-1
369 ✗ Delta_star_inv_.resize(nf_) ;
370 ✗ for(index_t f=0; f<nf_; ++f) {
371 double S = 0.0;
372 ✗ for(index_t c: mesh_.facets.corners(f)) {
373 ✗ S += Delta_inv_[c];
374 }
375 ✗ Delta_star_inv_[f] = 1.0 / S;
376 }
377
378 // 1.2) J* = J2.Delta^-1.J1^t
379 ✗ nlSparseMatrixZero(&J_star_);
380 ✗ for(index_t j=0; j<nalpha_; j++) {
381 ✗ const NLRowColumn& Cj = J2_.column[j];
382 ✗ for(NLuint ii=0; ii<Cj.size; ++ii) {
383 ✗ const NLCoeff& c = Cj.coeff[ii] ;
384 ✗ nlSparseMatrixAdd(
385 &J_star_,
386 ✗ NLuint(c.index), NLuint(c_to_f(j)),
387 ✗ c.value * Delta_inv_[j]
388 );
389 }
390 }
391 // Note: J** does not need to be built, it is directly added to M.
392
393 // 2) Right hand side: b1* and b2*
394
395 // 2.1) b1* = J1.Delta^-1.b1 - b2[1..nf]
396 ✗ b1_star_.resize(nf_);
397
398 ✗ for(index_t f=0; f<nf_; ++f) {
399 ✗ b1_star_[f] = -b2_[f];
400 ✗ for(index_t c: mesh_.facets.corners(f)) {
401 ✗ b1_star_[f] += Delta_inv_[c] * b1_[c];
402 }
403 }
404
405 // 2.2) b2* = J2.Delta^-1.b1 - b2[nf+1 .. nf+2.nint-1]
406 ✗ b2_star_.assign(2*nint_,0.0);
407 ✗ add_J_D_x(b2_star_, J2_, Delta_inv_, b1_);
408 ✗ for(index_t i=0; i<2*nint_; i++) {
409 ✗ b2_star_[i] -= b2_[nf_+i];
410 }
411
412
413 // 3) create final linear system
414
415 // 3.1) M = J*.Delta*^-1.J*^t - J**
416 // where J** = J2.Delta^-1.J2^t
417 ✗ nlSparseMatrixZero(&M_);
418 ✗ add_J_D_Jt(M_, J_star_, Delta_star_inv_) ;
419 ✗ sub_J_D_Jt(M_, J2_, Delta_inv_) ;
420
421 // 3.2) r = J*.Delta*^-1.b1* - b2*
422 ✗ r_.assign(2*nint_,0.0) ;
423 ✗ add_J_D_x(r_, J_star_, Delta_star_inv_, b1_star_);
424
425 geo_debug_assert(r_.size() == b2_star_.size());
426 ✗ for(index_t i=0; i<r_.size(); ++i) {
427 ✗ r_[i] -= b2_star_[i];
428 }
429
430 ✗ if(verbose_) {
431 ✗ Logger::out("ABF++") << "Solving linear system..." << std::endl;
432 }
433 NLMatrix Minv =
434 ✗ nlMatrixFactorize((NLMatrix)&M_, NL_PERM_SUPERLU_EXT);
435 ✗ nlMultMatrixVector(Minv, r_.data(), dlambda2_.data());
436 ✗ nlDeleteMatrix(Minv);
437 ✗ if(verbose_) {
438 ✗ Logger::out("ABF++") << "Solved" << std::endl;
439 }
440
441 // 4) compute dlambda1 and dalpha in function of dlambda2
442
443 // 4.1) dlambda1 = Delta*^-1 ( b1* - J*^t dlambda2 )
444 ✗ mult_transpose(J_star_, dlambda2_, dlambda1_) ;
445 ✗ for(index_t f=0; f<nf_; ++f) {
446 ✗ dlambda1_[f] =
447 ✗ Delta_star_inv_[f] * (b1_star_[f] - dlambda1_[f]) ;
448 }
449
450 // 4.2) Compute dalpha in function of dlambda:
451 // dalpha = Delta^-1( b1 - J^t.dlambda )
452 // = Delta^-1( b1 - (J1^t.dlambda1 + J2^t.dlambda2) )
453 ✗ mult_transpose(J2_, dlambda2_, dalpha_) ;
454
455 ✗ for(index_t f=0; f<nf_; ++f) {
456 ✗ for(index_t c: mesh_.facets.corners(f)) {
457 ✗ dalpha_[c] += dlambda1_[f];
458 }
459 }
460
461 ✗ for(index_t i=0; i<nalpha_; i++) {
462 ✗ dalpha_[i] = Delta_inv_[i] * (b1_[i] - dalpha_[i]) ;
463 }
464 ✗ }
465
466
467 ✗ double compute_errx_and_update_x(double s) {
468 double result = 0 ;
469
470 // alpha += s * dalpha
471 ✗ for(index_t i=0; i<nalpha_; i++) {
472 ✗ double dai = s * dalpha_[i];
473 ✗ alpha_[i] += dai ;
474 ✗ result += ::fabs(dai) ;
475 }
476
477 // lambda += s * dlambda
478 ✗ for(index_t i=0; i<nf_; i++) {
479 ✗ double dai = s * dlambda1_[i];
480 ✗ lambda_[i] += dai ;
481 ✗ result += ::fabs(dai) ;
482 }
483
484 ✗ for(index_t i=0; i<2*nint_; i++) {
485 ✗ double dai = s * dlambda2_[i];
486 ✗ lambda_[nf_+i] += dai ;
487 ✗ result += ::fabs(dai) ;
488 }
489 ✗ return result ;
490 }
491
492
493 // --------------------- Jacobian ----------------------------
494
495 ✗ void add_JC2() {
496 index_t i = 0 ;
497 ✗ for(index_t v: mesh_.vertices) {
498 ✗ if(v_on_border_[v]) {
499 ✗ continue ;
500 }
501 ✗ index_t c = v_to_c_[v];
502 do {
503 ✗ nlSparseMatrixAdd(&J2_, NLuint(i), NLuint(c), 1.0);
504 ✗ c = next_c_around_v_[c];
505 ✗ } while(c != NO_CORNER);
506 ✗ i++ ;
507 }
508 ✗ }
509
510 ✗ void add_JC3() {
511 ✗ index_t i = nint_ ;
512 ✗ for(index_t v: mesh_.vertices) {
513 ✗ if(v_on_border_[v]) {
514 ✗ continue ;
515 }
516 double prod_prev_sin ;
517 double prod_next_sin ;
518 ✗ compute_product_sin_angles(v, prod_prev_sin, prod_next_sin) ;
519 ✗ index_t c = v_to_c_[v];
520 do {
521 index_t f = c_to_f(c);
522 index_t next_c = mesh_.facets.next_corner_around_facet(f,c);
523 ✗ nlSparseMatrixAdd(
524 &J2_, NLuint(i), NLuint(next_c),
525 ✗ prod_next_sin * cos(alpha_[next_c])/sin(alpha_[next_c])
526 );
527 ✗ index_t prev_c = mesh_.facets.prev_corner_around_facet(f,c);
528 ✗ nlSparseMatrixAdd(
529 &J2_, NLuint(i), NLuint(prev_c),
530 ✗ -prod_prev_sin * cos(alpha_[prev_c])/sin(alpha_[prev_c])
531 );
532 ✗ c = next_c_around_v_[c];
533 ✗ } while(c != NO_CORNER);
534 ✗ i++ ;
535 }
536 ✗ }
537
538 // --------------------- Right-hand side ---------------------
539
540 ✗ void sub_grad_F() {
541 ✗ for(index_t i = 0; i < nalpha_; ++i) {
542 ✗ b1_[i] -= 2.0 * w_[i] * ( alpha_[i] - beta_[i] );
543 }
544 ✗ }
545
546 // For each facet: sum angles - PI * (nb_vertices(f)-2)
547 ✗ void sub_grad_C1() {
548 ✗ for(index_t f=0; f < nf_; ++f) {
549 ✗ for(index_t c: mesh_.facets.corners(f)) {
550 ✗ b1_[c] -= lambda_[f];
551 }
552 }
553 ✗ for(index_t f=0; f < nf_; ++f) {
554 ✗ b2_[f] += M_PI * double(mesh_.facets.nb_vertices(f)-2);
555 ✗ for(index_t c: mesh_.facets.corners(f)) {
556 ✗ b2_[f] -= alpha_[c];
557 }
558 }
559 ✗ }
560
561 ✗ void sub_grad_C2() {
562 ✗ index_t i = nf_ ;
563 ✗ for(index_t v: mesh_.vertices) {
564 ✗ if(v_on_border_[v]) {
565 ✗ continue ;
566 }
567 ✗ index_t c = v_to_c_[v];
568 do {
569 ✗ b2_[i] -= alpha_[c];
570 ✗ b1_[c] -= lambda_[i];
571 ✗ c = next_c_around_v_[c];
572 ✗ } while(c != NO_CORNER) ;
573 ✗ b2_[i] += 2.0 * M_PI ;
574 ✗ ++i;
575 }
576 ✗ }
577
578
579 // For each vertex: prod sin(next angle) - prod sin(prev angle)
580 ✗ void sub_grad_C3() {
581 ✗ index_t i = nf_ + nint_ ;
582 ✗ for(index_t v: mesh_.vertices) {
583 ✗ if(v_on_border_[v]) {
584 ✗ continue ;
585 }
586
587 double prod_prev_sin ;
588 double prod_next_sin ;
589 ✗ compute_product_sin_angles(v, prod_prev_sin, prod_next_sin) ;
590
591 ✗ b2_[i] -= prod_next_sin - prod_prev_sin ;
592
593 ✗ index_t c = v_to_c_[v];
594 do {
595 index_t f = c_to_f(c);
596 index_t next_c = mesh_.facets.next_corner_around_facet(f,c);
597 ✗ b1_[next_c] -=
598 ✗ lambda_[i] * prod_next_sin *
599 ✗ cos(alpha_[next_c]) / sin(alpha_[next_c]) ;
600
601 ✗ index_t prev_c = mesh_.facets.prev_corner_around_facet(f,c);
602 ✗ b1_[prev_c] +=
603 ✗ lambda_[i] * prod_prev_sin *
604 ✗ cos(alpha_[prev_c]) / sin(alpha_[prev_c]) ;
605
606 ✗ c = next_c_around_v_[c];
607 ✗ } while(c != NO_CORNER);
608 ✗ i++ ;
609 }
610 ✗ }
611
612 // -------------------------------------------------------
613
614 ✗ void compute_product_sin_angles(
615 index_t v, double& prod_prev_sin, double& prod_next_sin
616 ) {
617 ✗ prod_prev_sin = 1.0 ;
618 ✗ prod_next_sin = 1.0 ;
619 ✗ index_t c = v_to_c_[v];
620 do {
621 index_t f = c_to_f(c);
622 ✗ index_t prev_c = mesh_.facets.prev_corner_around_facet(f,c);
623 index_t next_c = mesh_.facets.next_corner_around_facet(f,c);
624 ✗ prod_prev_sin *= sin(alpha_[prev_c]);
625 ✗ prod_next_sin *= sin(alpha_[next_c]);
626 ✗ c = next_c_around_v_[c];
627 ✗ } while(c != NO_CORNER) ;
628 ✗ }
629
630 // --------------------- Convergence control -----------------
631
632 ✗ double compute_step_length_and_update_weights() {
633 ✗ double ratio = 1.0 ;
634 ✗ for(index_t i=0; i<nalpha_; i++) {
635 ✗ if(alpha_[i] + dalpha_[i] < 10.0 * epsilon_) {
636 ✗ double r1 = -.5 * (alpha_[i] - 10.0 * epsilon_)/dalpha_[i];
637 ✗ ratio = std::min(ratio, r1) ;
638 ✗ w_[i] *= positive_angle_ro_ ;
639 ✗ } else if(alpha_[i] + dalpha_[i] > M_PI - 10.0 * epsilon_) {
640 // double r1 =
641 // .5*(M_PI - alpha_[i]+10.0 * epsilon_)/dalpha_[i];
642 // ratio = ogf_min(ratio, r1) ;
643 // two previous lines commented-out, I'm unsure why, to be
644 // tested.
645 ✗ w_[i] *= positive_angle_ro_ ;
646 }
647 }
648 ✗ return ratio ;
649 }
650
651 ✗ double errf() const {
652 double result = 0 ;
653 ✗ for(index_t i=0; i<nalpha_; ++i) {
654 ✗ result += ::fabs(b1_[i]);
655 }
656 ✗ for(index_t i=0; i<nlambda_; ++i) {
657 ✗ result += ::fabs(b2_[i]);
658 }
659 ✗ return result ;
660 }
661
662
663 // -------------------- Matrix utilities ---------------------
664
665 ✗ static void mult_transpose(
666 const NLSparseMatrix& M,
667 const vector<double>& x,
668 vector<double>& y
669 ) {
670 geo_debug_assert(y.size() == M.n);
671 geo_debug_assert(x.size() == M.m);
672 ✗ if((M.storage & NL_MATRIX_STORE_COLUMNS) != 0) {
673 ✗ for(NLuint j=0; j<M.n; ++j) {
674 ✗ y[j] = 0.0;
675 ✗ const NLRowColumn& Cj = M.column[j];
676 ✗ for(NLuint ii=0; ii<Cj.size; ++ii) {
677 ✗ double a = Cj.coeff[ii].value;
678 ✗ index_t i = Cj.coeff[ii].index;
679 ✗ y[j] += a * x[i];
680 }
681 }
682 } else {
683 ✗ geo_assert((M.storage & NL_MATRIX_STORE_ROWS) != 0);
684 ✗ y.assign(y.size(), 0.0);
685 ✗ for(NLuint i=0; i<M.m; ++i) {
686 ✗ const NLRowColumn& Ci = M.row[i];
687 ✗ for(NLuint jj=0; jj<Ci.size; ++jj) {
688 ✗ double a = Ci.coeff[jj].value;
689 ✗ index_t j = Ci.coeff[jj].index;
690 ✗ y[j] += a * x[i];
691 }
692 }
693 }
694 ✗ }
695
696 ✗ static void add_J_D_x(
697 vector<double>& y,
698 const NLSparseMatrix& J,
699 const vector<double>& D,
700 const vector<double>& x
701 ) {
702 geo_debug_assert(y.size() == J.m) ;
703 geo_debug_assert(D.size() == J.n) ;
704 geo_debug_assert(x.size() == J.n) ;
705
706 ✗ for(NLuint j=0; j<D.size(); ++j) {
707 ✗ const NLRowColumn& Cj = J.column[j] ;
708 ✗ for(NLuint ii=0; ii<Cj.size; ++ii) {
709 ✗ const NLCoeff& c = Cj.coeff[ii] ;
710 ✗ y[c.index] += c.value * D[j] * x[j] ;
711 }
712 }
713 ✗ }
714
715 ✗ static void add_J_D_Jt(
716 NLSparseMatrix& M,
717 const NLSparseMatrix& J,
718 const vector<double>& D
719 ) {
720 geo_debug_assert(M.m == J.m) ;
721 geo_debug_assert(M.n == J.m) ;
722 geo_debug_assert(D.size() == J.n) ;
723
724 ✗ for(NLuint j=0; j<D.size(); j++) {
725 ✗ const NLRowColumn& Cj = J.column[j];
726 ✗ for(NLuint ii1=0; ii1<Cj.size; ii1++) {
727 ✗ for(NLuint ii2=0; ii2<Cj.size; ii2++) {
728 ✗ nlSparseMatrixAdd(
729 &M,
730 ✗ NLuint(Cj.coeff[ii1].index),
731 ✗ NLuint(Cj.coeff[ii2].index),
732 ✗ Cj.coeff[ii1].value * Cj.coeff[ii2].value * D[j]
733 );
734 }
735 }
736 }
737 ✗ }
738
739 ✗ static void sub_J_D_Jt(
740 NLSparseMatrix& M,
741 const NLSparseMatrix& J,
742 const vector<double>& D
743 ) {
744 geo_debug_assert(M.m == J.m) ;
745 geo_debug_assert(M.n == J.m) ;
746 geo_debug_assert(D.size() == J.n) ;
747
748 ✗ for(NLuint j=0; j<D.size(); j++) {
749 ✗ const NLRowColumn& Cj = J.column[j];
750 ✗ for(NLuint ii1=0; ii1<Cj.size; ii1++) {
751 ✗ for(NLuint ii2=0; ii2<Cj.size; ii2++) {
752 ✗ nlSparseMatrixAdd(
753 &M,
754 ✗ NLuint(Cj.coeff[ii1].index),
755 ✗ NLuint(Cj.coeff[ii2].index),
756 ✗ -Cj.coeff[ii1].value * Cj.coeff[ii2].value * D[j]
757 );
758 }
759 }
760 }
761 ✗ }
762
763 private:
764 Mesh& mesh_;
765 Attribute<double> angle_;
766 vector<bool> v_on_border_;
767 vector<index_t> v_to_c_;
768 vector<index_t> next_c_around_v_;
769 vector<index_t> c_to_f_;
770
771 // ------ Solver parameters -----------------------------------
772 double epsilon_; // Threshold for small angles
773 double newton_tolf_; // threshold for gradient norm (rhs)
774 double newton_tolx_;
775 index_t max_newton_iter_;
776 double positive_angle_ro_;
777 double step_length_factor_;
778
779 // ------ Sizes -----------------------------------------------
780 index_t nf_ ; // Number of facets
781 index_t nalpha_ ; // Number of angles
782 index_t nint_ ; // Number of interior nodes
783 index_t nlambda_ ; // Number of constraints (= nf+2.nint)
784 index_t ntot_ ; // Total number of unknowns (= nalpha + nlamda)
785
786 // ------ ABF variables & Lagrange multipliers ----------------
787 vector<double> alpha_ ; // Unknown angles. size = nalpha
788 vector<double> lambda_ ; // Lagrange multipliers. size = nlambda
789 vector<double> beta_ ; // Optimum angles. size = nalpha
790 vector<double> w_ ; // Weights. size = nalpha
791
792 // ------ Step vectors ----------------------------------------
793 vector<double> dalpha_ ; // size = nalpha ; angles
794 vector<double> dlambda1_ ; // size = nf ; C1 part
795 vector<double> dlambda2_ ; // size = 2.nint ; C2 and C3 part
796
797 // ------ Right-hand side ( - gradients ) ---------------------
798 vector<double> b1_ ; // size = nalpha
799 vector<double> b2_ ; // size = nlambda
800
801 // ------ Jacobian of the constraints -------------------------
802 // J1 (Jacobian of constraint 1) is not stored, it is implicit
803 NLSparseMatrix J2_ ; // size = 2.nint * nalpha
804
805
806 // ------ ABF++ variables -------------------------------------
807 vector<double> Delta_inv_ ; // size = nalpha
808 vector<double> Delta_star_inv_ ; // size = nf ;
809 NLSparseMatrix J_star_ ; // size = 2.nint * nf
810 vector<double> b1_star_ ; // size = nf
811 vector<double> b2_star_ ; // size = 2.nint
812
813 // ------ Final linear system ---------------------------------
814 NLSparseMatrix M_ ; // size = 2.nint * 2.nint
815 vector<double> r_ ; // size = 2.nint
816
817 bool verbose_;
818 };
819
820 }
821
822 namespace GEO {
823
824 ✗ void mesh_compute_ABF_plus_plus(
825 Mesh& M, const std::string& attribute_name, bool verbose
826 ) {
827 // Normally, the ABFPlusPlus class can handle non-triangulated
828 // surfaces, but:
829 // 1) it does not seem to converge (to be checked, I may have a bug)
830 // 2) the angle-to-uv algorithm in LSCM needs to be
831 // adapted.
832 // (for now, we still require triangulated surfaces)
833 ✗ geo_assert(M.facets.are_simplices());
834
835 ✗ ABFPlusPlus ABF(M);
836 ABF.set_verbose(verbose);
837 ✗ ABF.parameterize(); // This computes the "angle" attribute.
838 // Now use LSCM to retrieve (u,v) coordinates from the angles.
839 ✗ mesh_compute_LSCM(M, attribute_name, false, "angle", verbose);
840 ✗ M.facet_corners.attributes().delete_attribute_store("angle");
841 ✗ }
842
843 }
844