1 /*
2 * This program is free software; you can redistribute it and/or
3 * modify it under the terms of the GNU General Public License
4 * as published by the Free Software Foundation; either version 2
5 * of the License, or (at your option) any later version.
6 *
7 * This program is distributed in the hope that it will be useful,
8 * but WITHOUT ANY WARRANTY; without even the implied warranty of
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10 * GNU General Public License for more details.
11 *
12 * You should have received a copy of the GNU General Public License
13 * along with this program; if not, write to the Free Software Foundation,
14 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15 */
16
17 /** \file
18 * \ingroup eduv
19 */
20
21 #include "MEM_guardedalloc.h"
22
23 #include "BLI_boxpack_2d.h"
24 #include "BLI_convexhull_2d.h"
25 #include "BLI_heap.h"
26 #include "BLI_math.h"
27 #include "BLI_memarena.h"
28 #include "BLI_polyfill_2d.h"
29 #include "BLI_polyfill_2d_beautify.h"
30 #include "BLI_rand.h"
31 #include "BLI_utildefines.h"
32
33 #include "uvedit_parametrizer.h"
34
35 #include <math.h>
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <string.h>
39
40 #include "BLI_sys_types.h" /* for intptr_t support */
41
42 #include "eigen_capi.h"
43
44 /* Utils */
45
46 #define param_assert(condition) \
47 if (!(condition)) { /*printf("Assertion %s:%d\n", __FILE__, __LINE__); abort();*/ \
48 } \
49 (void)0
50 #define param_warning(message) \
51 {/*printf("Warning %s:%d: %s\n", __FILE__, __LINE__, message);*/}(void)0
52
53 typedef enum PBool {
54 P_TRUE = 1,
55 P_FALSE = 0,
56 } PBool;
57
58 /* Special Purpose Hash */
59
60 typedef intptr_t PHashKey;
61
62 typedef struct PHashLink {
63 struct PHashLink *next;
64 PHashKey key;
65 } PHashLink;
66
67 typedef struct PHash {
68 PHashLink **list;
69 PHashLink **buckets;
70 int size, cursize, cursize_id;
71 } PHash;
72
73 struct PChart;
74 struct PEdge;
75 struct PFace;
76 struct PHandle;
77 struct PVert;
78
79 /* Simplices */
80
81 typedef struct PVert {
82 struct PVert *nextlink;
83
84 union PVertUnion {
85 PHashKey key; /* construct */
86 int id; /* abf/lscm matrix index */
87 float distortion; /* area smoothing */
88 HeapNode *heaplink; /* edge collapsing */
89 } u;
90
91 struct PEdge *edge;
92 float co[3];
93 float uv[2];
94 uchar flag;
95
96 } PVert;
97
98 typedef struct PEdge {
99 struct PEdge *nextlink;
100
101 union PEdgeUnion {
102 PHashKey key; /* construct */
103 int id; /* abf matrix index */
104 HeapNode *heaplink; /* fill holes */
105 struct PEdge *nextcollapse; /* simplification */
106 } u;
107
108 struct PVert *vert;
109 struct PEdge *pair;
110 struct PEdge *next;
111 struct PFace *face;
112 float *orig_uv, old_uv[2];
113 ushort flag;
114
115 } PEdge;
116
117 typedef struct PFace {
118 struct PFace *nextlink;
119
120 union PFaceUnion {
121 PHashKey key; /* construct */
122 int chart; /* construct splitting*/
123 float area3d; /* stretch */
124 int id; /* abf matrix index */
125 } u;
126
127 struct PEdge *edge;
128 uchar flag;
129 } PFace;
130
131 enum PVertFlag {
132 PVERT_PIN = 1,
133 PVERT_SELECT = 2,
134 PVERT_INTERIOR = 4,
135 PVERT_COLLAPSE = 8,
136 PVERT_SPLIT = 16,
137 };
138
139 enum PEdgeFlag {
140 PEDGE_SEAM = 1,
141 PEDGE_VERTEX_SPLIT = 2,
142 PEDGE_PIN = 4,
143 PEDGE_SELECT = 8,
144 PEDGE_DONE = 16,
145 PEDGE_FILLED = 32,
146 PEDGE_COLLAPSE = 64,
147 PEDGE_COLLAPSE_EDGE = 128,
148 PEDGE_COLLAPSE_PAIR = 256,
149 };
150
151 /* for flipping faces */
152 #define PEDGE_VERTEX_FLAGS (PEDGE_PIN)
153
154 enum PFaceFlag {
155 PFACE_CONNECTED = 1,
156 PFACE_FILLED = 2,
157 PFACE_COLLAPSE = 4,
158 };
159
160 /* Chart */
161
162 typedef struct PChart {
163 PVert *verts;
164 PEdge *edges;
165 PFace *faces;
166 int nverts, nedges, nfaces;
167
168 PVert *collapsed_verts;
169 PEdge *collapsed_edges;
170 PFace *collapsed_faces;
171
172 union PChartUnion {
173 struct PChartLscm {
174 LinearSolver *context;
175 float *abf_alpha;
176 PVert *pin1, *pin2;
177 PVert *single_pin;
178 float single_pin_area;
179 float single_pin_uv[2];
180 } lscm;
181 struct PChartPack {
182 float rescale, area;
183 float size[2] /* , trans[2] */;
184 } pack;
185 } u;
186
187 uchar flag;
188 struct PHandle *handle;
189 } PChart;
190
191 enum PChartFlag {
192 PCHART_HAS_PINS = 1,
193 };
194
195 enum PHandleState {
196 PHANDLE_STATE_ALLOCATED,
197 PHANDLE_STATE_CONSTRUCTED,
198 PHANDLE_STATE_LSCM,
199 PHANDLE_STATE_STRETCH,
200 };
201
202 typedef struct PHandle {
203 enum PHandleState state;
204 MemArena *arena;
205 MemArena *polyfill_arena;
206 Heap *polyfill_heap;
207
208 PChart *construction_chart;
209 PHash *hash_verts;
210 PHash *hash_edges;
211 PHash *hash_faces;
212
213 PChart **charts;
214 int ncharts;
215
216 float aspx, aspy;
217
218 RNG *rng;
219 float blend;
220 char do_aspect;
221 } PHandle;
222
223 /* PHash
224 * - special purpose hash that keeps all its elements in a single linked list.
225 * - after construction, this hash is thrown away, and the list remains.
226 * - removing elements is not possible efficiently.
227 */
228
229 static int PHashSizes[] = {
230 1, 3, 5, 11, 17, 37, 67, 131, 257, 521,
231 1031, 2053, 4099, 8209, 16411, 32771, 65537, 131101, 262147, 524309,
232 1048583, 2097169, 4194319, 8388617, 16777259, 33554467, 67108879, 134217757, 268435459,
233 };
234
235 #define PHASH_hash(ph, item) (((uintptr_t)(item)) % ((uint)(ph)->cursize))
236 #define PHASH_edge(v1, v2) (((v1) < (v2)) ? ((v1)*39) ^ ((v2)*31) : ((v1)*31) ^ ((v2)*39))
237
phash_new(PHashLink ** list,int sizehint)238 static PHash *phash_new(PHashLink **list, int sizehint)
239 {
240 PHash *ph = (PHash *)MEM_callocN(sizeof(PHash), "PHash");
241 ph->size = 0;
242 ph->cursize_id = 0;
243 ph->list = list;
244
245 while (PHashSizes[ph->cursize_id] < sizehint) {
246 ph->cursize_id++;
247 }
248
249 ph->cursize = PHashSizes[ph->cursize_id];
250 ph->buckets = (PHashLink **)MEM_callocN(ph->cursize * sizeof(*ph->buckets), "PHashBuckets");
251
252 return ph;
253 }
254
phash_delete(PHash * ph)255 static void phash_delete(PHash *ph)
256 {
257 MEM_freeN(ph->buckets);
258 MEM_freeN(ph);
259 }
260
phash_size(PHash * ph)261 static int phash_size(PHash *ph)
262 {
263 return ph->size;
264 }
265
phash_insert(PHash * ph,PHashLink * link)266 static void phash_insert(PHash *ph, PHashLink *link)
267 {
268 int size = ph->cursize;
269 uintptr_t hash = PHASH_hash(ph, link->key);
270 PHashLink *lookup = ph->buckets[hash];
271
272 if (lookup == NULL) {
273 /* insert in front of the list */
274 ph->buckets[hash] = link;
275 link->next = *(ph->list);
276 *(ph->list) = link;
277 }
278 else {
279 /* insert after existing element */
280 link->next = lookup->next;
281 lookup->next = link;
282 }
283
284 ph->size++;
285
286 if (ph->size > (size * 3)) {
287 PHashLink *next = NULL, *first = *(ph->list);
288
289 ph->cursize = PHashSizes[++ph->cursize_id];
290 MEM_freeN(ph->buckets);
291 ph->buckets = (PHashLink **)MEM_callocN(ph->cursize * sizeof(*ph->buckets), "PHashBuckets");
292 ph->size = 0;
293 *(ph->list) = NULL;
294
295 for (link = first; link; link = next) {
296 next = link->next;
297 phash_insert(ph, link);
298 }
299 }
300 }
301
phash_lookup(PHash * ph,PHashKey key)302 static PHashLink *phash_lookup(PHash *ph, PHashKey key)
303 {
304 PHashLink *link;
305 uintptr_t hash = PHASH_hash(ph, key);
306
307 for (link = ph->buckets[hash]; link; link = link->next) {
308 if (link->key == key) {
309 return link;
310 }
311 if (PHASH_hash(ph, link->key) != hash) {
312 return NULL;
313 }
314 }
315
316 return link;
317 }
318
phash_next(PHash * ph,PHashKey key,PHashLink * link)319 static PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
320 {
321 uintptr_t hash = PHASH_hash(ph, key);
322
323 for (link = link->next; link; link = link->next) {
324 if (link->key == key) {
325 return link;
326 }
327 if (PHASH_hash(ph, link->key) != hash) {
328 return NULL;
329 }
330 }
331
332 return link;
333 }
334
335 /* Geometry */
336
p_vec_angle_cos(const float v1[3],const float v2[3],const float v3[3])337 static float p_vec_angle_cos(const float v1[3], const float v2[3], const float v3[3])
338 {
339 float d1[3], d2[3];
340
341 d1[0] = v1[0] - v2[0];
342 d1[1] = v1[1] - v2[1];
343 d1[2] = v1[2] - v2[2];
344
345 d2[0] = v3[0] - v2[0];
346 d2[1] = v3[1] - v2[1];
347 d2[2] = v3[2] - v2[2];
348
349 normalize_v3(d1);
350 normalize_v3(d2);
351
352 return d1[0] * d2[0] + d1[1] * d2[1] + d1[2] * d2[2];
353 }
354
p_vec_angle(const float v1[3],const float v2[3],const float v3[3])355 static float p_vec_angle(const float v1[3], const float v2[3], const float v3[3])
356 {
357 float dot = p_vec_angle_cos(v1, v2, v3);
358
359 if (dot <= -1.0f) {
360 return (float)M_PI;
361 }
362 if (dot >= 1.0f) {
363 return 0.0f;
364 }
365 return acosf(dot);
366 }
367
p_vec2_angle(const float v1[2],const float v2[2],const float v3[2])368 static float p_vec2_angle(const float v1[2], const float v2[2], const float v3[2])
369 {
370 float u1[3], u2[3], u3[3];
371
372 u1[0] = v1[0];
373 u1[1] = v1[1];
374 u1[2] = 0.0f;
375 u2[0] = v2[0];
376 u2[1] = v2[1];
377 u2[2] = 0.0f;
378 u3[0] = v3[0];
379 u3[1] = v3[1];
380 u3[2] = 0.0f;
381
382 return p_vec_angle(u1, u2, u3);
383 }
384
p_triangle_angles(const float v1[3],const float v2[3],const float v3[3],float * r_a1,float * r_a2,float * r_a3)385 static void p_triangle_angles(
386 const float v1[3], const float v2[3], const float v3[3], float *r_a1, float *r_a2, float *r_a3)
387 {
388 *r_a1 = p_vec_angle(v3, v1, v2);
389 *r_a2 = p_vec_angle(v1, v2, v3);
390 *r_a3 = (float)M_PI - *r_a2 - *r_a1;
391 }
392
p_face_angles(PFace * f,float * r_a1,float * r_a2,float * r_a3)393 static void p_face_angles(PFace *f, float *r_a1, float *r_a2, float *r_a3)
394 {
395 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
396 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
397
398 p_triangle_angles(v1->co, v2->co, v3->co, r_a1, r_a2, r_a3);
399 }
400
p_face_area(PFace * f)401 static float p_face_area(PFace *f)
402 {
403 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
404 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
405
406 return area_tri_v3(v1->co, v2->co, v3->co);
407 }
408
p_area_signed(const float v1[2],const float v2[2],const float v3[2])409 static float p_area_signed(const float v1[2], const float v2[2], const float v3[2])
410 {
411 return 0.5f * (((v2[0] - v1[0]) * (v3[1] - v1[1])) - ((v3[0] - v1[0]) * (v2[1] - v1[1])));
412 }
413
p_face_uv_area_signed(PFace * f)414 static float p_face_uv_area_signed(PFace *f)
415 {
416 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
417 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
418
419 return 0.5f * (((v2->uv[0] - v1->uv[0]) * (v3->uv[1] - v1->uv[1])) -
420 ((v3->uv[0] - v1->uv[0]) * (v2->uv[1] - v1->uv[1])));
421 }
422
p_edge_length(PEdge * e)423 static float p_edge_length(PEdge *e)
424 {
425 PVert *v1 = e->vert, *v2 = e->next->vert;
426 float d[3];
427
428 d[0] = v2->co[0] - v1->co[0];
429 d[1] = v2->co[1] - v1->co[1];
430 d[2] = v2->co[2] - v1->co[2];
431
432 return sqrtf(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
433 }
434
p_edge_uv_length(PEdge * e)435 static float p_edge_uv_length(PEdge *e)
436 {
437 PVert *v1 = e->vert, *v2 = e->next->vert;
438 float d[3];
439
440 d[0] = v2->uv[0] - v1->uv[0];
441 d[1] = v2->uv[1] - v1->uv[1];
442
443 return sqrtf(d[0] * d[0] + d[1] * d[1]);
444 }
445
p_chart_uv_bbox(PChart * chart,float minv[2],float maxv[2])446 static void p_chart_uv_bbox(PChart *chart, float minv[2], float maxv[2])
447 {
448 PVert *v;
449
450 INIT_MINMAX2(minv, maxv);
451
452 for (v = chart->verts; v; v = v->nextlink) {
453 minmax_v2v2_v2(minv, maxv, v->uv);
454 }
455 }
456
p_chart_uv_area(PChart * chart)457 static float p_chart_uv_area(PChart *chart)
458 {
459 float area = 0.0f;
460
461 for (PFace *f = chart->faces; f; f = f->nextlink) {
462 area += fabsf(p_face_uv_area_signed(f));
463 }
464
465 return area;
466 }
467
p_chart_uv_scale(PChart * chart,float scale)468 static void p_chart_uv_scale(PChart *chart, float scale)
469 {
470 PVert *v;
471
472 for (v = chart->verts; v; v = v->nextlink) {
473 v->uv[0] *= scale;
474 v->uv[1] *= scale;
475 }
476 }
477
p_chart_uv_scale_xy(PChart * chart,float x,float y)478 static void p_chart_uv_scale_xy(PChart *chart, float x, float y)
479 {
480 PVert *v;
481
482 for (v = chart->verts; v; v = v->nextlink) {
483 v->uv[0] *= x;
484 v->uv[1] *= y;
485 }
486 }
487
p_chart_uv_translate(PChart * chart,const float trans[2])488 static void p_chart_uv_translate(PChart *chart, const float trans[2])
489 {
490 PVert *v;
491
492 for (v = chart->verts; v; v = v->nextlink) {
493 v->uv[0] += trans[0];
494 v->uv[1] += trans[1];
495 }
496 }
497
p_chart_uv_transform(PChart * chart,const float mat[2][2])498 static void p_chart_uv_transform(PChart *chart, const float mat[2][2])
499 {
500 PVert *v;
501
502 for (v = chart->verts; v; v = v->nextlink) {
503 mul_m2_v2(mat, v->uv);
504 }
505 }
506
p_chart_uv_to_array(PChart * chart,float (* points)[2])507 static void p_chart_uv_to_array(PChart *chart, float (*points)[2])
508 {
509 PVert *v;
510 uint i = 0;
511
512 for (v = chart->verts; v; v = v->nextlink) {
513 copy_v2_v2(points[i++], v->uv);
514 }
515 }
516
UNUSED_FUNCTION(p_chart_uv_from_array)517 static void UNUSED_FUNCTION(p_chart_uv_from_array)(PChart *chart, float (*points)[2])
518 {
519 PVert *v;
520 uint i = 0;
521
522 for (v = chart->verts; v; v = v->nextlink) {
523 copy_v2_v2(v->uv, points[i++]);
524 }
525 }
526
p_intersect_line_2d_dir(const float v1[2],const float dir1[2],const float v2[2],const float dir2[2],float r_isect[2])527 static PBool p_intersect_line_2d_dir(const float v1[2],
528 const float dir1[2],
529 const float v2[2],
530 const float dir2[2],
531 float r_isect[2])
532 {
533 float lmbda, div;
534
535 div = dir2[0] * dir1[1] - dir2[1] * dir1[0];
536
537 if (div == 0.0f) {
538 return P_FALSE;
539 }
540
541 lmbda = ((v1[1] - v2[1]) * dir1[0] - (v1[0] - v2[0]) * dir1[1]) / div;
542 r_isect[0] = v1[0] + lmbda * dir2[0];
543 r_isect[1] = v1[1] + lmbda * dir2[1];
544
545 return P_TRUE;
546 }
547
548 #if 0
549 static PBool p_intersect_line_2d(const float v1[2],
550 const float v2[2],
551 const float v3[2],
552 const float v4[2],
553 const float r_isect[2])
554 {
555 float dir1[2], dir2[2];
556
557 dir1[0] = v4[0] - v3[0];
558 dir1[1] = v4[1] - v3[1];
559
560 dir2[0] = v2[0] - v1[0];
561 dir2[1] = v2[1] - v1[1];
562
563 if (!p_intersect_line_2d_dir(v1, dir1, v2, dir2, isect)) {
564 /* parallel - should never happen in theory for polygon kernel, but
565 * let's give a point nearby in case things go wrong */
566 isect[0] = (v1[0] + v2[0]) * 0.5f;
567 isect[1] = (v1[1] + v2[1]) * 0.5f;
568 return P_FALSE;
569 }
570
571 return P_TRUE;
572 }
573 #endif
574
575 /* Topological Utilities */
576
p_wheel_edge_next(PEdge * e)577 static PEdge *p_wheel_edge_next(PEdge *e)
578 {
579 return e->next->next->pair;
580 }
581
p_wheel_edge_prev(PEdge * e)582 static PEdge *p_wheel_edge_prev(PEdge *e)
583 {
584 return (e->pair) ? e->pair->next : NULL;
585 }
586
p_boundary_edge_next(PEdge * e)587 static PEdge *p_boundary_edge_next(PEdge *e)
588 {
589 return e->next->vert->edge;
590 }
591
p_boundary_edge_prev(PEdge * e)592 static PEdge *p_boundary_edge_prev(PEdge *e)
593 {
594 PEdge *we = e, *last;
595
596 do {
597 last = we;
598 we = p_wheel_edge_next(we);
599 } while (we && (we != e));
600
601 return last->next->next;
602 }
603
p_vert_interior(PVert * v)604 static PBool p_vert_interior(PVert *v)
605 {
606 return (v->edge->pair != NULL);
607 }
608
p_face_flip(PFace * f)609 static void p_face_flip(PFace *f)
610 {
611 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
612 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
613 int f1 = e1->flag, f2 = e2->flag, f3 = e3->flag;
614 float *orig_uv1 = e1->orig_uv, *orig_uv2 = e2->orig_uv, *orig_uv3 = e3->orig_uv;
615
616 e1->vert = v2;
617 e1->next = e3;
618 e1->orig_uv = orig_uv2;
619 e1->flag = (f1 & ~PEDGE_VERTEX_FLAGS) | (f2 & PEDGE_VERTEX_FLAGS);
620
621 e2->vert = v3;
622 e2->next = e1;
623 e2->orig_uv = orig_uv3;
624 e2->flag = (f2 & ~PEDGE_VERTEX_FLAGS) | (f3 & PEDGE_VERTEX_FLAGS);
625
626 e3->vert = v1;
627 e3->next = e2;
628 e3->orig_uv = orig_uv1;
629 e3->flag = (f3 & ~PEDGE_VERTEX_FLAGS) | (f1 & PEDGE_VERTEX_FLAGS);
630 }
631
632 #if 0
633 static void p_chart_topological_sanity_check(PChart *chart)
634 {
635 PVert *v;
636 PEdge *e;
637
638 for (v = chart->verts; v; v = v->nextlink) {
639 param_test_equals_ptr("v->edge->vert", v, v->edge->vert);
640 }
641
642 for (e = chart->edges; e; e = e->nextlink) {
643 if (e->pair) {
644 param_test_equals_ptr("e->pair->pair", e, e->pair->pair);
645 param_test_equals_ptr("pair->vert", e->vert, e->pair->next->vert);
646 param_test_equals_ptr("pair->next->vert", e->next->vert, e->pair->vert);
647 }
648 }
649 }
650 #endif
651
652 /* Loading / Flushing */
653
p_vert_load_pin_select_uvs(PHandle * handle,PVert * v)654 static void p_vert_load_pin_select_uvs(PHandle *handle, PVert *v)
655 {
656 PEdge *e;
657 int nedges = 0, npins = 0;
658 float pinuv[2];
659
660 v->uv[0] = v->uv[1] = 0.0f;
661 pinuv[0] = pinuv[1] = 0.0f;
662 e = v->edge;
663 do {
664 if (e->orig_uv) {
665 if (e->flag & PEDGE_SELECT) {
666 v->flag |= PVERT_SELECT;
667 }
668
669 if (e->flag & PEDGE_PIN) {
670 pinuv[0] += e->orig_uv[0] * handle->aspx;
671 pinuv[1] += e->orig_uv[1] * handle->aspy;
672 npins++;
673 }
674 else {
675 v->uv[0] += e->orig_uv[0] * handle->aspx;
676 v->uv[1] += e->orig_uv[1] * handle->aspy;
677 }
678
679 nedges++;
680 }
681
682 e = p_wheel_edge_next(e);
683 } while (e && e != (v->edge));
684
685 if (npins > 0) {
686 v->uv[0] = pinuv[0] / npins;
687 v->uv[1] = pinuv[1] / npins;
688 v->flag |= PVERT_PIN;
689 }
690 else if (nedges > 0) {
691 v->uv[0] /= nedges;
692 v->uv[1] /= nedges;
693 }
694 }
695
p_flush_uvs(PHandle * handle,PChart * chart)696 static void p_flush_uvs(PHandle *handle, PChart *chart)
697 {
698 PEdge *e;
699
700 for (e = chart->edges; e; e = e->nextlink) {
701 if (e->orig_uv) {
702 e->orig_uv[0] = e->vert->uv[0] / handle->aspx;
703 e->orig_uv[1] = e->vert->uv[1] / handle->aspy;
704 }
705 }
706 }
707
p_flush_uvs_blend(PHandle * handle,PChart * chart,float blend)708 static void p_flush_uvs_blend(PHandle *handle, PChart *chart, float blend)
709 {
710 PEdge *e;
711 float invblend = 1.0f - blend;
712
713 for (e = chart->edges; e; e = e->nextlink) {
714 if (e->orig_uv) {
715 e->orig_uv[0] = blend * e->old_uv[0] + invblend * e->vert->uv[0] / handle->aspx;
716 e->orig_uv[1] = blend * e->old_uv[1] + invblend * e->vert->uv[1] / handle->aspy;
717 }
718 }
719 }
720
p_face_backup_uvs(PFace * f)721 static void p_face_backup_uvs(PFace *f)
722 {
723 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
724
725 if (e1->orig_uv) {
726 e1->old_uv[0] = e1->orig_uv[0];
727 e1->old_uv[1] = e1->orig_uv[1];
728 }
729 if (e2->orig_uv) {
730 e2->old_uv[0] = e2->orig_uv[0];
731 e2->old_uv[1] = e2->orig_uv[1];
732 }
733 if (e3->orig_uv) {
734 e3->old_uv[0] = e3->orig_uv[0];
735 e3->old_uv[1] = e3->orig_uv[1];
736 }
737 }
738
p_face_restore_uvs(PFace * f)739 static void p_face_restore_uvs(PFace *f)
740 {
741 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
742
743 if (e1->orig_uv) {
744 e1->orig_uv[0] = e1->old_uv[0];
745 e1->orig_uv[1] = e1->old_uv[1];
746 }
747 if (e2->orig_uv) {
748 e2->orig_uv[0] = e2->old_uv[0];
749 e2->orig_uv[1] = e2->old_uv[1];
750 }
751 if (e3->orig_uv) {
752 e3->orig_uv[0] = e3->old_uv[0];
753 e3->orig_uv[1] = e3->old_uv[1];
754 }
755 }
756
757 /* Construction (use only during construction, relies on u.key being set */
758
p_vert_add(PHandle * handle,PHashKey key,const float co[3],PEdge * e)759 static PVert *p_vert_add(PHandle *handle, PHashKey key, const float co[3], PEdge *e)
760 {
761 PVert *v = (PVert *)BLI_memarena_alloc(handle->arena, sizeof(*v));
762 copy_v3_v3(v->co, co);
763
764 /* Sanity check, a single nan/inf point causes the entire result to be invalid.
765 * Note that values within the calculation may _become_ non-finite,
766 * so the rest of the code still needs to take this possibility into account. */
767 for (int i = 0; i < 3; i++) {
768 if (UNLIKELY(!isfinite(v->co[i]))) {
769 v->co[i] = 0.0f;
770 }
771 }
772
773 v->u.key = key;
774 v->edge = e;
775 v->flag = 0;
776
777 phash_insert(handle->hash_verts, (PHashLink *)v);
778
779 return v;
780 }
781
p_vert_lookup(PHandle * handle,PHashKey key,const float co[3],PEdge * e)782 static PVert *p_vert_lookup(PHandle *handle, PHashKey key, const float co[3], PEdge *e)
783 {
784 PVert *v = (PVert *)phash_lookup(handle->hash_verts, key);
785
786 if (v) {
787 return v;
788 }
789 return p_vert_add(handle, key, co, e);
790 }
791
p_vert_copy(PChart * chart,PVert * v)792 static PVert *p_vert_copy(PChart *chart, PVert *v)
793 {
794 PVert *nv = (PVert *)BLI_memarena_alloc(chart->handle->arena, sizeof(*nv));
795
796 copy_v3_v3(nv->co, v->co);
797 nv->uv[0] = v->uv[0];
798 nv->uv[1] = v->uv[1];
799 nv->u.key = v->u.key;
800 nv->edge = v->edge;
801 nv->flag = v->flag;
802
803 return nv;
804 }
805
p_edge_lookup(PHandle * handle,const PHashKey * vkeys)806 static PEdge *p_edge_lookup(PHandle *handle, const PHashKey *vkeys)
807 {
808 PHashKey key = PHASH_edge(vkeys[0], vkeys[1]);
809 PEdge *e = (PEdge *)phash_lookup(handle->hash_edges, key);
810
811 while (e) {
812 if ((e->vert->u.key == vkeys[0]) && (e->next->vert->u.key == vkeys[1])) {
813 return e;
814 }
815 if ((e->vert->u.key == vkeys[1]) && (e->next->vert->u.key == vkeys[0])) {
816 return e;
817 }
818
819 e = (PEdge *)phash_next(handle->hash_edges, key, (PHashLink *)e);
820 }
821
822 return NULL;
823 }
824
p_face_exists(ParamHandle * phandle,ParamKey * pvkeys,int i1,int i2,int i3)825 static int p_face_exists(ParamHandle *phandle, ParamKey *pvkeys, int i1, int i2, int i3)
826 {
827 PHandle *handle = (PHandle *)phandle;
828 PHashKey *vkeys = (PHashKey *)pvkeys;
829 PHashKey key = PHASH_edge(vkeys[i1], vkeys[i2]);
830 PEdge *e = (PEdge *)phash_lookup(handle->hash_edges, key);
831
832 while (e) {
833 if ((e->vert->u.key == vkeys[i1]) && (e->next->vert->u.key == vkeys[i2])) {
834 if (e->next->next->vert->u.key == vkeys[i3]) {
835 return P_TRUE;
836 }
837 }
838 else if ((e->vert->u.key == vkeys[i2]) && (e->next->vert->u.key == vkeys[i1])) {
839 if (e->next->next->vert->u.key == vkeys[i3]) {
840 return P_TRUE;
841 }
842 }
843
844 e = (PEdge *)phash_next(handle->hash_edges, key, (PHashLink *)e);
845 }
846
847 return P_FALSE;
848 }
849
p_chart_new(PHandle * handle)850 static PChart *p_chart_new(PHandle *handle)
851 {
852 PChart *chart = (PChart *)MEM_callocN(sizeof(*chart), "PChart");
853 chart->handle = handle;
854
855 return chart;
856 }
857
p_chart_delete(PChart * chart)858 static void p_chart_delete(PChart *chart)
859 {
860 /* the actual links are free by memarena */
861 MEM_freeN(chart);
862 }
863
p_edge_implicit_seam(PEdge * e,PEdge * ep)864 static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
865 {
866 float *uv1, *uv2, *uvp1, *uvp2;
867 float limit[2];
868
869 limit[0] = 0.00001;
870 limit[1] = 0.00001;
871
872 uv1 = e->orig_uv;
873 uv2 = e->next->orig_uv;
874
875 if (e->vert->u.key == ep->vert->u.key) {
876 uvp1 = ep->orig_uv;
877 uvp2 = ep->next->orig_uv;
878 }
879 else {
880 uvp1 = ep->next->orig_uv;
881 uvp2 = ep->orig_uv;
882 }
883
884 if ((fabsf(uv1[0] - uvp1[0]) > limit[0]) || (fabsf(uv1[1] - uvp1[1]) > limit[1])) {
885 e->flag |= PEDGE_SEAM;
886 ep->flag |= PEDGE_SEAM;
887 return P_TRUE;
888 }
889 if ((fabsf(uv2[0] - uvp2[0]) > limit[0]) || (fabsf(uv2[1] - uvp2[1]) > limit[1])) {
890 e->flag |= PEDGE_SEAM;
891 ep->flag |= PEDGE_SEAM;
892 return P_TRUE;
893 }
894
895 return P_FALSE;
896 }
897
p_edge_has_pair(PHandle * handle,PEdge * e,PBool impl,PEdge ** r_pair)898 static PBool p_edge_has_pair(PHandle *handle, PEdge *e, PBool impl, PEdge **r_pair)
899 {
900 PHashKey key;
901 PEdge *pe;
902 PVert *v1, *v2;
903 PHashKey key1 = e->vert->u.key;
904 PHashKey key2 = e->next->vert->u.key;
905
906 if (e->flag & PEDGE_SEAM) {
907 return P_FALSE;
908 }
909
910 key = PHASH_edge(key1, key2);
911 pe = (PEdge *)phash_lookup(handle->hash_edges, key);
912 *r_pair = NULL;
913
914 while (pe) {
915 if (pe != e) {
916 v1 = pe->vert;
917 v2 = pe->next->vert;
918
919 if (((v1->u.key == key1) && (v2->u.key == key2)) ||
920 ((v1->u.key == key2) && (v2->u.key == key1))) {
921
922 /* don't connect seams and t-junctions */
923 if ((pe->flag & PEDGE_SEAM) || *r_pair || (impl && p_edge_implicit_seam(e, pe))) {
924 *r_pair = NULL;
925 return P_FALSE;
926 }
927
928 *r_pair = pe;
929 }
930 }
931
932 pe = (PEdge *)phash_next(handle->hash_edges, key, (PHashLink *)pe);
933 }
934
935 if (*r_pair && (e->vert == (*r_pair)->vert)) {
936 if ((*r_pair)->next->pair || (*r_pair)->next->next->pair) {
937 /* non unfoldable, maybe mobius ring or klein bottle */
938 *r_pair = NULL;
939 return P_FALSE;
940 }
941 }
942
943 return (*r_pair != NULL);
944 }
945
p_edge_connect_pair(PHandle * handle,PEdge * e,PBool impl,PEdge *** stack)946 static PBool p_edge_connect_pair(PHandle *handle, PEdge *e, PBool impl, PEdge ***stack)
947 {
948 PEdge *pair = NULL;
949
950 if (!e->pair && p_edge_has_pair(handle, e, impl, &pair)) {
951 if (e->vert == pair->vert) {
952 p_face_flip(pair->face);
953 }
954
955 e->pair = pair;
956 pair->pair = e;
957
958 if (!(pair->face->flag & PFACE_CONNECTED)) {
959 **stack = pair;
960 (*stack)++;
961 }
962 }
963
964 return (e->pair != NULL);
965 }
966
p_connect_pairs(PHandle * handle,PBool impl)967 static int p_connect_pairs(PHandle *handle, PBool impl)
968 {
969 PEdge **stackbase = MEM_mallocN(sizeof(*stackbase) * phash_size(handle->hash_faces),
970 "Pstackbase");
971 PEdge **stack = stackbase;
972 PFace *f, *first;
973 PEdge *e, *e1, *e2;
974 PChart *chart = handle->construction_chart;
975 int ncharts = 0;
976
977 /* connect pairs, count edges, set vertex-edge pointer to a pairless edge */
978 for (first = chart->faces; first; first = first->nextlink) {
979 if (first->flag & PFACE_CONNECTED) {
980 continue;
981 }
982
983 *stack = first->edge;
984 stack++;
985
986 while (stack != stackbase) {
987 stack--;
988 e = *stack;
989 e1 = e->next;
990 e2 = e1->next;
991
992 f = e->face;
993 f->flag |= PFACE_CONNECTED;
994
995 /* assign verts to charts so we can sort them later */
996 f->u.chart = ncharts;
997
998 if (!p_edge_connect_pair(handle, e, impl, &stack)) {
999 e->vert->edge = e;
1000 }
1001 if (!p_edge_connect_pair(handle, e1, impl, &stack)) {
1002 e1->vert->edge = e1;
1003 }
1004 if (!p_edge_connect_pair(handle, e2, impl, &stack)) {
1005 e2->vert->edge = e2;
1006 }
1007 }
1008
1009 ncharts++;
1010 }
1011
1012 MEM_freeN(stackbase);
1013
1014 return ncharts;
1015 }
1016
p_split_vert(PChart * chart,PEdge * e)1017 static void p_split_vert(PChart *chart, PEdge *e)
1018 {
1019 PEdge *we, *lastwe = NULL;
1020 PVert *v = e->vert;
1021 PBool copy = P_TRUE;
1022
1023 if (e->flag & PEDGE_PIN) {
1024 chart->flag |= PCHART_HAS_PINS;
1025 }
1026
1027 if (e->flag & PEDGE_VERTEX_SPLIT) {
1028 return;
1029 }
1030
1031 /* rewind to start */
1032 lastwe = e;
1033 for (we = p_wheel_edge_prev(e); we && (we != e); we = p_wheel_edge_prev(we)) {
1034 lastwe = we;
1035 }
1036
1037 /* go over all edges in wheel */
1038 for (we = lastwe; we; we = p_wheel_edge_next(we)) {
1039 if (we->flag & PEDGE_VERTEX_SPLIT) {
1040 break;
1041 }
1042
1043 we->flag |= PEDGE_VERTEX_SPLIT;
1044
1045 if (we == v->edge) {
1046 /* found it, no need to copy */
1047 copy = P_FALSE;
1048 v->nextlink = chart->verts;
1049 chart->verts = v;
1050 chart->nverts++;
1051 }
1052 }
1053
1054 if (copy) {
1055 /* not found, copying */
1056 v->flag |= PVERT_SPLIT;
1057 v = p_vert_copy(chart, v);
1058 v->flag |= PVERT_SPLIT;
1059
1060 v->nextlink = chart->verts;
1061 chart->verts = v;
1062 chart->nverts++;
1063
1064 v->edge = lastwe;
1065
1066 we = lastwe;
1067 do {
1068 we->vert = v;
1069 we = p_wheel_edge_next(we);
1070 } while (we && (we != lastwe));
1071 }
1072 }
1073
p_split_charts(PHandle * handle,PChart * chart,int ncharts)1074 static PChart **p_split_charts(PHandle *handle, PChart *chart, int ncharts)
1075 {
1076 PChart **charts = MEM_mallocN(sizeof(*charts) * ncharts, "PCharts"), *nchart;
1077 PFace *f, *nextf;
1078 int i;
1079
1080 for (i = 0; i < ncharts; i++) {
1081 charts[i] = p_chart_new(handle);
1082 }
1083
1084 f = chart->faces;
1085 while (f) {
1086 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1087 nextf = f->nextlink;
1088
1089 nchart = charts[f->u.chart];
1090
1091 f->nextlink = nchart->faces;
1092 nchart->faces = f;
1093 e1->nextlink = nchart->edges;
1094 nchart->edges = e1;
1095 e2->nextlink = nchart->edges;
1096 nchart->edges = e2;
1097 e3->nextlink = nchart->edges;
1098 nchart->edges = e3;
1099
1100 nchart->nfaces++;
1101 nchart->nedges += 3;
1102
1103 p_split_vert(nchart, e1);
1104 p_split_vert(nchart, e2);
1105 p_split_vert(nchart, e3);
1106
1107 f = nextf;
1108 }
1109
1110 return charts;
1111 }
1112
p_face_add(PHandle * handle)1113 static PFace *p_face_add(PHandle *handle)
1114 {
1115 PFace *f;
1116 PEdge *e1, *e2, *e3;
1117
1118 /* allocate */
1119 f = (PFace *)BLI_memarena_alloc(handle->arena, sizeof(*f));
1120 f->flag = 0; /* init ! */
1121
1122 e1 = (PEdge *)BLI_memarena_alloc(handle->arena, sizeof(*e1));
1123 e2 = (PEdge *)BLI_memarena_alloc(handle->arena, sizeof(*e2));
1124 e3 = (PEdge *)BLI_memarena_alloc(handle->arena, sizeof(*e3));
1125
1126 /* set up edges */
1127 f->edge = e1;
1128 e1->face = e2->face = e3->face = f;
1129
1130 e1->next = e2;
1131 e2->next = e3;
1132 e3->next = e1;
1133
1134 e1->pair = NULL;
1135 e2->pair = NULL;
1136 e3->pair = NULL;
1137
1138 e1->flag = 0;
1139 e2->flag = 0;
1140 e3->flag = 0;
1141
1142 return f;
1143 }
1144
p_face_add_construct(PHandle * handle,ParamKey key,const ParamKey * vkeys,float * co[4],float * uv[4],int i1,int i2,int i3,const ParamBool * pin,const ParamBool * select)1145 static PFace *p_face_add_construct(PHandle *handle,
1146 ParamKey key,
1147 const ParamKey *vkeys,
1148 float *co[4],
1149 float *uv[4],
1150 int i1,
1151 int i2,
1152 int i3,
1153 const ParamBool *pin,
1154 const ParamBool *select)
1155 {
1156 PFace *f = p_face_add(handle);
1157 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1158
1159 e1->vert = p_vert_lookup(handle, vkeys[i1], co[i1], e1);
1160 e2->vert = p_vert_lookup(handle, vkeys[i2], co[i2], e2);
1161 e3->vert = p_vert_lookup(handle, vkeys[i3], co[i3], e3);
1162
1163 e1->orig_uv = uv[i1];
1164 e2->orig_uv = uv[i2];
1165 e3->orig_uv = uv[i3];
1166
1167 if (pin) {
1168 if (pin[i1]) {
1169 e1->flag |= PEDGE_PIN;
1170 }
1171 if (pin[i2]) {
1172 e2->flag |= PEDGE_PIN;
1173 }
1174 if (pin[i3]) {
1175 e3->flag |= PEDGE_PIN;
1176 }
1177 }
1178
1179 if (select) {
1180 if (select[i1]) {
1181 e1->flag |= PEDGE_SELECT;
1182 }
1183 if (select[i2]) {
1184 e2->flag |= PEDGE_SELECT;
1185 }
1186 if (select[i3]) {
1187 e3->flag |= PEDGE_SELECT;
1188 }
1189 }
1190
1191 /* insert into hash */
1192 f->u.key = key;
1193 phash_insert(handle->hash_faces, (PHashLink *)f);
1194
1195 e1->u.key = PHASH_edge(vkeys[i1], vkeys[i2]);
1196 e2->u.key = PHASH_edge(vkeys[i2], vkeys[i3]);
1197 e3->u.key = PHASH_edge(vkeys[i3], vkeys[i1]);
1198
1199 phash_insert(handle->hash_edges, (PHashLink *)e1);
1200 phash_insert(handle->hash_edges, (PHashLink *)e2);
1201 phash_insert(handle->hash_edges, (PHashLink *)e3);
1202
1203 return f;
1204 }
1205
p_face_add_fill(PChart * chart,PVert * v1,PVert * v2,PVert * v3)1206 static PFace *p_face_add_fill(PChart *chart, PVert *v1, PVert *v2, PVert *v3)
1207 {
1208 PFace *f = p_face_add(chart->handle);
1209 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1210
1211 e1->vert = v1;
1212 e2->vert = v2;
1213 e3->vert = v3;
1214
1215 e1->orig_uv = e2->orig_uv = e3->orig_uv = NULL;
1216
1217 f->nextlink = chart->faces;
1218 chart->faces = f;
1219 e1->nextlink = chart->edges;
1220 chart->edges = e1;
1221 e2->nextlink = chart->edges;
1222 chart->edges = e2;
1223 e3->nextlink = chart->edges;
1224 chart->edges = e3;
1225
1226 chart->nfaces++;
1227 chart->nedges += 3;
1228
1229 return f;
1230 }
1231
p_quad_split_direction(PHandle * handle,float ** co,PHashKey * vkeys)1232 static PBool p_quad_split_direction(PHandle *handle, float **co, PHashKey *vkeys)
1233 {
1234 /* Slight bias to prefer one edge over the other in case they are equal, so
1235 * that in symmetric models we choose the same split direction instead of
1236 * depending on floating point errors to decide. */
1237 float bias = 1.0f + 1e-6f;
1238 float fac = len_v3v3(co[0], co[2]) * bias - len_v3v3(co[1], co[3]);
1239 PBool dir = (fac <= 0.0f);
1240
1241 /* The face exists check is there because of a special case:
1242 * when two quads share three vertices, they can each be split into two triangles,
1243 * resulting in two identical triangles. For example in Suzanne's nose. */
1244 if (dir) {
1245 if (p_face_exists(handle, vkeys, 0, 1, 2) || p_face_exists(handle, vkeys, 0, 2, 3)) {
1246 return !dir;
1247 }
1248 }
1249 else {
1250 if (p_face_exists(handle, vkeys, 0, 1, 3) || p_face_exists(handle, vkeys, 1, 2, 3)) {
1251 return !dir;
1252 }
1253 }
1254
1255 return dir;
1256 }
1257
1258 /* Construction: boundary filling */
1259
p_chart_boundaries(PChart * chart,int * r_nboundaries,PEdge ** r_outer)1260 static void p_chart_boundaries(PChart *chart, int *r_nboundaries, PEdge **r_outer)
1261 {
1262 PEdge *e, *be;
1263 float len, maxlen = -1.0;
1264
1265 if (r_nboundaries) {
1266 *r_nboundaries = 0;
1267 }
1268 if (r_outer) {
1269 *r_outer = NULL;
1270 }
1271
1272 for (e = chart->edges; e; e = e->nextlink) {
1273 if (e->pair || (e->flag & PEDGE_DONE)) {
1274 continue;
1275 }
1276
1277 if (r_nboundaries) {
1278 (*r_nboundaries)++;
1279 }
1280
1281 len = 0.0f;
1282
1283 be = e;
1284 do {
1285 be->flag |= PEDGE_DONE;
1286 len += p_edge_length(be);
1287 be = be->next->vert->edge;
1288 } while (be != e);
1289
1290 if (r_outer && (len > maxlen)) {
1291 *r_outer = e;
1292 maxlen = len;
1293 }
1294 }
1295
1296 for (e = chart->edges; e; e = e->nextlink) {
1297 e->flag &= ~PEDGE_DONE;
1298 }
1299 }
1300
p_edge_boundary_angle(PEdge * e)1301 static float p_edge_boundary_angle(PEdge *e)
1302 {
1303 PEdge *we;
1304 PVert *v, *v1, *v2;
1305 float angle;
1306 int n = 0;
1307
1308 v = e->vert;
1309
1310 /* concave angle check -- could be better */
1311 angle = M_PI;
1312
1313 we = v->edge;
1314 do {
1315 v1 = we->next->vert;
1316 v2 = we->next->next->vert;
1317 angle -= p_vec_angle(v1->co, v->co, v2->co);
1318
1319 we = we->next->next->pair;
1320 n++;
1321 } while (we && (we != v->edge));
1322
1323 return angle;
1324 }
1325
p_chart_fill_boundary(PChart * chart,PEdge * be,int nedges)1326 static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
1327 {
1328 PEdge *e, *e1, *e2;
1329
1330 PFace *f;
1331 struct Heap *heap = BLI_heap_new();
1332 float angle;
1333
1334 e = be;
1335 do {
1336 angle = p_edge_boundary_angle(e);
1337 e->u.heaplink = BLI_heap_insert(heap, angle, e);
1338
1339 e = p_boundary_edge_next(e);
1340 } while (e != be);
1341
1342 if (nedges == 2) {
1343 /* no real boundary, but an isolated seam */
1344 e = be->next->vert->edge;
1345 e->pair = be;
1346 be->pair = e;
1347
1348 BLI_heap_remove(heap, e->u.heaplink);
1349 BLI_heap_remove(heap, be->u.heaplink);
1350 }
1351 else {
1352 while (nedges > 2) {
1353 PEdge *ne, *ne1, *ne2;
1354
1355 e = (PEdge *)BLI_heap_pop_min(heap);
1356
1357 e1 = p_boundary_edge_prev(e);
1358 e2 = p_boundary_edge_next(e);
1359
1360 BLI_heap_remove(heap, e1->u.heaplink);
1361 BLI_heap_remove(heap, e2->u.heaplink);
1362 e->u.heaplink = e1->u.heaplink = e2->u.heaplink = NULL;
1363
1364 e->flag |= PEDGE_FILLED;
1365 e1->flag |= PEDGE_FILLED;
1366
1367 f = p_face_add_fill(chart, e->vert, e1->vert, e2->vert);
1368 f->flag |= PFACE_FILLED;
1369
1370 ne = f->edge->next->next;
1371 ne1 = f->edge;
1372 ne2 = f->edge->next;
1373
1374 ne->flag = ne1->flag = ne2->flag = PEDGE_FILLED;
1375
1376 e->pair = ne;
1377 ne->pair = e;
1378 e1->pair = ne1;
1379 ne1->pair = e1;
1380
1381 ne->vert = e2->vert;
1382 ne1->vert = e->vert;
1383 ne2->vert = e1->vert;
1384
1385 if (nedges == 3) {
1386 e2->pair = ne2;
1387 ne2->pair = e2;
1388 }
1389 else {
1390 ne2->vert->edge = ne2;
1391
1392 ne2->u.heaplink = BLI_heap_insert(heap, p_edge_boundary_angle(ne2), ne2);
1393 e2->u.heaplink = BLI_heap_insert(heap, p_edge_boundary_angle(e2), e2);
1394 }
1395
1396 nedges--;
1397 }
1398 }
1399
1400 BLI_heap_free(heap, NULL);
1401 }
1402
p_chart_fill_boundaries(PChart * chart,PEdge * outer)1403 static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
1404 {
1405 PEdge *e, *be; /* *enext - as yet unused */
1406 int nedges;
1407
1408 for (e = chart->edges; e; e = e->nextlink) {
1409 /* enext = e->nextlink; - as yet unused */
1410
1411 if (e->pair || (e->flag & PEDGE_FILLED)) {
1412 continue;
1413 }
1414
1415 nedges = 0;
1416 be = e;
1417 do {
1418 be->flag |= PEDGE_FILLED;
1419 be = be->next->vert->edge;
1420 nedges++;
1421 } while (be != e);
1422
1423 if (e != outer) {
1424 p_chart_fill_boundary(chart, e, nedges);
1425 }
1426 }
1427 }
1428
1429 #if 0
1430 /* Polygon kernel for inserting uv's non overlapping */
1431
1432 static int p_polygon_point_in(const float cp1[2], const float cp2[2], const float p[2])
1433 {
1434 if ((cp1[0] == p[0]) && (cp1[1] == p[1])) {
1435 return 2;
1436 }
1437 else if ((cp2[0] == p[0]) && (cp2[1] == p[1])) {
1438 return 3;
1439 }
1440 else {
1441 return (p_area_signed(cp1, cp2, p) >= 0.0f);
1442 }
1443 }
1444
1445 static void p_polygon_kernel_clip(float (*oldpoints)[2],
1446 int noldpoints,
1447 float (*newpoints)[2],
1448 int *r_nnewpoints,
1449 const float cp1[2],
1450 const float cp2[2])
1451 {
1452 float *p2, *p1, isect[2];
1453 int i, p2in, p1in;
1454
1455 p1 = oldpoints[noldpoints - 1];
1456 p1in = p_polygon_point_in(cp1, cp2, p1);
1457 *r_nnewpoints = 0;
1458
1459 for (i = 0; i < noldpoints; i++) {
1460 p2 = oldpoints[i];
1461 p2in = p_polygon_point_in(cp1, cp2, p2);
1462
1463 if ((p2in >= 2) || (p1in && p2in)) {
1464 newpoints[*r_nnewpoints][0] = p2[0];
1465 newpoints[*r_nnewpoints][1] = p2[1];
1466 (*r_nnewpoints)++;
1467 }
1468 else if (p1in && !p2in) {
1469 if (p1in != 3) {
1470 p_intersect_line_2d(p1, p2, cp1, cp2, isect);
1471 newpoints[*r_nnewpoints][0] = isect[0];
1472 newpoints[*r_nnewpoints][1] = isect[1];
1473 (*r_nnewpoints)++;
1474 }
1475 }
1476 else if (!p1in && p2in) {
1477 p_intersect_line_2d(p1, p2, cp1, cp2, isect);
1478 newpoints[*r_nnewpoints][0] = isect[0];
1479 newpoints[*r_nnewpoints][1] = isect[1];
1480 (*r_nnewpoints)++;
1481
1482 newpoints[*r_nnewpoints][0] = p2[0];
1483 newpoints[*r_nnewpoints][1] = p2[1];
1484 (*r_nnewpoints)++;
1485 }
1486
1487 p1in = p2in;
1488 p1 = p2;
1489 }
1490 }
1491
1492 static void p_polygon_kernel_center(float (*points)[2], int npoints, float *center)
1493 {
1494 int i, size, nnewpoints = npoints;
1495 float(*oldpoints)[2], (*newpoints)[2], *p1, *p2;
1496
1497 size = npoints * 3;
1498 oldpoints = MEM_mallocN(sizeof(float[2]) * size, "PPolygonOldPoints");
1499 newpoints = MEM_mallocN(sizeof(float[2]) * size, "PPolygonNewPoints");
1500
1501 memcpy(oldpoints, points, sizeof(float[2]) * npoints);
1502
1503 for (i = 0; i < npoints; i++) {
1504 p1 = points[i];
1505 p2 = points[(i + 1) % npoints];
1506 p_polygon_kernel_clip(oldpoints, nnewpoints, newpoints, &nnewpoints, p1, p2);
1507
1508 if (nnewpoints == 0) {
1509 /* degenerate case, use center of original polygon */
1510 memcpy(oldpoints, points, sizeof(float[2]) * npoints);
1511 nnewpoints = npoints;
1512 break;
1513 }
1514 else if (nnewpoints == 1) {
1515 /* degenerate case, use remaining point */
1516 center[0] = newpoints[0][0];
1517 center[1] = newpoints[0][1];
1518
1519 MEM_freeN(oldpoints);
1520 MEM_freeN(newpoints);
1521
1522 return;
1523 }
1524
1525 if (nnewpoints * 2 > size) {
1526 size *= 2;
1527 MEM_freeN(oldpoints);
1528 oldpoints = MEM_mallocN(sizeof(float[2]) * size, "oldpoints");
1529 memcpy(oldpoints, newpoints, sizeof(float[2]) * nnewpoints);
1530 MEM_freeN(newpoints);
1531 newpoints = MEM_mallocN(sizeof(float[2]) * size, "newpoints");
1532 }
1533 else {
1534 float(*sw_points)[2] = oldpoints;
1535 oldpoints = newpoints;
1536 newpoints = sw_points;
1537 }
1538 }
1539
1540 center[0] = center[1] = 0.0f;
1541
1542 for (i = 0; i < nnewpoints; i++) {
1543 center[0] += oldpoints[i][0];
1544 center[1] += oldpoints[i][1];
1545 }
1546
1547 center[0] /= nnewpoints;
1548 center[1] /= nnewpoints;
1549
1550 MEM_freeN(oldpoints);
1551 MEM_freeN(newpoints);
1552 }
1553 #endif
1554
1555 #if 0
1556 /* Edge Collapser */
1557
1558 int NCOLLAPSE = 1;
1559 int NCOLLAPSEX = 0;
1560
1561 static float p_vert_cotan(const float v1[3], const float v2[3], const float v3[3])
1562 {
1563 float a[3], b[3], c[3], clen;
1564
1565 sub_v3_v3v3(a, v2, v1);
1566 sub_v3_v3v3(b, v3, v1);
1567 cross_v3_v3v3(c, a, b);
1568
1569 clen = len_v3(c);
1570
1571 if (clen == 0.0f) {
1572 return 0.0f;
1573 }
1574
1575 return dot_v3v3(a, b) / clen;
1576 }
1577
1578 static PBool p_vert_flipped_wheel_triangle(PVert *v)
1579 {
1580 PEdge *e = v->edge;
1581
1582 do {
1583 if (p_face_uv_area_signed(e->face) < 0.0f) {
1584 return P_TRUE;
1585 }
1586
1587 e = p_wheel_edge_next(e);
1588 } while (e && (e != v->edge));
1589
1590 return P_FALSE;
1591 }
1592
1593 static PBool p_vert_map_harmonic_weights(PVert *v)
1594 {
1595 float weightsum, positionsum[2], olduv[2];
1596
1597 weightsum = 0.0f;
1598 positionsum[0] = positionsum[1] = 0.0f;
1599
1600 if (p_vert_interior(v)) {
1601 PEdge *e = v->edge;
1602
1603 do {
1604 float t1, t2, weight;
1605 PVert *v1, *v2;
1606
1607 v1 = e->next->vert;
1608 v2 = e->next->next->vert;
1609 t1 = p_vert_cotan(v2->co, e->vert->co, v1->co);
1610
1611 v1 = e->pair->next->vert;
1612 v2 = e->pair->next->next->vert;
1613 t2 = p_vert_cotan(v2->co, e->pair->vert->co, v1->co);
1614
1615 weight = 0.5f * (t1 + t2);
1616 weightsum += weight;
1617 positionsum[0] += weight * e->pair->vert->uv[0];
1618 positionsum[1] += weight * e->pair->vert->uv[1];
1619
1620 e = p_wheel_edge_next(e);
1621 } while (e && (e != v->edge));
1622 }
1623 else {
1624 PEdge *e = v->edge;
1625
1626 do {
1627 float t1, t2;
1628 PVert *v1, *v2;
1629
1630 v2 = e->next->vert;
1631 v1 = e->next->next->vert;
1632
1633 t1 = p_vert_cotan(v1->co, v->co, v2->co);
1634 t2 = p_vert_cotan(v2->co, v->co, v1->co);
1635
1636 weightsum += t1 + t2;
1637 positionsum[0] += (v2->uv[1] - v1->uv[1]) + (t1 * v2->uv[0] + t2 * v1->uv[0]);
1638 positionsum[1] += (v1->uv[0] - v2->uv[0]) + (t1 * v2->uv[1] + t2 * v1->uv[1]);
1639
1640 e = p_wheel_edge_next(e);
1641 } while (e && (e != v->edge));
1642 }
1643
1644 if (weightsum != 0.0f) {
1645 weightsum = 1.0f / weightsum;
1646 positionsum[0] *= weightsum;
1647 positionsum[1] *= weightsum;
1648 }
1649
1650 olduv[0] = v->uv[0];
1651 olduv[1] = v->uv[1];
1652 v->uv[0] = positionsum[0];
1653 v->uv[1] = positionsum[1];
1654
1655 if (p_vert_flipped_wheel_triangle(v)) {
1656 v->uv[0] = olduv[0];
1657 v->uv[1] = olduv[1];
1658
1659 return P_FALSE;
1660 }
1661
1662 return P_TRUE;
1663 }
1664
1665 static void p_vert_harmonic_insert(PVert *v)
1666 {
1667 PEdge *e;
1668
1669 if (!p_vert_map_harmonic_weights(v)) {
1670 /* do polygon kernel center insertion: this is quite slow, but should
1671 * only be needed for 0.01 % of verts or so, when insert with harmonic
1672 * weights fails */
1673
1674 int npoints = 0, i;
1675 float(*points)[2];
1676
1677 e = v->edge;
1678 do {
1679 npoints++;
1680 e = p_wheel_edge_next(e);
1681 } while (e && (e != v->edge));
1682
1683 if (e == NULL) {
1684 npoints++;
1685 }
1686
1687 points = MEM_mallocN(sizeof(float[2]) * npoints, "PHarmonicPoints");
1688
1689 e = v->edge;
1690 i = 0;
1691 do {
1692 PEdge *nexte = p_wheel_edge_next(e);
1693
1694 points[i][0] = e->next->vert->uv[0];
1695 points[i][1] = e->next->vert->uv[1];
1696
1697 if (nexte == NULL) {
1698 i++;
1699 points[i][0] = e->next->next->vert->uv[0];
1700 points[i][1] = e->next->next->vert->uv[1];
1701 break;
1702 }
1703
1704 e = nexte;
1705 i++;
1706 } while (e != v->edge);
1707
1708 p_polygon_kernel_center(points, npoints, v->uv);
1709
1710 MEM_freeN(points);
1711 }
1712
1713 e = v->edge;
1714 do {
1715 if (!(e->next->vert->flag & PVERT_PIN)) {
1716 p_vert_map_harmonic_weights(e->next->vert);
1717 }
1718 e = p_wheel_edge_next(e);
1719 } while (e && (e != v->edge));
1720
1721 p_vert_map_harmonic_weights(v);
1722 }
1723
1724 static void p_vert_fix_edge_pointer(PVert *v)
1725 {
1726 PEdge *start = v->edge;
1727
1728 /* set v->edge pointer to the edge with no pair, if there is one */
1729 while (v->edge->pair) {
1730 v->edge = p_wheel_edge_prev(v->edge);
1731
1732 if (v->edge == start) {
1733 break;
1734 }
1735 }
1736 }
1737
1738 static void p_collapsing_verts(PEdge *edge, PEdge *pair, PVert **r_newv, PVert **r_keepv)
1739 {
1740 /* the two vertices that are involved in the collapse */
1741 if (edge) {
1742 *r_newv = edge->vert;
1743 *r_keepv = edge->next->vert;
1744 }
1745 else {
1746 *r_newv = pair->next->vert;
1747 *r_keepv = pair->vert;
1748 }
1749 }
1750
1751 static void p_collapse_edge(PEdge *edge, PEdge *pair)
1752 {
1753 PVert *oldv, *keepv;
1754 PEdge *e;
1755
1756 p_collapsing_verts(edge, pair, &oldv, &keepv);
1757
1758 /* change e->vert pointers from old vertex to the target vertex */
1759 e = oldv->edge;
1760 do {
1761 if ((e != edge) && !(pair && pair->next == e)) {
1762 e->vert = keepv;
1763 }
1764
1765 e = p_wheel_edge_next(e);
1766 } while (e && (e != oldv->edge));
1767
1768 /* set keepv->edge pointer */
1769 if ((edge && (keepv->edge == edge->next)) || (keepv->edge == pair)) {
1770 if (edge && edge->next->pair) {
1771 keepv->edge = edge->next->pair->next;
1772 }
1773 else if (pair && pair->next->next->pair) {
1774 keepv->edge = pair->next->next->pair;
1775 }
1776 else if (edge && edge->next->next->pair) {
1777 keepv->edge = edge->next->next->pair;
1778 }
1779 else {
1780 keepv->edge = pair->next->pair->next;
1781 }
1782 }
1783
1784 /* update pairs and v->edge pointers */
1785 if (edge) {
1786 PEdge *e1 = edge->next, *e2 = e1->next;
1787
1788 if (e1->pair) {
1789 e1->pair->pair = e2->pair;
1790 }
1791
1792 if (e2->pair) {
1793 e2->pair->pair = e1->pair;
1794 e2->vert->edge = p_wheel_edge_prev(e2);
1795 }
1796 else {
1797 e2->vert->edge = p_wheel_edge_next(e2);
1798 }
1799
1800 p_vert_fix_edge_pointer(e2->vert);
1801 }
1802
1803 if (pair) {
1804 PEdge *e1 = pair->next, *e2 = e1->next;
1805
1806 if (e1->pair) {
1807 e1->pair->pair = e2->pair;
1808 }
1809
1810 if (e2->pair) {
1811 e2->pair->pair = e1->pair;
1812 e2->vert->edge = p_wheel_edge_prev(e2);
1813 }
1814 else {
1815 e2->vert->edge = p_wheel_edge_next(e2);
1816 }
1817
1818 p_vert_fix_edge_pointer(e2->vert);
1819 }
1820
1821 p_vert_fix_edge_pointer(keepv);
1822
1823 /* mark for move to collapsed list later */
1824 oldv->flag |= PVERT_COLLAPSE;
1825
1826 if (edge) {
1827 PFace *f = edge->face;
1828 PEdge *e1 = edge->next, *e2 = e1->next;
1829
1830 f->flag |= PFACE_COLLAPSE;
1831 edge->flag |= PEDGE_COLLAPSE;
1832 e1->flag |= PEDGE_COLLAPSE;
1833 e2->flag |= PEDGE_COLLAPSE;
1834 }
1835
1836 if (pair) {
1837 PFace *f = pair->face;
1838 PEdge *e1 = pair->next, *e2 = e1->next;
1839
1840 f->flag |= PFACE_COLLAPSE;
1841 pair->flag |= PEDGE_COLLAPSE;
1842 e1->flag |= PEDGE_COLLAPSE;
1843 e2->flag |= PEDGE_COLLAPSE;
1844 }
1845 }
1846
1847 static void p_split_vertex(PEdge *edge, PEdge *pair)
1848 {
1849 PVert *newv, *keepv;
1850 PEdge *e;
1851
1852 p_collapsing_verts(edge, pair, &newv, &keepv);
1853
1854 /* update edge pairs */
1855 if (edge) {
1856 PEdge *e1 = edge->next, *e2 = e1->next;
1857
1858 if (e1->pair) {
1859 e1->pair->pair = e1;
1860 }
1861 if (e2->pair) {
1862 e2->pair->pair = e2;
1863 }
1864
1865 e2->vert->edge = e2;
1866 p_vert_fix_edge_pointer(e2->vert);
1867 keepv->edge = e1;
1868 }
1869
1870 if (pair) {
1871 PEdge *e1 = pair->next, *e2 = e1->next;
1872
1873 if (e1->pair) {
1874 e1->pair->pair = e1;
1875 }
1876 if (e2->pair) {
1877 e2->pair->pair = e2;
1878 }
1879
1880 e2->vert->edge = e2;
1881 p_vert_fix_edge_pointer(e2->vert);
1882 keepv->edge = pair;
1883 }
1884
1885 p_vert_fix_edge_pointer(keepv);
1886
1887 /* set e->vert pointers to restored vertex */
1888 e = newv->edge;
1889 do {
1890 e->vert = newv;
1891 e = p_wheel_edge_next(e);
1892 } while (e && (e != newv->edge));
1893 }
1894
1895 static PBool p_collapse_allowed_topologic(PEdge *edge, PEdge *pair)
1896 {
1897 PVert *oldv, *keepv;
1898
1899 p_collapsing_verts(edge, pair, &oldv, &keepv);
1900
1901 /* boundary edges */
1902 if (!edge || !pair) {
1903 /* avoid collapsing chart into an edge */
1904 if (edge && !edge->next->pair && !edge->next->next->pair) {
1905 return P_FALSE;
1906 }
1907 else if (pair && !pair->next->pair && !pair->next->next->pair) {
1908 return P_FALSE;
1909 }
1910 }
1911 /* avoid merging two boundaries (oldv and keepv are on the 'other side' of
1912 * the chart) */
1913 else if (!p_vert_interior(oldv) && !p_vert_interior(keepv)) {
1914 return P_FALSE;
1915 }
1916
1917 return P_TRUE;
1918 }
1919
1920 static PBool p_collapse_normal_flipped(float *v1, float *v2, float *vold, float *vnew)
1921 {
1922 float nold[3], nnew[3], sub1[3], sub2[3];
1923
1924 sub_v3_v3v3(sub1, vold, v1);
1925 sub_v3_v3v3(sub2, vold, v2);
1926 cross_v3_v3v3(nold, sub1, sub2);
1927
1928 sub_v3_v3v3(sub1, vnew, v1);
1929 sub_v3_v3v3(sub2, vnew, v2);
1930 cross_v3_v3v3(nnew, sub1, sub2);
1931
1932 return (dot_v3v3(nold, nnew) <= 0.0f);
1933 }
1934
1935 static PBool p_collapse_allowed_geometric(PEdge *edge, PEdge *pair)
1936 {
1937 PVert *oldv, *keepv;
1938 PEdge *e;
1939 float angulardefect, angle;
1940
1941 p_collapsing_verts(edge, pair, &oldv, &keepv);
1942
1943 angulardefect = 2 * M_PI;
1944
1945 e = oldv->edge;
1946 do {
1947 float a[3], b[3], minangle, maxangle;
1948 PEdge *e1 = e->next, *e2 = e1->next;
1949 PVert *v1 = e1->vert, *v2 = e2->vert;
1950 int i;
1951
1952 angle = p_vec_angle(v1->co, oldv->co, v2->co);
1953 angulardefect -= angle;
1954
1955 /* skip collapsing faces */
1956 if (v1 == keepv || v2 == keepv) {
1957 e = p_wheel_edge_next(e);
1958 continue;
1959 }
1960
1961 if (p_collapse_normal_flipped(v1->co, v2->co, oldv->co, keepv->co)) {
1962 return P_FALSE;
1963 }
1964
1965 a[0] = angle;
1966 a[1] = p_vec_angle(v2->co, v1->co, oldv->co);
1967 a[2] = M_PI - a[0] - a[1];
1968
1969 b[0] = p_vec_angle(v1->co, keepv->co, v2->co);
1970 b[1] = p_vec_angle(v2->co, v1->co, keepv->co);
1971 b[2] = M_PI - b[0] - b[1];
1972
1973 /* abf criterion 1: avoid sharp and obtuse angles */
1974 minangle = 15.0f * M_PI / 180.0f;
1975 maxangle = M_PI - minangle;
1976
1977 for (i = 0; i < 3; i++) {
1978 if ((b[i] < a[i]) && (b[i] < minangle)) {
1979 return P_FALSE;
1980 }
1981 else if ((b[i] > a[i]) && (b[i] > maxangle)) {
1982 return P_FALSE;
1983 }
1984 }
1985
1986 e = p_wheel_edge_next(e);
1987 } while (e && (e != oldv->edge));
1988
1989 if (p_vert_interior(oldv)) {
1990 /* hlscm criterion: angular defect smaller than threshold */
1991 if (fabsf(angulardefect) > (float)(M_PI * 30.0 / 180.0)) {
1992 return P_FALSE;
1993 }
1994 }
1995 else {
1996 PVert *v1 = p_boundary_edge_next(oldv->edge)->vert;
1997 PVert *v2 = p_boundary_edge_prev(oldv->edge)->vert;
1998
1999 /* abf++ criterion 2: avoid collapsing verts inwards */
2000 if (p_vert_interior(keepv)) {
2001 return P_FALSE;
2002 }
2003
2004 /* don't collapse significant boundary changes */
2005 angle = p_vec_angle(v1->co, oldv->co, v2->co);
2006 if (angle < (M_PI * 160.0 / 180.0)) {
2007 return P_FALSE;
2008 }
2009 }
2010
2011 return P_TRUE;
2012 }
2013
2014 static PBool p_collapse_allowed(PEdge *edge, PEdge *pair)
2015 {
2016 PVert *oldv, *keepv;
2017
2018 p_collapsing_verts(edge, pair, &oldv, &keepv);
2019
2020 if (oldv->flag & PVERT_PIN) {
2021 return P_FALSE;
2022 }
2023
2024 return (p_collapse_allowed_topologic(edge, pair) && p_collapse_allowed_geometric(edge, pair));
2025 }
2026
2027 static float p_collapse_cost(PEdge *edge, PEdge *pair)
2028 {
2029 /* based on volume and boundary optimization from:
2030 * "Fast and Memory Efficient Polygonal Simplification" P. Lindstrom, G. Turk */
2031
2032 PVert *oldv, *keepv;
2033 PEdge *e;
2034 PFace *oldf1, *oldf2;
2035 float volumecost = 0.0f, areacost = 0.0f, edgevec[3], cost, weight, elen;
2036 float shapecost = 0.0f;
2037 float shapeold = 0.0f, shapenew = 0.0f;
2038 int nshapeold = 0, nshapenew = 0;
2039
2040 p_collapsing_verts(edge, pair, &oldv, &keepv);
2041 oldf1 = (edge) ? edge->face : NULL;
2042 oldf2 = (pair) ? pair->face : NULL;
2043
2044 sub_v3_v3v3(edgevec, keepv->co, oldv->co);
2045
2046 e = oldv->edge;
2047 do {
2048 float a1, a2, a3;
2049 float *co1 = e->next->vert->co;
2050 float *co2 = e->next->next->vert->co;
2051
2052 if ((e->face != oldf1) && (e->face != oldf2)) {
2053 float tetrav2[3], tetrav3[3];
2054
2055 /* tetrahedron volume = (1/3!)*|a.(b x c)| */
2056 sub_v3_v3v3(tetrav2, co1, oldv->co);
2057 sub_v3_v3v3(tetrav3, co2, oldv->co);
2058 volumecost += fabsf(volume_tri_tetrahedron_signed_v3(tetrav2, tetrav3, edgevec));
2059
2060 # if 0
2061 shapecost += dot_v3v3(co1, keepv->co);
2062
2063 if (p_wheel_edge_next(e) == NULL) {
2064 shapecost += dot_v3v3(co2, keepv->co);
2065 }
2066 # endif
2067
2068 p_triangle_angles(oldv->co, co1, co2, &a1, &a2, &a3);
2069 a1 = a1 - M_PI / 3.0;
2070 a2 = a2 - M_PI / 3.0;
2071 a3 = a3 - M_PI / 3.0;
2072 shapeold = (a1 * a1 + a2 * a2 + a3 * a3) / ((M_PI / 2) * (M_PI / 2));
2073
2074 nshapeold++;
2075 }
2076 else {
2077 p_triangle_angles(keepv->co, co1, co2, &a1, &a2, &a3);
2078 a1 = a1 - M_PI / 3.0;
2079 a2 = a2 - M_PI / 3.0;
2080 a3 = a3 - M_PI / 3.0;
2081 shapenew = (a1 * a1 + a2 * a2 + a3 * a3) / ((M_PI / 2) * (M_PI / 2));
2082
2083 nshapenew++;
2084 }
2085
2086 e = p_wheel_edge_next(e);
2087 } while (e && (e != oldv->edge));
2088
2089 if (!p_vert_interior(oldv)) {
2090 PVert *v1 = p_boundary_edge_prev(oldv->edge)->vert;
2091 PVert *v2 = p_boundary_edge_next(oldv->edge)->vert;
2092
2093 areacost = area_tri_v3(oldv->co, v1->co, v2->co);
2094 }
2095
2096 elen = len_v3(edgevec);
2097 weight = 1.0f; /* 0.2f */
2098 cost = weight * volumecost * volumecost + elen * elen * areacost * areacost;
2099 # if 0
2100 cost += shapecost;
2101 # else
2102 shapeold /= nshapeold;
2103 shapenew /= nshapenew;
2104 shapecost = (shapeold + 0.00001) / (shapenew + 0.00001);
2105
2106 cost *= shapecost;
2107 # endif
2108
2109 return cost;
2110 }
2111
2112 static void p_collapse_cost_vertex(PVert *vert, float *r_mincost, PEdge **r_mine)
2113 {
2114 PEdge *e, *enext, *pair;
2115
2116 *r_mine = NULL;
2117 *r_mincost = 0.0f;
2118 e = vert->edge;
2119 do {
2120 if (p_collapse_allowed(e, e->pair)) {
2121 float cost = p_collapse_cost(e, e->pair);
2122
2123 if ((*r_mine == NULL) || (cost < *r_mincost)) {
2124 *r_mincost = cost;
2125 *r_mine = e;
2126 }
2127 }
2128
2129 enext = p_wheel_edge_next(e);
2130
2131 if (enext == NULL) {
2132 /* the other boundary edge, where we only have the pair halfedge */
2133 pair = e->next->next;
2134
2135 if (p_collapse_allowed(NULL, pair)) {
2136 float cost = p_collapse_cost(NULL, pair);
2137
2138 if ((*r_mine == NULL) || (cost < *r_mincost)) {
2139 *r_mincost = cost;
2140 *r_mine = pair;
2141 }
2142 }
2143
2144 break;
2145 }
2146
2147 e = enext;
2148 } while (e != vert->edge);
2149 }
2150
2151 static void p_chart_post_collapse_flush(PChart *chart, PEdge *collapsed)
2152 {
2153 /* move to collapsed_ */
2154
2155 PVert *v, *nextv = NULL, *verts = chart->verts;
2156 PEdge *e, *nexte = NULL, *edges = chart->edges, *laste = NULL;
2157 PFace *f, *nextf = NULL, *faces = chart->faces;
2158
2159 chart->verts = chart->collapsed_verts = NULL;
2160 chart->edges = chart->collapsed_edges = NULL;
2161 chart->faces = chart->collapsed_faces = NULL;
2162
2163 chart->nverts = chart->nedges = chart->nfaces = 0;
2164
2165 for (v = verts; v; v = nextv) {
2166 nextv = v->nextlink;
2167
2168 if (v->flag & PVERT_COLLAPSE) {
2169 v->nextlink = chart->collapsed_verts;
2170 chart->collapsed_verts = v;
2171 }
2172 else {
2173 v->nextlink = chart->verts;
2174 chart->verts = v;
2175 chart->nverts++;
2176 }
2177 }
2178
2179 for (e = edges; e; e = nexte) {
2180 nexte = e->nextlink;
2181
2182 if (!collapsed || !(e->flag & PEDGE_COLLAPSE_EDGE)) {
2183 if (e->flag & PEDGE_COLLAPSE) {
2184 e->nextlink = chart->collapsed_edges;
2185 chart->collapsed_edges = e;
2186 }
2187 else {
2188 e->nextlink = chart->edges;
2189 chart->edges = e;
2190 chart->nedges++;
2191 }
2192 }
2193 }
2194
2195 /* these are added last so they can be popped of in the right order
2196 * for splitting */
2197 for (e = collapsed; e; e = e->nextlink) {
2198 e->nextlink = e->u.nextcollapse;
2199 laste = e;
2200 }
2201 if (laste) {
2202 laste->nextlink = chart->collapsed_edges;
2203 chart->collapsed_edges = collapsed;
2204 }
2205
2206 for (f = faces; f; f = nextf) {
2207 nextf = f->nextlink;
2208
2209 if (f->flag & PFACE_COLLAPSE) {
2210 f->nextlink = chart->collapsed_faces;
2211 chart->collapsed_faces = f;
2212 }
2213 else {
2214 f->nextlink = chart->faces;
2215 chart->faces = f;
2216 chart->nfaces++;
2217 }
2218 }
2219 }
2220
2221 static void p_chart_post_split_flush(PChart *chart)
2222 {
2223 /* move from collapsed_ */
2224
2225 PVert *v, *nextv = NULL;
2226 PEdge *e, *nexte = NULL;
2227 PFace *f, *nextf = NULL;
2228
2229 for (v = chart->collapsed_verts; v; v = nextv) {
2230 nextv = v->nextlink;
2231 v->nextlink = chart->verts;
2232 chart->verts = v;
2233 chart->nverts++;
2234 }
2235
2236 for (e = chart->collapsed_edges; e; e = nexte) {
2237 nexte = e->nextlink;
2238 e->nextlink = chart->edges;
2239 chart->edges = e;
2240 chart->nedges++;
2241 }
2242
2243 for (f = chart->collapsed_faces; f; f = nextf) {
2244 nextf = f->nextlink;
2245 f->nextlink = chart->faces;
2246 chart->faces = f;
2247 chart->nfaces++;
2248 }
2249
2250 chart->collapsed_verts = NULL;
2251 chart->collapsed_edges = NULL;
2252 chart->collapsed_faces = NULL;
2253 }
2254
2255 static void p_chart_simplify_compute(PChart *chart)
2256 {
2257 /* Computes a list of edge collapses / vertex splits. The collapsed
2258 * simplices go in the chart->collapsed_* lists, The original and
2259 * collapsed may then be view as stacks, where the next collapse/split
2260 * is at the top of the respective lists. */
2261
2262 Heap *heap = BLI_heap_new();
2263 PVert *v, **wheelverts;
2264 PEdge *collapsededges = NULL, *e;
2265 int nwheelverts, i, ncollapsed = 0;
2266
2267 wheelverts = MEM_mallocN(sizeof(PVert *) * chart->nverts, "PChartWheelVerts");
2268
2269 /* insert all potential collapses into heap */
2270 for (v = chart->verts; v; v = v->nextlink) {
2271 float cost;
2272 PEdge *e = NULL;
2273
2274 p_collapse_cost_vertex(v, &cost, &e);
2275
2276 if (e) {
2277 v->u.heaplink = BLI_heap_insert(heap, cost, e);
2278 }
2279 else {
2280 v->u.heaplink = NULL;
2281 }
2282 }
2283
2284 for (e = chart->edges; e; e = e->nextlink) {
2285 e->u.nextcollapse = NULL;
2286 }
2287
2288 /* pop edge collapse out of heap one by one */
2289 while (!BLI_heap_is_empty(heap)) {
2290 if (ncollapsed == NCOLLAPSE) {
2291 break;
2292 }
2293
2294 HeapNode *link = BLI_heap_top(heap);
2295 PEdge *edge = (PEdge *)BLI_heap_pop_min(heap), *pair = edge->pair;
2296 PVert *oldv, *keepv;
2297 PEdge *wheele, *nexte;
2298
2299 /* remember the edges we collapsed */
2300 edge->u.nextcollapse = collapsededges;
2301 collapsededges = edge;
2302
2303 if (edge->vert->u.heaplink != link) {
2304 edge->flag |= (PEDGE_COLLAPSE_EDGE | PEDGE_COLLAPSE_PAIR);
2305 edge->next->vert->u.heaplink = NULL;
2306 SWAP(PEdge *, edge, pair);
2307 }
2308 else {
2309 edge->flag |= PEDGE_COLLAPSE_EDGE;
2310 edge->vert->u.heaplink = NULL;
2311 }
2312
2313 p_collapsing_verts(edge, pair, &oldv, &keepv);
2314
2315 /* gather all wheel verts and remember them before collapse */
2316 nwheelverts = 0;
2317 wheele = oldv->edge;
2318
2319 do {
2320 wheelverts[nwheelverts++] = wheele->next->vert;
2321 nexte = p_wheel_edge_next(wheele);
2322
2323 if (nexte == NULL) {
2324 wheelverts[nwheelverts++] = wheele->next->next->vert;
2325 }
2326
2327 wheele = nexte;
2328 } while (wheele && (wheele != oldv->edge));
2329
2330 /* collapse */
2331 p_collapse_edge(edge, pair);
2332
2333 for (i = 0; i < nwheelverts; i++) {
2334 float cost;
2335 PEdge *collapse = NULL;
2336
2337 v = wheelverts[i];
2338
2339 if (v->u.heaplink) {
2340 BLI_heap_remove(heap, v->u.heaplink);
2341 v->u.heaplink = NULL;
2342 }
2343
2344 p_collapse_cost_vertex(v, &cost, &collapse);
2345
2346 if (collapse) {
2347 v->u.heaplink = BLI_heap_insert(heap, cost, collapse);
2348 }
2349 }
2350
2351 ncollapsed++;
2352 }
2353
2354 MEM_freeN(wheelverts);
2355 BLI_heap_free(heap, NULL);
2356
2357 p_chart_post_collapse_flush(chart, collapsededges);
2358 }
2359
2360 static void p_chart_complexify(PChart *chart)
2361 {
2362 PEdge *e, *pair, *edge;
2363 PVert *newv, *keepv;
2364 int x = 0;
2365
2366 for (e = chart->collapsed_edges; e; e = e->nextlink) {
2367 if (!(e->flag & PEDGE_COLLAPSE_EDGE)) {
2368 break;
2369 }
2370
2371 edge = e;
2372 pair = e->pair;
2373
2374 if (edge->flag & PEDGE_COLLAPSE_PAIR) {
2375 SWAP(PEdge *, edge, pair);
2376 }
2377
2378 p_split_vertex(edge, pair);
2379 p_collapsing_verts(edge, pair, &newv, &keepv);
2380
2381 if (x >= NCOLLAPSEX) {
2382 newv->uv[0] = keepv->uv[0];
2383 newv->uv[1] = keepv->uv[1];
2384 }
2385 else {
2386 p_vert_harmonic_insert(newv);
2387 x++;
2388 }
2389 }
2390
2391 p_chart_post_split_flush(chart);
2392 }
2393
2394 # if 0
2395 static void p_chart_simplify(PChart *chart)
2396 {
2397 /* Not implemented, needs proper reordering in split_flush. */
2398 }
2399 # endif
2400 #endif
2401
2402 /* ABF */
2403
2404 #define ABF_MAX_ITER 20
2405
2406 typedef struct PAbfSystem {
2407 int ninterior, nfaces, nangles;
2408 float *alpha, *beta, *sine, *cosine, *weight;
2409 float *bAlpha, *bTriangle, *bInterior;
2410 float *lambdaTriangle, *lambdaPlanar, *lambdaLength;
2411 float (*J2dt)[3], *bstar, *dstar;
2412 float minangle, maxangle;
2413 } PAbfSystem;
2414
p_abf_setup_system(PAbfSystem * sys)2415 static void p_abf_setup_system(PAbfSystem *sys)
2416 {
2417 int i;
2418
2419 sys->alpha = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFalpha");
2420 sys->beta = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFbeta");
2421 sys->sine = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFsine");
2422 sys->cosine = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFcosine");
2423 sys->weight = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFweight");
2424
2425 sys->bAlpha = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFbalpha");
2426 sys->bTriangle = (float *)MEM_mallocN(sizeof(float) * sys->nfaces, "ABFbtriangle");
2427 sys->bInterior = (float *)MEM_mallocN(sizeof(float[2]) * sys->ninterior, "ABFbinterior");
2428
2429 sys->lambdaTriangle = (float *)MEM_callocN(sizeof(float) * sys->nfaces, "ABFlambdatri");
2430 sys->lambdaPlanar = (float *)MEM_callocN(sizeof(float) * sys->ninterior, "ABFlamdaplane");
2431 sys->lambdaLength = (float *)MEM_mallocN(sizeof(float) * sys->ninterior, "ABFlambdalen");
2432
2433 sys->J2dt = MEM_mallocN(sizeof(float) * sys->nangles * 3, "ABFj2dt");
2434 sys->bstar = (float *)MEM_mallocN(sizeof(float) * sys->nfaces, "ABFbstar");
2435 sys->dstar = (float *)MEM_mallocN(sizeof(float) * sys->nfaces, "ABFdstar");
2436
2437 for (i = 0; i < sys->ninterior; i++) {
2438 sys->lambdaLength[i] = 1.0;
2439 }
2440
2441 sys->minangle = 1.0 * M_PI / 180.0;
2442 sys->maxangle = (float)M_PI - sys->minangle;
2443 }
2444
p_abf_free_system(PAbfSystem * sys)2445 static void p_abf_free_system(PAbfSystem *sys)
2446 {
2447 MEM_freeN(sys->alpha);
2448 MEM_freeN(sys->beta);
2449 MEM_freeN(sys->sine);
2450 MEM_freeN(sys->cosine);
2451 MEM_freeN(sys->weight);
2452 MEM_freeN(sys->bAlpha);
2453 MEM_freeN(sys->bTriangle);
2454 MEM_freeN(sys->bInterior);
2455 MEM_freeN(sys->lambdaTriangle);
2456 MEM_freeN(sys->lambdaPlanar);
2457 MEM_freeN(sys->lambdaLength);
2458 MEM_freeN(sys->J2dt);
2459 MEM_freeN(sys->bstar);
2460 MEM_freeN(sys->dstar);
2461 }
2462
p_abf_compute_sines(PAbfSystem * sys)2463 static void p_abf_compute_sines(PAbfSystem *sys)
2464 {
2465 int i;
2466 float *sine = sys->sine, *cosine = sys->cosine, *alpha = sys->alpha;
2467
2468 for (i = 0; i < sys->nangles; i++, sine++, cosine++, alpha++) {
2469 *sine = sinf(*alpha);
2470 *cosine = cosf(*alpha);
2471 }
2472 }
2473
p_abf_compute_sin_product(PAbfSystem * sys,PVert * v,int aid)2474 static float p_abf_compute_sin_product(PAbfSystem *sys, PVert *v, int aid)
2475 {
2476 PEdge *e, *e1, *e2;
2477 float sin1, sin2;
2478
2479 sin1 = sin2 = 1.0;
2480
2481 e = v->edge;
2482 do {
2483 e1 = e->next;
2484 e2 = e->next->next;
2485
2486 if (aid == e1->u.id) {
2487 /* we are computing a derivative for this angle,
2488 * so we use cos and drop the other part */
2489 sin1 *= sys->cosine[e1->u.id];
2490 sin2 = 0.0;
2491 }
2492 else {
2493 sin1 *= sys->sine[e1->u.id];
2494 }
2495
2496 if (aid == e2->u.id) {
2497 /* see above */
2498 sin1 = 0.0;
2499 sin2 *= sys->cosine[e2->u.id];
2500 }
2501 else {
2502 sin2 *= sys->sine[e2->u.id];
2503 }
2504
2505 e = e->next->next->pair;
2506 } while (e && (e != v->edge));
2507
2508 return (sin1 - sin2);
2509 }
2510
p_abf_compute_grad_alpha(PAbfSystem * sys,PFace * f,PEdge * e)2511 static float p_abf_compute_grad_alpha(PAbfSystem *sys, PFace *f, PEdge *e)
2512 {
2513 PVert *v = e->vert, *v1 = e->next->vert, *v2 = e->next->next->vert;
2514 float deriv;
2515
2516 deriv = (sys->alpha[e->u.id] - sys->beta[e->u.id]) * sys->weight[e->u.id];
2517 deriv += sys->lambdaTriangle[f->u.id];
2518
2519 if (v->flag & PVERT_INTERIOR) {
2520 deriv += sys->lambdaPlanar[v->u.id];
2521 }
2522
2523 if (v1->flag & PVERT_INTERIOR) {
2524 float product = p_abf_compute_sin_product(sys, v1, e->u.id);
2525 deriv += sys->lambdaLength[v1->u.id] * product;
2526 }
2527
2528 if (v2->flag & PVERT_INTERIOR) {
2529 float product = p_abf_compute_sin_product(sys, v2, e->u.id);
2530 deriv += sys->lambdaLength[v2->u.id] * product;
2531 }
2532
2533 return deriv;
2534 }
2535
p_abf_compute_gradient(PAbfSystem * sys,PChart * chart)2536 static float p_abf_compute_gradient(PAbfSystem *sys, PChart *chart)
2537 {
2538 PFace *f;
2539 PEdge *e;
2540 PVert *v;
2541 float norm = 0.0;
2542
2543 for (f = chart->faces; f; f = f->nextlink) {
2544 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
2545 float gtriangle, galpha1, galpha2, galpha3;
2546
2547 galpha1 = p_abf_compute_grad_alpha(sys, f, e1);
2548 galpha2 = p_abf_compute_grad_alpha(sys, f, e2);
2549 galpha3 = p_abf_compute_grad_alpha(sys, f, e3);
2550
2551 sys->bAlpha[e1->u.id] = -galpha1;
2552 sys->bAlpha[e2->u.id] = -galpha2;
2553 sys->bAlpha[e3->u.id] = -galpha3;
2554
2555 norm += galpha1 * galpha1 + galpha2 * galpha2 + galpha3 * galpha3;
2556
2557 gtriangle = sys->alpha[e1->u.id] + sys->alpha[e2->u.id] + sys->alpha[e3->u.id] - (float)M_PI;
2558 sys->bTriangle[f->u.id] = -gtriangle;
2559 norm += gtriangle * gtriangle;
2560 }
2561
2562 for (v = chart->verts; v; v = v->nextlink) {
2563 if (v->flag & PVERT_INTERIOR) {
2564 float gplanar = -2 * M_PI, glength;
2565
2566 e = v->edge;
2567 do {
2568 gplanar += sys->alpha[e->u.id];
2569 e = e->next->next->pair;
2570 } while (e && (e != v->edge));
2571
2572 sys->bInterior[v->u.id] = -gplanar;
2573 norm += gplanar * gplanar;
2574
2575 glength = p_abf_compute_sin_product(sys, v, -1);
2576 sys->bInterior[sys->ninterior + v->u.id] = -glength;
2577 norm += glength * glength;
2578 }
2579 }
2580
2581 return norm;
2582 }
2583
p_abf_matrix_invert(PAbfSystem * sys,PChart * chart)2584 static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
2585 {
2586 PFace *f;
2587 PEdge *e;
2588 int i, j, ninterior = sys->ninterior, nvar = 2 * sys->ninterior;
2589 PBool success;
2590 LinearSolver *context;
2591
2592 context = EIG_linear_solver_new(0, nvar, 1);
2593
2594 for (i = 0; i < nvar; i++) {
2595 EIG_linear_solver_right_hand_side_add(context, 0, i, sys->bInterior[i]);
2596 }
2597
2598 for (f = chart->faces; f; f = f->nextlink) {
2599 float wi1, wi2, wi3, b, si, beta[3], j2[3][3], W[3][3];
2600 float row1[6], row2[6], row3[6];
2601 int vid[6];
2602 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
2603 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
2604
2605 wi1 = 1.0f / sys->weight[e1->u.id];
2606 wi2 = 1.0f / sys->weight[e2->u.id];
2607 wi3 = 1.0f / sys->weight[e3->u.id];
2608
2609 /* bstar1 = (J1*dInv*bAlpha - bTriangle) */
2610 b = sys->bAlpha[e1->u.id] * wi1;
2611 b += sys->bAlpha[e2->u.id] * wi2;
2612 b += sys->bAlpha[e3->u.id] * wi3;
2613 b -= sys->bTriangle[f->u.id];
2614
2615 /* si = J1*d*J1t */
2616 si = 1.0f / (wi1 + wi2 + wi3);
2617
2618 /* J1t*si*bstar1 - bAlpha */
2619 beta[0] = b * si - sys->bAlpha[e1->u.id];
2620 beta[1] = b * si - sys->bAlpha[e2->u.id];
2621 beta[2] = b * si - sys->bAlpha[e3->u.id];
2622
2623 /* use this later for computing other lambda's */
2624 sys->bstar[f->u.id] = b;
2625 sys->dstar[f->u.id] = si;
2626
2627 /* set matrix */
2628 W[0][0] = si - sys->weight[e1->u.id];
2629 W[0][1] = si;
2630 W[0][2] = si;
2631 W[1][0] = si;
2632 W[1][1] = si - sys->weight[e2->u.id];
2633 W[1][2] = si;
2634 W[2][0] = si;
2635 W[2][1] = si;
2636 W[2][2] = si - sys->weight[e3->u.id];
2637
2638 vid[0] = vid[1] = vid[2] = vid[3] = vid[4] = vid[5] = -1;
2639
2640 if (v1->flag & PVERT_INTERIOR) {
2641 vid[0] = v1->u.id;
2642 vid[3] = ninterior + v1->u.id;
2643
2644 sys->J2dt[e1->u.id][0] = j2[0][0] = 1.0f * wi1;
2645 sys->J2dt[e2->u.id][0] = j2[1][0] = p_abf_compute_sin_product(sys, v1, e2->u.id) * wi2;
2646 sys->J2dt[e3->u.id][0] = j2[2][0] = p_abf_compute_sin_product(sys, v1, e3->u.id) * wi3;
2647
2648 EIG_linear_solver_right_hand_side_add(context, 0, v1->u.id, j2[0][0] * beta[0]);
2649 EIG_linear_solver_right_hand_side_add(
2650 context, 0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
2651
2652 row1[0] = j2[0][0] * W[0][0];
2653 row2[0] = j2[0][0] * W[1][0];
2654 row3[0] = j2[0][0] * W[2][0];
2655
2656 row1[3] = j2[1][0] * W[0][1] + j2[2][0] * W[0][2];
2657 row2[3] = j2[1][0] * W[1][1] + j2[2][0] * W[1][2];
2658 row3[3] = j2[1][0] * W[2][1] + j2[2][0] * W[2][2];
2659 }
2660
2661 if (v2->flag & PVERT_INTERIOR) {
2662 vid[1] = v2->u.id;
2663 vid[4] = ninterior + v2->u.id;
2664
2665 sys->J2dt[e1->u.id][1] = j2[0][1] = p_abf_compute_sin_product(sys, v2, e1->u.id) * wi1;
2666 sys->J2dt[e2->u.id][1] = j2[1][1] = 1.0f * wi2;
2667 sys->J2dt[e3->u.id][1] = j2[2][1] = p_abf_compute_sin_product(sys, v2, e3->u.id) * wi3;
2668
2669 EIG_linear_solver_right_hand_side_add(context, 0, v2->u.id, j2[1][1] * beta[1]);
2670 EIG_linear_solver_right_hand_side_add(
2671 context, 0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
2672
2673 row1[1] = j2[1][1] * W[0][1];
2674 row2[1] = j2[1][1] * W[1][1];
2675 row3[1] = j2[1][1] * W[2][1];
2676
2677 row1[4] = j2[0][1] * W[0][0] + j2[2][1] * W[0][2];
2678 row2[4] = j2[0][1] * W[1][0] + j2[2][1] * W[1][2];
2679 row3[4] = j2[0][1] * W[2][0] + j2[2][1] * W[2][2];
2680 }
2681
2682 if (v3->flag & PVERT_INTERIOR) {
2683 vid[2] = v3->u.id;
2684 vid[5] = ninterior + v3->u.id;
2685
2686 sys->J2dt[e1->u.id][2] = j2[0][2] = p_abf_compute_sin_product(sys, v3, e1->u.id) * wi1;
2687 sys->J2dt[e2->u.id][2] = j2[1][2] = p_abf_compute_sin_product(sys, v3, e2->u.id) * wi2;
2688 sys->J2dt[e3->u.id][2] = j2[2][2] = 1.0f * wi3;
2689
2690 EIG_linear_solver_right_hand_side_add(context, 0, v3->u.id, j2[2][2] * beta[2]);
2691 EIG_linear_solver_right_hand_side_add(
2692 context, 0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
2693
2694 row1[2] = j2[2][2] * W[0][2];
2695 row2[2] = j2[2][2] * W[1][2];
2696 row3[2] = j2[2][2] * W[2][2];
2697
2698 row1[5] = j2[0][2] * W[0][0] + j2[1][2] * W[0][1];
2699 row2[5] = j2[0][2] * W[1][0] + j2[1][2] * W[1][1];
2700 row3[5] = j2[0][2] * W[2][0] + j2[1][2] * W[2][1];
2701 }
2702
2703 for (i = 0; i < 3; i++) {
2704 int r = vid[i];
2705
2706 if (r == -1) {
2707 continue;
2708 }
2709
2710 for (j = 0; j < 6; j++) {
2711 int c = vid[j];
2712
2713 if (c == -1) {
2714 continue;
2715 }
2716
2717 if (i == 0) {
2718 EIG_linear_solver_matrix_add(context, r, c, j2[0][i] * row1[j]);
2719 }
2720 else {
2721 EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[0][i] * row1[j]);
2722 }
2723
2724 if (i == 1) {
2725 EIG_linear_solver_matrix_add(context, r, c, j2[1][i] * row2[j]);
2726 }
2727 else {
2728 EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[1][i] * row2[j]);
2729 }
2730
2731 if (i == 2) {
2732 EIG_linear_solver_matrix_add(context, r, c, j2[2][i] * row3[j]);
2733 }
2734 else {
2735 EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[2][i] * row3[j]);
2736 }
2737 }
2738 }
2739 }
2740
2741 success = EIG_linear_solver_solve(context);
2742
2743 if (success) {
2744 for (f = chart->faces; f; f = f->nextlink) {
2745 float dlambda1, pre[3], dalpha;
2746 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
2747 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
2748
2749 pre[0] = pre[1] = pre[2] = 0.0;
2750
2751 if (v1->flag & PVERT_INTERIOR) {
2752 float x = EIG_linear_solver_variable_get(context, 0, v1->u.id);
2753 float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v1->u.id);
2754 pre[0] += sys->J2dt[e1->u.id][0] * x;
2755 pre[1] += sys->J2dt[e2->u.id][0] * x2;
2756 pre[2] += sys->J2dt[e3->u.id][0] * x2;
2757 }
2758
2759 if (v2->flag & PVERT_INTERIOR) {
2760 float x = EIG_linear_solver_variable_get(context, 0, v2->u.id);
2761 float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v2->u.id);
2762 pre[0] += sys->J2dt[e1->u.id][1] * x2;
2763 pre[1] += sys->J2dt[e2->u.id][1] * x;
2764 pre[2] += sys->J2dt[e3->u.id][1] * x2;
2765 }
2766
2767 if (v3->flag & PVERT_INTERIOR) {
2768 float x = EIG_linear_solver_variable_get(context, 0, v3->u.id);
2769 float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v3->u.id);
2770 pre[0] += sys->J2dt[e1->u.id][2] * x2;
2771 pre[1] += sys->J2dt[e2->u.id][2] * x2;
2772 pre[2] += sys->J2dt[e3->u.id][2] * x;
2773 }
2774
2775 dlambda1 = pre[0] + pre[1] + pre[2];
2776 dlambda1 = sys->dstar[f->u.id] * (sys->bstar[f->u.id] - dlambda1);
2777
2778 sys->lambdaTriangle[f->u.id] += dlambda1;
2779
2780 dalpha = (sys->bAlpha[e1->u.id] - dlambda1);
2781 sys->alpha[e1->u.id] += dalpha / sys->weight[e1->u.id] - pre[0];
2782
2783 dalpha = (sys->bAlpha[e2->u.id] - dlambda1);
2784 sys->alpha[e2->u.id] += dalpha / sys->weight[e2->u.id] - pre[1];
2785
2786 dalpha = (sys->bAlpha[e3->u.id] - dlambda1);
2787 sys->alpha[e3->u.id] += dalpha / sys->weight[e3->u.id] - pre[2];
2788
2789 /* clamp */
2790 e = f->edge;
2791 do {
2792 if (sys->alpha[e->u.id] > (float)M_PI) {
2793 sys->alpha[e->u.id] = (float)M_PI;
2794 }
2795 else if (sys->alpha[e->u.id] < 0.0f) {
2796 sys->alpha[e->u.id] = 0.0f;
2797 }
2798 } while (e != f->edge);
2799 }
2800
2801 for (i = 0; i < ninterior; i++) {
2802 sys->lambdaPlanar[i] += (float)EIG_linear_solver_variable_get(context, 0, i);
2803 sys->lambdaLength[i] += (float)EIG_linear_solver_variable_get(context, 0, ninterior + i);
2804 }
2805 }
2806
2807 EIG_linear_solver_delete(context);
2808
2809 return success;
2810 }
2811
p_chart_abf_solve(PChart * chart)2812 static PBool p_chart_abf_solve(PChart *chart)
2813 {
2814 PVert *v;
2815 PFace *f;
2816 PEdge *e, *e1, *e2, *e3;
2817 PAbfSystem sys;
2818 int i;
2819 float /* lastnorm, */ /* UNUSED */ limit = (chart->nfaces > 100) ? 1.0f : 0.001f;
2820
2821 /* setup id's */
2822 sys.ninterior = sys.nfaces = sys.nangles = 0;
2823
2824 for (v = chart->verts; v; v = v->nextlink) {
2825 if (p_vert_interior(v)) {
2826 v->flag |= PVERT_INTERIOR;
2827 v->u.id = sys.ninterior++;
2828 }
2829 else {
2830 v->flag &= ~PVERT_INTERIOR;
2831 }
2832 }
2833
2834 for (f = chart->faces; f; f = f->nextlink) {
2835 e1 = f->edge;
2836 e2 = e1->next;
2837 e3 = e2->next;
2838 f->u.id = sys.nfaces++;
2839
2840 /* angle id's are conveniently stored in half edges */
2841 e1->u.id = sys.nangles++;
2842 e2->u.id = sys.nangles++;
2843 e3->u.id = sys.nangles++;
2844 }
2845
2846 p_abf_setup_system(&sys);
2847
2848 /* compute initial angles */
2849 for (f = chart->faces; f; f = f->nextlink) {
2850 float a1, a2, a3;
2851
2852 e1 = f->edge;
2853 e2 = e1->next;
2854 e3 = e2->next;
2855 p_face_angles(f, &a1, &a2, &a3);
2856
2857 if (a1 < sys.minangle) {
2858 a1 = sys.minangle;
2859 }
2860 else if (a1 > sys.maxangle) {
2861 a1 = sys.maxangle;
2862 }
2863 if (a2 < sys.minangle) {
2864 a2 = sys.minangle;
2865 }
2866 else if (a2 > sys.maxangle) {
2867 a2 = sys.maxangle;
2868 }
2869 if (a3 < sys.minangle) {
2870 a3 = sys.minangle;
2871 }
2872 else if (a3 > sys.maxangle) {
2873 a3 = sys.maxangle;
2874 }
2875
2876 sys.alpha[e1->u.id] = sys.beta[e1->u.id] = a1;
2877 sys.alpha[e2->u.id] = sys.beta[e2->u.id] = a2;
2878 sys.alpha[e3->u.id] = sys.beta[e3->u.id] = a3;
2879
2880 sys.weight[e1->u.id] = 2.0f / (a1 * a1);
2881 sys.weight[e2->u.id] = 2.0f / (a2 * a2);
2882 sys.weight[e3->u.id] = 2.0f / (a3 * a3);
2883 }
2884
2885 for (v = chart->verts; v; v = v->nextlink) {
2886 if (v->flag & PVERT_INTERIOR) {
2887 float anglesum = 0.0, scale;
2888
2889 e = v->edge;
2890 do {
2891 anglesum += sys.beta[e->u.id];
2892 e = e->next->next->pair;
2893 } while (e && (e != v->edge));
2894
2895 scale = (anglesum == 0.0f) ? 0.0f : 2.0f * (float)M_PI / anglesum;
2896
2897 e = v->edge;
2898 do {
2899 sys.beta[e->u.id] = sys.alpha[e->u.id] = sys.beta[e->u.id] * scale;
2900 e = e->next->next->pair;
2901 } while (e && (e != v->edge));
2902 }
2903 }
2904
2905 if (sys.ninterior > 0) {
2906 p_abf_compute_sines(&sys);
2907
2908 /* iteration */
2909 /* lastnorm = 1e10; */ /* UNUSED */
2910
2911 for (i = 0; i < ABF_MAX_ITER; i++) {
2912 float norm = p_abf_compute_gradient(&sys, chart);
2913
2914 /* lastnorm = norm; */ /* UNUSED */
2915
2916 if (norm < limit) {
2917 break;
2918 }
2919
2920 if (!p_abf_matrix_invert(&sys, chart)) {
2921 param_warning("ABF failed to invert matrix");
2922 p_abf_free_system(&sys);
2923 return P_FALSE;
2924 }
2925
2926 p_abf_compute_sines(&sys);
2927 }
2928
2929 if (i == ABF_MAX_ITER) {
2930 param_warning("ABF maximum iterations reached");
2931 p_abf_free_system(&sys);
2932 return P_FALSE;
2933 }
2934 }
2935
2936 chart->u.lscm.abf_alpha = MEM_dupallocN(sys.alpha);
2937 p_abf_free_system(&sys);
2938
2939 return P_TRUE;
2940 }
2941
2942 /* Least Squares Conformal Maps */
2943
p_chart_pin_positions(PChart * chart,PVert ** pin1,PVert ** pin2)2944 static void p_chart_pin_positions(PChart *chart, PVert **pin1, PVert **pin2)
2945 {
2946 if (!*pin1 || !*pin2 || *pin1 == *pin2) {
2947 /* degenerate case */
2948 PFace *f = chart->faces;
2949 *pin1 = f->edge->vert;
2950 *pin2 = f->edge->next->vert;
2951
2952 (*pin1)->uv[0] = 0.0f;
2953 (*pin1)->uv[1] = 0.5f;
2954 (*pin2)->uv[0] = 1.0f;
2955 (*pin2)->uv[1] = 0.5f;
2956 }
2957 else {
2958 int diru, dirv, dirx, diry;
2959 float sub[3];
2960
2961 sub_v3_v3v3(sub, (*pin1)->co, (*pin2)->co);
2962 sub[0] = fabsf(sub[0]);
2963 sub[1] = fabsf(sub[1]);
2964 sub[2] = fabsf(sub[2]);
2965
2966 if ((sub[0] > sub[1]) && (sub[0] > sub[2])) {
2967 dirx = 0;
2968 diry = (sub[1] > sub[2]) ? 1 : 2;
2969 }
2970 else if ((sub[1] > sub[0]) && (sub[1] > sub[2])) {
2971 dirx = 1;
2972 diry = (sub[0] > sub[2]) ? 0 : 2;
2973 }
2974 else {
2975 dirx = 2;
2976 diry = (sub[0] > sub[1]) ? 0 : 1;
2977 }
2978
2979 if (dirx == 2) {
2980 diru = 1;
2981 dirv = 0;
2982 }
2983 else {
2984 diru = 0;
2985 dirv = 1;
2986 }
2987
2988 (*pin1)->uv[diru] = (*pin1)->co[dirx];
2989 (*pin1)->uv[dirv] = (*pin1)->co[diry];
2990 (*pin2)->uv[diru] = (*pin2)->co[dirx];
2991 (*pin2)->uv[dirv] = (*pin2)->co[diry];
2992 }
2993 }
2994
p_chart_symmetry_pins(PChart * chart,PEdge * outer,PVert ** pin1,PVert ** pin2)2995 static PBool p_chart_symmetry_pins(PChart *chart, PEdge *outer, PVert **pin1, PVert **pin2)
2996 {
2997 PEdge *be, *lastbe = NULL, *maxe1 = NULL, *maxe2 = NULL, *be1, *be2;
2998 PEdge *cure = NULL, *firste1 = NULL, *firste2 = NULL, *nextbe;
2999 float maxlen = 0.0f, curlen = 0.0f, totlen = 0.0f, firstlen = 0.0f;
3000 float len1, len2;
3001
3002 /* find longest series of verts split in the chart itself, these are
3003 * marked during construction */
3004 be = outer;
3005 lastbe = p_boundary_edge_prev(be);
3006 do {
3007 float len = p_edge_length(be);
3008 totlen += len;
3009
3010 nextbe = p_boundary_edge_next(be);
3011
3012 if ((be->vert->flag & PVERT_SPLIT) ||
3013 (lastbe->vert->flag & nextbe->vert->flag & PVERT_SPLIT)) {
3014 if (!cure) {
3015 if (be == outer) {
3016 firste1 = be;
3017 }
3018 cure = be;
3019 }
3020 else {
3021 curlen += p_edge_length(lastbe);
3022 }
3023 }
3024 else if (cure) {
3025 if (curlen > maxlen) {
3026 maxlen = curlen;
3027 maxe1 = cure;
3028 maxe2 = lastbe;
3029 }
3030
3031 if (firste1 == cure) {
3032 firstlen = curlen;
3033 firste2 = lastbe;
3034 }
3035
3036 curlen = 0.0f;
3037 cure = NULL;
3038 }
3039
3040 lastbe = be;
3041 be = nextbe;
3042 } while (be != outer);
3043
3044 /* make sure we also count a series of splits over the starting point */
3045 if (cure && (cure != outer)) {
3046 firstlen += curlen + p_edge_length(be);
3047
3048 if (firstlen > maxlen) {
3049 maxlen = firstlen;
3050 maxe1 = cure;
3051 maxe2 = firste2;
3052 }
3053 }
3054
3055 if (!maxe1 || !maxe2 || (maxlen < 0.5f * totlen)) {
3056 return P_FALSE;
3057 }
3058
3059 /* find pin1 in the split vertices */
3060 be1 = maxe1;
3061 be2 = maxe2;
3062 len1 = 0.0f;
3063 len2 = 0.0f;
3064
3065 do {
3066 if (len1 < len2) {
3067 len1 += p_edge_length(be1);
3068 be1 = p_boundary_edge_next(be1);
3069 }
3070 else {
3071 be2 = p_boundary_edge_prev(be2);
3072 len2 += p_edge_length(be2);
3073 }
3074 } while (be1 != be2);
3075
3076 *pin1 = be1->vert;
3077
3078 /* find pin2 outside the split vertices */
3079 be1 = maxe1;
3080 be2 = maxe2;
3081 len1 = 0.0f;
3082 len2 = 0.0f;
3083
3084 do {
3085 if (len1 < len2) {
3086 be1 = p_boundary_edge_prev(be1);
3087 len1 += p_edge_length(be1);
3088 }
3089 else {
3090 len2 += p_edge_length(be2);
3091 be2 = p_boundary_edge_next(be2);
3092 }
3093 } while (be1 != be2);
3094
3095 *pin2 = be1->vert;
3096
3097 p_chart_pin_positions(chart, pin1, pin2);
3098
3099 return !equals_v3v3((*pin1)->co, (*pin2)->co);
3100 }
3101
p_chart_extrema_verts(PChart * chart,PVert ** pin1,PVert ** pin2)3102 static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **pin2)
3103 {
3104 float minv[3], maxv[3], dirlen;
3105 PVert *v, *minvert[3], *maxvert[3];
3106 int i, dir;
3107
3108 /* find minimum and maximum verts over x/y/z axes */
3109 minv[0] = minv[1] = minv[2] = 1e20;
3110 maxv[0] = maxv[1] = maxv[2] = -1e20;
3111
3112 minvert[0] = minvert[1] = minvert[2] = NULL;
3113 maxvert[0] = maxvert[1] = maxvert[2] = NULL;
3114
3115 for (v = chart->verts; v; v = v->nextlink) {
3116 for (i = 0; i < 3; i++) {
3117 if (v->co[i] < minv[i]) {
3118 minv[i] = v->co[i];
3119 minvert[i] = v;
3120 }
3121 if (v->co[i] > maxv[i]) {
3122 maxv[i] = v->co[i];
3123 maxvert[i] = v;
3124 }
3125 }
3126 }
3127
3128 /* find axes with longest distance */
3129 dir = 0;
3130 dirlen = -1.0;
3131
3132 for (i = 0; i < 3; i++) {
3133 if (maxv[i] - minv[i] > dirlen) {
3134 dir = i;
3135 dirlen = maxv[i] - minv[i];
3136 }
3137 }
3138
3139 *pin1 = minvert[dir];
3140 *pin2 = maxvert[dir];
3141
3142 p_chart_pin_positions(chart, pin1, pin2);
3143 }
3144
p_chart_lscm_load_solution(PChart * chart)3145 static void p_chart_lscm_load_solution(PChart *chart)
3146 {
3147 LinearSolver *context = chart->u.lscm.context;
3148 PVert *v;
3149
3150 for (v = chart->verts; v; v = v->nextlink) {
3151 v->uv[0] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id);
3152 v->uv[1] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id + 1);
3153 }
3154 }
3155
p_chart_lscm_begin(PChart * chart,PBool live,PBool abf)3156 static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
3157 {
3158 PVert *v, *pin1, *pin2;
3159 PBool select = P_FALSE, deselect = P_FALSE;
3160 int npins = 0, id = 0;
3161
3162 /* give vertices matrix indices and count pins */
3163 for (v = chart->verts; v; v = v->nextlink) {
3164 if (v->flag & PVERT_PIN) {
3165 npins++;
3166 if (v->flag & PVERT_SELECT) {
3167 select = P_TRUE;
3168 }
3169 }
3170
3171 if (!(v->flag & PVERT_SELECT)) {
3172 deselect = P_TRUE;
3173 }
3174 }
3175
3176 if ((live && (!select || !deselect))) {
3177 chart->u.lscm.context = NULL;
3178 }
3179 else {
3180 #if 0
3181 p_chart_simplify_compute(chart);
3182 p_chart_topological_sanity_check(chart);
3183 #endif
3184
3185 if (npins == 1) {
3186 chart->u.lscm.single_pin_area = p_chart_uv_area(chart);
3187 for (v = chart->verts; v; v = v->nextlink) {
3188 if (v->flag & PVERT_PIN) {
3189 chart->u.lscm.single_pin = v;
3190 break;
3191 }
3192 }
3193 }
3194
3195 if (abf) {
3196 if (!p_chart_abf_solve(chart)) {
3197 param_warning("ABF solving failed: falling back to LSCM.\n");
3198 }
3199 }
3200
3201 if (npins <= 1) {
3202 /* No pins, let's find some ourself. */
3203 PEdge *outer;
3204
3205 p_chart_boundaries(chart, NULL, &outer);
3206
3207 /* Outer can be NULL with non-finite coords. */
3208 if (!(outer && p_chart_symmetry_pins(chart, outer, &pin1, &pin2))) {
3209 p_chart_extrema_verts(chart, &pin1, &pin2);
3210 }
3211
3212 chart->u.lscm.pin1 = pin1;
3213 chart->u.lscm.pin2 = pin2;
3214 }
3215
3216 for (v = chart->verts; v; v = v->nextlink) {
3217 v->u.id = id++;
3218 }
3219
3220 chart->u.lscm.context = EIG_linear_least_squares_solver_new(
3221 2 * chart->nfaces, 2 * chart->nverts, 1);
3222 }
3223 }
3224
p_chart_lscm_solve(PHandle * handle,PChart * chart)3225 static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
3226 {
3227 LinearSolver *context = chart->u.lscm.context;
3228 PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
3229 PFace *f;
3230 const float *alpha = chart->u.lscm.abf_alpha;
3231 float area_pinned_up, area_pinned_down;
3232 bool flip_faces;
3233 int row;
3234
3235 #if 0
3236 /* TODO: make loading pins work for simplify/complexify. */
3237 #endif
3238
3239 for (v = chart->verts; v; v = v->nextlink) {
3240 if (v->flag & PVERT_PIN) {
3241 p_vert_load_pin_select_uvs(handle, v); /* reload for live */
3242 }
3243 }
3244
3245 if (chart->u.lscm.single_pin) {
3246 /* If only one pin, save area and pin for transform later. */
3247 copy_v2_v2(chart->u.lscm.single_pin_uv, chart->u.lscm.single_pin->uv);
3248 }
3249
3250 if (chart->u.lscm.pin1) {
3251 EIG_linear_solver_variable_lock(context, 2 * pin1->u.id);
3252 EIG_linear_solver_variable_lock(context, 2 * pin1->u.id + 1);
3253 EIG_linear_solver_variable_lock(context, 2 * pin2->u.id);
3254 EIG_linear_solver_variable_lock(context, 2 * pin2->u.id + 1);
3255
3256 EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id, pin1->uv[0]);
3257 EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id + 1, pin1->uv[1]);
3258 EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id, pin2->uv[0]);
3259 EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id + 1, pin2->uv[1]);
3260 }
3261 else {
3262 /* set and lock the pins */
3263 for (v = chart->verts; v; v = v->nextlink) {
3264 if (v->flag & PVERT_PIN) {
3265 EIG_linear_solver_variable_lock(context, 2 * v->u.id);
3266 EIG_linear_solver_variable_lock(context, 2 * v->u.id + 1);
3267
3268 EIG_linear_solver_variable_set(context, 0, 2 * v->u.id, v->uv[0]);
3269 EIG_linear_solver_variable_set(context, 0, 2 * v->u.id + 1, v->uv[1]);
3270 }
3271 }
3272 }
3273
3274 /* detect up direction based on pinned vertices */
3275 area_pinned_up = 0.0f;
3276 area_pinned_down = 0.0f;
3277
3278 for (f = chart->faces; f; f = f->nextlink) {
3279 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
3280 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
3281
3282 if ((v1->flag & PVERT_PIN) && (v2->flag & PVERT_PIN) && (v3->flag & PVERT_PIN)) {
3283 float area = p_face_uv_area_signed(f);
3284
3285 if (area > 0.0f) {
3286 area_pinned_up += area;
3287 }
3288 else {
3289 area_pinned_down -= area;
3290 }
3291 }
3292 }
3293
3294 flip_faces = (area_pinned_down > area_pinned_up);
3295
3296 /* construct matrix */
3297
3298 row = 0;
3299 for (f = chart->faces; f; f = f->nextlink) {
3300 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
3301 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
3302 float a1, a2, a3, ratio, cosine, sine;
3303 float sina1, sina2, sina3, sinmax;
3304
3305 if (alpha) {
3306 /* use abf angles if passed on */
3307 a1 = *(alpha++);
3308 a2 = *(alpha++);
3309 a3 = *(alpha++);
3310 }
3311 else {
3312 p_face_angles(f, &a1, &a2, &a3);
3313 }
3314
3315 if (flip_faces) {
3316 SWAP(float, a2, a3);
3317 SWAP(PEdge *, e2, e3);
3318 SWAP(PVert *, v2, v3);
3319 }
3320
3321 sina1 = sinf(a1);
3322 sina2 = sinf(a2);
3323 sina3 = sinf(a3);
3324
3325 sinmax = max_fff(sina1, sina2, sina3);
3326
3327 /* shift vertices to find most stable order */
3328 if (sina3 != sinmax) {
3329 SHIFT3(PVert *, v1, v2, v3);
3330 SHIFT3(float, a1, a2, a3);
3331 SHIFT3(float, sina1, sina2, sina3);
3332
3333 if (sina2 == sinmax) {
3334 SHIFT3(PVert *, v1, v2, v3);
3335 SHIFT3(float, a1, a2, a3);
3336 SHIFT3(float, sina1, sina2, sina3);
3337 }
3338 }
3339
3340 /* angle based lscm formulation */
3341 ratio = (sina3 == 0.0f) ? 1.0f : sina2 / sina3;
3342 cosine = cosf(a1) * ratio;
3343 sine = sina1 * ratio;
3344
3345 EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id, cosine - 1.0f);
3346 EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, -sine);
3347 EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id, -cosine);
3348 EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, sine);
3349 EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id, 1.0);
3350 row++;
3351
3352 EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id, sine);
3353 EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, cosine - 1.0f);
3354 EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id, -sine);
3355 EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, -cosine);
3356 EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id + 1, 1.0);
3357 row++;
3358 }
3359
3360 if (EIG_linear_solver_solve(context)) {
3361 p_chart_lscm_load_solution(chart);
3362 return P_TRUE;
3363 }
3364
3365 for (v = chart->verts; v; v = v->nextlink) {
3366 v->uv[0] = 0.0f;
3367 v->uv[1] = 0.0f;
3368 }
3369
3370 return P_FALSE;
3371 }
3372
p_chart_lscm_transform_single_pin(PChart * chart)3373 static void p_chart_lscm_transform_single_pin(PChart *chart)
3374 {
3375 PVert *pin = chart->u.lscm.single_pin;
3376
3377 /* If only one pin, keep UV area the same. */
3378 const float new_area = p_chart_uv_area(chart);
3379 if (new_area > 0.0f) {
3380 const float scale = chart->u.lscm.single_pin_area / new_area;
3381 if (scale > 0.0f) {
3382 p_chart_uv_scale(chart, sqrtf(scale));
3383 }
3384 }
3385
3386 /* Translate to keep the pinned vertex in place. */
3387 float offset[2];
3388 sub_v2_v2v2(offset, chart->u.lscm.single_pin_uv, pin->uv);
3389 p_chart_uv_translate(chart, offset);
3390 }
3391
p_chart_lscm_end(PChart * chart)3392 static void p_chart_lscm_end(PChart *chart)
3393 {
3394 if (chart->u.lscm.context) {
3395 EIG_linear_solver_delete(chart->u.lscm.context);
3396 }
3397
3398 if (chart->u.lscm.abf_alpha) {
3399 MEM_freeN(chart->u.lscm.abf_alpha);
3400 chart->u.lscm.abf_alpha = NULL;
3401 }
3402
3403 chart->u.lscm.context = NULL;
3404 chart->u.lscm.pin1 = NULL;
3405 chart->u.lscm.pin2 = NULL;
3406 chart->u.lscm.single_pin = NULL;
3407 chart->u.lscm.single_pin_area = 0.0f;
3408 }
3409
3410 /* Stretch */
3411
3412 #define P_STRETCH_ITER 20
3413
p_stretch_pin_boundary(PChart * chart)3414 static void p_stretch_pin_boundary(PChart *chart)
3415 {
3416 PVert *v;
3417
3418 for (v = chart->verts; v; v = v->nextlink) {
3419 if (v->edge->pair == NULL) {
3420 v->flag |= PVERT_PIN;
3421 }
3422 else {
3423 v->flag &= ~PVERT_PIN;
3424 }
3425 }
3426 }
3427
p_face_stretch(PFace * f)3428 static float p_face_stretch(PFace *f)
3429 {
3430 float T, w, tmp[3];
3431 float Ps[3], Pt[3];
3432 float a, c, area;
3433 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
3434 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
3435
3436 area = p_face_uv_area_signed(f);
3437
3438 if (area <= 0.0f) { /* flipped face -> infinite stretch */
3439 return 1e10f;
3440 }
3441
3442 w = 1.0f / (2.0f * area);
3443
3444 /* compute derivatives */
3445 copy_v3_v3(Ps, v1->co);
3446 mul_v3_fl(Ps, (v2->uv[1] - v3->uv[1]));
3447
3448 copy_v3_v3(tmp, v2->co);
3449 mul_v3_fl(tmp, (v3->uv[1] - v1->uv[1]));
3450 add_v3_v3(Ps, tmp);
3451
3452 copy_v3_v3(tmp, v3->co);
3453 mul_v3_fl(tmp, (v1->uv[1] - v2->uv[1]));
3454 add_v3_v3(Ps, tmp);
3455
3456 mul_v3_fl(Ps, w);
3457
3458 copy_v3_v3(Pt, v1->co);
3459 mul_v3_fl(Pt, (v3->uv[0] - v2->uv[0]));
3460
3461 copy_v3_v3(tmp, v2->co);
3462 mul_v3_fl(tmp, (v1->uv[0] - v3->uv[0]));
3463 add_v3_v3(Pt, tmp);
3464
3465 copy_v3_v3(tmp, v3->co);
3466 mul_v3_fl(tmp, (v2->uv[0] - v1->uv[0]));
3467 add_v3_v3(Pt, tmp);
3468
3469 mul_v3_fl(Pt, w);
3470
3471 /* Sander Tensor */
3472 a = dot_v3v3(Ps, Ps);
3473 c = dot_v3v3(Pt, Pt);
3474
3475 T = sqrtf(0.5f * (a + c));
3476 if (f->flag & PFACE_FILLED) {
3477 T *= 0.2f;
3478 }
3479
3480 return T;
3481 }
3482
p_stretch_compute_vertex(PVert * v)3483 static float p_stretch_compute_vertex(PVert *v)
3484 {
3485 PEdge *e = v->edge;
3486 float sum = 0.0f;
3487
3488 do {
3489 sum += p_face_stretch(e->face);
3490 e = p_wheel_edge_next(e);
3491 } while (e && e != (v->edge));
3492
3493 return sum;
3494 }
3495
p_chart_stretch_minimize(PChart * chart,RNG * rng)3496 static void p_chart_stretch_minimize(PChart *chart, RNG *rng)
3497 {
3498 PVert *v;
3499 PEdge *e;
3500 int j, nedges;
3501 float orig_stretch, low, stretch_low, high, stretch_high, mid, stretch;
3502 float orig_uv[2], dir[2], random_angle, trusted_radius;
3503
3504 for (v = chart->verts; v; v = v->nextlink) {
3505 if ((v->flag & PVERT_PIN) || !(v->flag & PVERT_SELECT)) {
3506 continue;
3507 }
3508
3509 orig_stretch = p_stretch_compute_vertex(v);
3510 orig_uv[0] = v->uv[0];
3511 orig_uv[1] = v->uv[1];
3512
3513 /* move vertex in a random direction */
3514 trusted_radius = 0.0f;
3515 nedges = 0;
3516 e = v->edge;
3517
3518 do {
3519 trusted_radius += p_edge_uv_length(e);
3520 nedges++;
3521
3522 e = p_wheel_edge_next(e);
3523 } while (e && e != (v->edge));
3524
3525 trusted_radius /= 2 * nedges;
3526
3527 random_angle = BLI_rng_get_float(rng) * 2.0f * (float)M_PI;
3528 dir[0] = trusted_radius * cosf(random_angle);
3529 dir[1] = trusted_radius * sinf(random_angle);
3530
3531 /* calculate old and new stretch */
3532 low = 0;
3533 stretch_low = orig_stretch;
3534
3535 add_v2_v2v2(v->uv, orig_uv, dir);
3536 high = 1;
3537 stretch = stretch_high = p_stretch_compute_vertex(v);
3538
3539 /* binary search for lowest stretch position */
3540 for (j = 0; j < P_STRETCH_ITER; j++) {
3541 mid = 0.5f * (low + high);
3542 v->uv[0] = orig_uv[0] + mid * dir[0];
3543 v->uv[1] = orig_uv[1] + mid * dir[1];
3544 stretch = p_stretch_compute_vertex(v);
3545
3546 if (stretch_low < stretch_high) {
3547 high = mid;
3548 stretch_high = stretch;
3549 }
3550 else {
3551 low = mid;
3552 stretch_low = stretch;
3553 }
3554 }
3555
3556 /* no luck, stretch has increased, reset to old values */
3557 if (stretch >= orig_stretch) {
3558 copy_v2_v2(v->uv, orig_uv);
3559 }
3560 }
3561 }
3562
3563 /* Minimum area enclosing rectangle for packing */
3564
p_compare_geometric_uv(const void * a,const void * b)3565 static int p_compare_geometric_uv(const void *a, const void *b)
3566 {
3567 const PVert *v1 = *(const PVert *const *)a;
3568 const PVert *v2 = *(const PVert *const *)b;
3569
3570 if (v1->uv[0] < v2->uv[0]) {
3571 return -1;
3572 }
3573 if (v1->uv[0] == v2->uv[0]) {
3574 if (v1->uv[1] < v2->uv[1]) {
3575 return -1;
3576 }
3577 if (v1->uv[1] == v2->uv[1]) {
3578 return 0;
3579 }
3580 return 1;
3581 }
3582 return 1;
3583 }
3584
p_chart_convex_hull(PChart * chart,PVert *** r_verts,int * r_nverts,int * r_right)3585 static PBool p_chart_convex_hull(PChart *chart, PVert ***r_verts, int *r_nverts, int *r_right)
3586 {
3587 /* Graham algorithm, taken from:
3588 * http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/117225 */
3589
3590 PEdge *be, *e;
3591 int npoints = 0, i, ulen, llen;
3592 PVert **U, **L, **points, **p;
3593
3594 p_chart_boundaries(chart, NULL, &be);
3595
3596 if (!be) {
3597 return P_FALSE;
3598 }
3599
3600 e = be;
3601 do {
3602 npoints++;
3603 e = p_boundary_edge_next(e);
3604 } while (e != be);
3605
3606 p = points = (PVert **)MEM_mallocN(sizeof(PVert *) * npoints * 2, "PCHullpoints");
3607 U = (PVert **)MEM_mallocN(sizeof(PVert *) * npoints, "PCHullU");
3608 L = (PVert **)MEM_mallocN(sizeof(PVert *) * npoints, "PCHullL");
3609
3610 e = be;
3611 do {
3612 *p = e->vert;
3613 p++;
3614 e = p_boundary_edge_next(e);
3615 } while (e != be);
3616
3617 qsort(points, npoints, sizeof(PVert *), p_compare_geometric_uv);
3618
3619 ulen = llen = 0;
3620 for (p = points, i = 0; i < npoints; i++, p++) {
3621 while ((ulen > 1) && (p_area_signed(U[ulen - 2]->uv, (*p)->uv, U[ulen - 1]->uv) <= 0)) {
3622 ulen--;
3623 }
3624 while ((llen > 1) && (p_area_signed(L[llen - 2]->uv, (*p)->uv, L[llen - 1]->uv) >= 0)) {
3625 llen--;
3626 }
3627
3628 U[ulen] = *p;
3629 ulen++;
3630 L[llen] = *p;
3631 llen++;
3632 }
3633
3634 npoints = 0;
3635 for (p = points, i = 0; i < ulen; i++, p++, npoints++) {
3636 *p = U[i];
3637 }
3638
3639 /* the first and last point in L are left out, since they are also in U */
3640 for (i = llen - 2; i > 0; i--, p++, npoints++) {
3641 *p = L[i];
3642 }
3643
3644 *r_verts = points;
3645 *r_nverts = npoints;
3646 *r_right = ulen - 1;
3647
3648 MEM_freeN(U);
3649 MEM_freeN(L);
3650
3651 return P_TRUE;
3652 }
3653
p_rectangle_area(float * p1,float * dir,float * p2,float * p3,float * p4)3654 static float p_rectangle_area(float *p1, float *dir, float *p2, float *p3, float *p4)
3655 {
3656 /* given 4 points on the rectangle edges and the direction of on edge,
3657 * compute the area of the rectangle */
3658
3659 float orthodir[2], corner1[2], corner2[2], corner3[2];
3660
3661 orthodir[0] = dir[1];
3662 orthodir[1] = -dir[0];
3663
3664 if (!p_intersect_line_2d_dir(p1, dir, p2, orthodir, corner1)) {
3665 return 1e10;
3666 }
3667
3668 if (!p_intersect_line_2d_dir(p1, dir, p4, orthodir, corner2)) {
3669 return 1e10;
3670 }
3671
3672 if (!p_intersect_line_2d_dir(p3, dir, p4, orthodir, corner3)) {
3673 return 1e10;
3674 }
3675
3676 return len_v2v2(corner1, corner2) * len_v2v2(corner2, corner3);
3677 }
3678
p_chart_minimum_area_angle(PChart * chart)3679 static float p_chart_minimum_area_angle(PChart *chart)
3680 {
3681 /* minimum area enclosing rectangle with rotating calipers, info:
3682 * http://cgm.cs.mcgill.ca/~orm/maer.html */
3683
3684 float rotated, minarea, minangle, area, len;
3685 float *angles, miny, maxy, v[2], a[4], mina;
3686 int npoints, right, i_min, i_max, i, idx[4], nextidx;
3687 PVert **points, *p1, *p2, *p3, *p4, *p1n;
3688
3689 /* compute convex hull */
3690 if (!p_chart_convex_hull(chart, &points, &npoints, &right)) {
3691 return 0.0;
3692 }
3693
3694 /* find left/top/right/bottom points, and compute angle for each point */
3695 angles = MEM_mallocN(sizeof(float) * npoints, "PMinAreaAngles");
3696
3697 i_min = i_max = 0;
3698 miny = 1e10;
3699 maxy = -1e10;
3700
3701 for (i = 0; i < npoints; i++) {
3702 p1 = (i == 0) ? points[npoints - 1] : points[i - 1];
3703 p2 = points[i];
3704 p3 = (i == npoints - 1) ? points[0] : points[i + 1];
3705
3706 angles[i] = (float)M_PI - p_vec2_angle(p1->uv, p2->uv, p3->uv);
3707
3708 if (points[i]->uv[1] < miny) {
3709 miny = points[i]->uv[1];
3710 i_min = i;
3711 }
3712 if (points[i]->uv[1] > maxy) {
3713 maxy = points[i]->uv[1];
3714 i_max = i;
3715 }
3716 }
3717
3718 /* left, top, right, bottom */
3719 idx[0] = 0;
3720 idx[1] = i_max;
3721 idx[2] = right;
3722 idx[3] = i_min;
3723
3724 v[0] = points[idx[0]]->uv[0];
3725 v[1] = points[idx[0]]->uv[1] + 1.0f;
3726 a[0] = p_vec2_angle(points[(idx[0] + 1) % npoints]->uv, points[idx[0]]->uv, v);
3727
3728 v[0] = points[idx[1]]->uv[0] + 1.0f;
3729 v[1] = points[idx[1]]->uv[1];
3730 a[1] = p_vec2_angle(points[(idx[1] + 1) % npoints]->uv, points[idx[1]]->uv, v);
3731
3732 v[0] = points[idx[2]]->uv[0];
3733 v[1] = points[idx[2]]->uv[1] - 1.0f;
3734 a[2] = p_vec2_angle(points[(idx[2] + 1) % npoints]->uv, points[idx[2]]->uv, v);
3735
3736 v[0] = points[idx[3]]->uv[0] - 1.0f;
3737 v[1] = points[idx[3]]->uv[1];
3738 a[3] = p_vec2_angle(points[(idx[3] + 1) % npoints]->uv, points[idx[3]]->uv, v);
3739
3740 /* 4 rotating calipers */
3741
3742 rotated = 0.0;
3743 minarea = 1e10;
3744 minangle = 0.0;
3745
3746 while (rotated <= (float)(M_PI / 2.0)) { /* INVESTIGATE: how far to rotate? */
3747 /* rotate with the smallest angle */
3748 i_min = 0;
3749 mina = 1e10;
3750
3751 for (i = 0; i < 4; i++) {
3752 if (a[i] < mina) {
3753 mina = a[i];
3754 i_min = i;
3755 }
3756 }
3757
3758 rotated += mina;
3759 nextidx = (idx[i_min] + 1) % npoints;
3760
3761 a[i_min] = angles[nextidx];
3762 a[(i_min + 1) % 4] = a[(i_min + 1) % 4] - mina;
3763 a[(i_min + 2) % 4] = a[(i_min + 2) % 4] - mina;
3764 a[(i_min + 3) % 4] = a[(i_min + 3) % 4] - mina;
3765
3766 /* compute area */
3767 p1 = points[idx[i_min]];
3768 p1n = points[nextidx];
3769 p2 = points[idx[(i_min + 1) % 4]];
3770 p3 = points[idx[(i_min + 2) % 4]];
3771 p4 = points[idx[(i_min + 3) % 4]];
3772
3773 len = len_v2v2(p1->uv, p1n->uv);
3774
3775 if (len > 0.0f) {
3776 len = 1.0f / len;
3777 v[0] = (p1n->uv[0] - p1->uv[0]) * len;
3778 v[1] = (p1n->uv[1] - p1->uv[1]) * len;
3779
3780 area = p_rectangle_area(p1->uv, v, p2->uv, p3->uv, p4->uv);
3781
3782 /* remember smallest area */
3783 if (area < minarea) {
3784 minarea = area;
3785 minangle = rotated;
3786 }
3787 }
3788
3789 idx[i_min] = nextidx;
3790 }
3791
3792 /* try keeping rotation as small as possible */
3793 if (minangle > (float)(M_PI / 4)) {
3794 minangle -= (float)(M_PI / 2.0);
3795 }
3796
3797 MEM_freeN(angles);
3798 MEM_freeN(points);
3799
3800 return minangle;
3801 }
3802
p_chart_rotate_minimum_area(PChart * chart)3803 static void p_chart_rotate_minimum_area(PChart *chart)
3804 {
3805 float angle = p_chart_minimum_area_angle(chart);
3806 float sine = sinf(angle);
3807 float cosine = cosf(angle);
3808 PVert *v;
3809
3810 for (v = chart->verts; v; v = v->nextlink) {
3811 float oldu = v->uv[0], oldv = v->uv[1];
3812 v->uv[0] = cosine * oldu - sine * oldv;
3813 v->uv[1] = sine * oldu + cosine * oldv;
3814 }
3815 }
3816
p_chart_rotate_fit_aabb(PChart * chart)3817 static void p_chart_rotate_fit_aabb(PChart *chart)
3818 {
3819 float(*points)[2] = MEM_mallocN(sizeof(*points) * chart->nverts, __func__);
3820
3821 p_chart_uv_to_array(chart, points);
3822
3823 float angle = BLI_convexhull_aabb_fit_points_2d(points, chart->nverts);
3824
3825 MEM_freeN(points);
3826
3827 if (angle != 0.0f) {
3828 float mat[2][2];
3829 angle_to_mat2(mat, angle);
3830 p_chart_uv_transform(chart, mat);
3831 }
3832 }
3833
3834 /* Area Smoothing */
3835
3836 /* 2d bsp tree for inverse mapping - that's a bit silly */
3837
3838 typedef struct SmoothTriangle {
3839 float co1[2], co2[2], co3[2];
3840 float oco1[2], oco2[2], oco3[2];
3841 } SmoothTriangle;
3842
3843 typedef struct SmoothNode {
3844 struct SmoothNode *c1, *c2;
3845 SmoothTriangle **tri;
3846 float split;
3847 int axis, ntri;
3848 } SmoothNode;
3849
p_barycentric_2d(const float v1[2],const float v2[2],const float v3[2],const float p[2],float b[3])3850 static void p_barycentric_2d(
3851 const float v1[2], const float v2[2], const float v3[2], const float p[2], float b[3])
3852 {
3853 float a[2], c[2], h[2], div;
3854
3855 a[0] = v2[0] - v1[0];
3856 a[1] = v2[1] - v1[1];
3857 c[0] = v3[0] - v1[0];
3858 c[1] = v3[1] - v1[1];
3859
3860 div = a[0] * c[1] - a[1] * c[0];
3861
3862 if (div == 0.0f) {
3863 b[0] = 1.0f / 3.0f;
3864 b[1] = 1.0f / 3.0f;
3865 b[2] = 1.0f / 3.0f;
3866 }
3867 else {
3868 h[0] = p[0] - v1[0];
3869 h[1] = p[1] - v1[1];
3870
3871 div = 1.0f / div;
3872
3873 b[1] = (h[0] * c[1] - h[1] * c[0]) * div;
3874 b[2] = (a[0] * h[1] - a[1] * h[0]) * div;
3875 b[0] = 1.0f - b[1] - b[2];
3876 }
3877 }
3878
p_triangle_inside(SmoothTriangle * t,float co[2])3879 static PBool p_triangle_inside(SmoothTriangle *t, float co[2])
3880 {
3881 float b[3];
3882
3883 p_barycentric_2d(t->co1, t->co2, t->co3, co, b);
3884
3885 if ((b[0] >= 0.0f) && (b[1] >= 0.0f) && (b[2] >= 0.0f)) {
3886 co[0] = t->oco1[0] * b[0] + t->oco2[0] * b[1] + t->oco3[0] * b[2];
3887 co[1] = t->oco1[1] * b[0] + t->oco2[1] * b[1] + t->oco3[1] * b[2];
3888 return P_TRUE;
3889 }
3890
3891 return P_FALSE;
3892 }
3893
p_node_new(MemArena * arena,SmoothTriangle ** tri,int ntri,float * bmin,float * bmax,int depth)3894 static SmoothNode *p_node_new(
3895 MemArena *arena, SmoothTriangle **tri, int ntri, float *bmin, float *bmax, int depth)
3896 {
3897 SmoothNode *node = BLI_memarena_alloc(arena, sizeof(*node));
3898 int axis, i, t1size = 0, t2size = 0;
3899 float split, /* mi, */ /* UNUSED */ mx;
3900 SmoothTriangle **t1, **t2, *t;
3901
3902 node->tri = tri;
3903 node->ntri = ntri;
3904
3905 if (ntri <= 10 || depth >= 15) {
3906 return node;
3907 }
3908
3909 t1 = MEM_mallocN(sizeof(*t1) * ntri, "PNodeTri1");
3910 t2 = MEM_mallocN(sizeof(*t2) * ntri, "PNodeTri1");
3911
3912 axis = (bmax[0] - bmin[0] > bmax[1] - bmin[1]) ? 0 : 1;
3913 split = 0.5f * (bmin[axis] + bmax[axis]);
3914
3915 for (i = 0; i < ntri; i++) {
3916 t = tri[i];
3917
3918 if ((t->co1[axis] <= split) || (t->co2[axis] <= split) || (t->co3[axis] <= split)) {
3919 t1[t1size] = t;
3920 t1size++;
3921 }
3922 if ((t->co1[axis] >= split) || (t->co2[axis] >= split) || (t->co3[axis] >= split)) {
3923 t2[t2size] = t;
3924 t2size++;
3925 }
3926 }
3927
3928 if ((t1size == t2size) && (t1size == ntri)) {
3929 MEM_freeN(t1);
3930 MEM_freeN(t2);
3931 return node;
3932 }
3933
3934 node->tri = NULL;
3935 node->ntri = 0;
3936 MEM_freeN(tri);
3937
3938 node->axis = axis;
3939 node->split = split;
3940
3941 /* mi = bmin[axis]; */ /* UNUSED */
3942 mx = bmax[axis];
3943 bmax[axis] = split;
3944 node->c1 = p_node_new(arena, t1, t1size, bmin, bmax, depth + 1);
3945
3946 bmin[axis] = bmax[axis];
3947 bmax[axis] = mx;
3948 node->c2 = p_node_new(arena, t2, t2size, bmin, bmax, depth + 1);
3949
3950 return node;
3951 }
3952
p_node_delete(SmoothNode * node)3953 static void p_node_delete(SmoothNode *node)
3954 {
3955 if (node->c1) {
3956 p_node_delete(node->c1);
3957 }
3958 if (node->c2) {
3959 p_node_delete(node->c2);
3960 }
3961 if (node->tri) {
3962 MEM_freeN(node->tri);
3963 }
3964 }
3965
p_node_intersect(SmoothNode * node,float co[2])3966 static PBool p_node_intersect(SmoothNode *node, float co[2])
3967 {
3968 int i;
3969
3970 if (node->tri) {
3971 for (i = 0; i < node->ntri; i++) {
3972 if (p_triangle_inside(node->tri[i], co)) {
3973 return P_TRUE;
3974 }
3975 }
3976
3977 return P_FALSE;
3978 }
3979
3980 if (co[node->axis] < node->split) {
3981 return p_node_intersect(node->c1, co);
3982 }
3983 return p_node_intersect(node->c2, co);
3984 }
3985
3986 /* smoothing */
3987
p_compare_float(const void * a_,const void * b_)3988 static int p_compare_float(const void *a_, const void *b_)
3989 {
3990 const float a = *(const float *)a_;
3991 const float b = *(const float *)b_;
3992
3993 if (a < b) {
3994 return -1;
3995 }
3996 if (a == b) {
3997 return 0;
3998 }
3999 return 1;
4000 }
4001
p_smooth_median_edge_length(PChart * chart)4002 static float p_smooth_median_edge_length(PChart *chart)
4003 {
4004 PEdge *e;
4005 float *lengths = MEM_mallocN(sizeof(chart->edges) * chart->nedges, "PMedianLength");
4006 float median;
4007 int i;
4008
4009 /* ok, so I'm lazy */
4010 for (i = 0, e = chart->edges; e; e = e->nextlink, i++) {
4011 lengths[i] = p_edge_length(e);
4012 }
4013
4014 qsort(lengths, i, sizeof(float), p_compare_float);
4015
4016 median = lengths[i / 2];
4017 MEM_freeN(lengths);
4018
4019 return median;
4020 }
4021
p_smooth_distortion(PEdge * e,float avg2d,float avg3d)4022 static float p_smooth_distortion(PEdge *e, float avg2d, float avg3d)
4023 {
4024 float len2d = p_edge_uv_length(e) * avg3d;
4025 float len3d = p_edge_length(e) * avg2d;
4026
4027 return (len3d == 0.0f) ? 0.0f : len2d / len3d;
4028 }
4029
p_smooth(PChart * chart)4030 static void p_smooth(PChart *chart)
4031 {
4032 PEdge *e;
4033 PVert *v;
4034 PFace *f;
4035 int j, it2, maxiter2, it;
4036 int nedges = chart->nedges, nwheel, gridx, gridy;
4037 int edgesx, edgesy, nsize, esize, i, x, y, maxiter, totiter;
4038 float minv[2], maxv[2], median, invmedian, avglen2d, avglen3d;
4039 float center[2], dx, dy, *nodes, dlimit, d, *oldnodesx, *oldnodesy;
4040 float *nodesx, *nodesy, *hedges, *vedges, climit, moved, padding;
4041 SmoothTriangle *triangles, *t, *t2, **tri, **trip;
4042 SmoothNode *root;
4043 MemArena *arena;
4044
4045 if (nedges == 0) {
4046 return;
4047 }
4048
4049 p_chart_uv_bbox(chart, minv, maxv);
4050 median = p_smooth_median_edge_length(chart) * 0.10f;
4051
4052 if (median == 0.0f) {
4053 return;
4054 }
4055
4056 invmedian = 1.0f / median;
4057
4058 /* compute edge distortion */
4059 avglen2d = avglen3d = 0.0;
4060
4061 for (e = chart->edges; e; e = e->nextlink) {
4062 avglen2d += p_edge_uv_length(e);
4063 avglen3d += p_edge_length(e);
4064 }
4065
4066 avglen2d /= nedges;
4067 avglen3d /= nedges;
4068
4069 for (v = chart->verts; v; v = v->nextlink) {
4070 v->u.distortion = 0.0;
4071 nwheel = 0;
4072
4073 e = v->edge;
4074 do {
4075 v->u.distortion += p_smooth_distortion(e, avglen2d, avglen3d);
4076 nwheel++;
4077
4078 e = e->next->next->pair;
4079 } while (e && (e != v->edge));
4080
4081 v->u.distortion /= nwheel;
4082 }
4083
4084 /* need to do excessive grid size checking still */
4085 center[0] = 0.5f * (minv[0] + maxv[0]);
4086 center[1] = 0.5f * (minv[1] + maxv[1]);
4087
4088 dx = 0.5f * (maxv[0] - minv[0]);
4089 dy = 0.5f * (maxv[1] - minv[1]);
4090
4091 padding = 0.15f;
4092 dx += padding * dx + 2.0f * median;
4093 dy += padding * dy + 2.0f * median;
4094
4095 gridx = (int)(dx * invmedian);
4096 gridy = (int)(dy * invmedian);
4097
4098 minv[0] = center[0] - median * gridx;
4099 minv[1] = center[1] - median * gridy;
4100 maxv[0] = center[0] + median * gridx;
4101 maxv[1] = center[1] + median * gridy;
4102
4103 /* create grid */
4104 gridx = gridx * 2 + 1;
4105 gridy = gridy * 2 + 1;
4106
4107 if ((gridx <= 2) || (gridy <= 2)) {
4108 return;
4109 }
4110
4111 edgesx = gridx - 1;
4112 edgesy = gridy - 1;
4113 nsize = gridx * gridy;
4114 esize = edgesx * edgesy;
4115
4116 nodes = MEM_mallocN(sizeof(float) * nsize, "PSmoothNodes");
4117 nodesx = MEM_mallocN(sizeof(float) * nsize, "PSmoothNodesX");
4118 nodesy = MEM_mallocN(sizeof(float) * nsize, "PSmoothNodesY");
4119 oldnodesx = MEM_mallocN(sizeof(float) * nsize, "PSmoothOldNodesX");
4120 oldnodesy = MEM_mallocN(sizeof(float) * nsize, "PSmoothOldNodesY");
4121 hedges = MEM_mallocN(sizeof(float) * esize, "PSmoothHEdges");
4122 vedges = MEM_mallocN(sizeof(float) * esize, "PSmoothVEdges");
4123
4124 if (!nodes || !nodesx || !nodesy || !oldnodesx || !oldnodesy || !hedges || !vedges) {
4125 if (nodes) {
4126 MEM_freeN(nodes);
4127 }
4128 if (nodesx) {
4129 MEM_freeN(nodesx);
4130 }
4131 if (nodesy) {
4132 MEM_freeN(nodesy);
4133 }
4134 if (oldnodesx) {
4135 MEM_freeN(oldnodesx);
4136 }
4137 if (oldnodesy) {
4138 MEM_freeN(oldnodesy);
4139 }
4140 if (hedges) {
4141 MEM_freeN(hedges);
4142 }
4143 if (vedges) {
4144 MEM_freeN(vedges);
4145 }
4146
4147 // printf("Not enough memory for area smoothing grid");
4148 return;
4149 }
4150
4151 for (x = 0; x < gridx; x++) {
4152 for (y = 0; y < gridy; y++) {
4153 i = x + y * gridx;
4154
4155 nodesx[i] = minv[0] + median * x;
4156 nodesy[i] = minv[1] + median * y;
4157
4158 nodes[i] = 1.0f;
4159 }
4160 }
4161
4162 /* embed in grid */
4163 for (f = chart->faces; f; f = f->nextlink) {
4164 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
4165 float fmin[2], fmax[2];
4166 int bx1, by1, bx2, by2;
4167
4168 INIT_MINMAX2(fmin, fmax);
4169
4170 minmax_v2v2_v2(fmin, fmax, e1->vert->uv);
4171 minmax_v2v2_v2(fmin, fmax, e2->vert->uv);
4172 minmax_v2v2_v2(fmin, fmax, e3->vert->uv);
4173
4174 bx1 = (int)((fmin[0] - minv[0]) * invmedian);
4175 by1 = (int)((fmin[1] - minv[1]) * invmedian);
4176 bx2 = (int)((fmax[0] - minv[0]) * invmedian + 2);
4177 by2 = (int)((fmax[1] - minv[1]) * invmedian + 2);
4178
4179 for (x = bx1; x < bx2; x++) {
4180 for (y = by1; y < by2; y++) {
4181 float p[2], b[3];
4182
4183 i = x + y * gridx;
4184
4185 p[0] = nodesx[i];
4186 p[1] = nodesy[i];
4187
4188 p_barycentric_2d(e1->vert->uv, e2->vert->uv, e3->vert->uv, p, b);
4189
4190 if ((b[0] > 0.0f) && (b[1] > 0.0f) && (b[2] > 0.0f)) {
4191 nodes[i] = e1->vert->u.distortion * b[0];
4192 nodes[i] += e2->vert->u.distortion * b[1];
4193 nodes[i] += e3->vert->u.distortion * b[2];
4194 }
4195 }
4196 }
4197 }
4198
4199 /* smooth the grid */
4200 maxiter = 10;
4201 totiter = 0;
4202 climit = 0.00001f * nsize;
4203
4204 for (it = 0; it < maxiter; it++) {
4205 moved = 0.0f;
4206
4207 for (x = 0; x < edgesx; x++) {
4208 for (y = 0; y < edgesy; y++) {
4209 i = x + y * gridx;
4210 j = x + y * edgesx;
4211
4212 hedges[j] = (nodes[i] + nodes[i + 1]) * 0.5f;
4213 vedges[j] = (nodes[i] + nodes[i + gridx]) * 0.5f;
4214
4215 /* we do *inverse* mapping */
4216 hedges[j] = 1.0f / hedges[j];
4217 vedges[j] = 1.0f / vedges[j];
4218 }
4219 }
4220
4221 maxiter2 = 50;
4222 dlimit = 0.0001f;
4223
4224 for (it2 = 0; it2 < maxiter2; it2++) {
4225 d = 0.0f;
4226 totiter += 1;
4227
4228 memcpy(oldnodesx, nodesx, sizeof(float) * nsize);
4229 memcpy(oldnodesy, nodesy, sizeof(float) * nsize);
4230
4231 for (x = 1; x < gridx - 1; x++) {
4232 for (y = 1; y < gridy - 1; y++) {
4233 float p[2], oldp[2], sum1, sum2, diff[2], length;
4234
4235 i = x + gridx * y;
4236 j = x + edgesx * y;
4237
4238 oldp[0] = oldnodesx[i];
4239 oldp[1] = oldnodesy[i];
4240
4241 sum1 = hedges[j - 1] * oldnodesx[i - 1];
4242 sum1 += hedges[j] * oldnodesx[i + 1];
4243 sum1 += vedges[j - edgesx] * oldnodesx[i - gridx];
4244 sum1 += vedges[j] * oldnodesx[i + gridx];
4245
4246 sum2 = hedges[j - 1];
4247 sum2 += hedges[j];
4248 sum2 += vedges[j - edgesx];
4249 sum2 += vedges[j];
4250
4251 nodesx[i] = sum1 / sum2;
4252
4253 sum1 = hedges[j - 1] * oldnodesy[i - 1];
4254 sum1 += hedges[j] * oldnodesy[i + 1];
4255 sum1 += vedges[j - edgesx] * oldnodesy[i - gridx];
4256 sum1 += vedges[j] * oldnodesy[i + gridx];
4257
4258 nodesy[i] = sum1 / sum2;
4259
4260 p[0] = nodesx[i];
4261 p[1] = nodesy[i];
4262
4263 diff[0] = p[0] - oldp[0];
4264 diff[1] = p[1] - oldp[1];
4265
4266 length = len_v2(diff);
4267 d = max_ff(d, length);
4268 moved += length;
4269 }
4270 }
4271
4272 if (d < dlimit) {
4273 break;
4274 }
4275 }
4276
4277 if (moved < climit) {
4278 break;
4279 }
4280 }
4281
4282 MEM_freeN(oldnodesx);
4283 MEM_freeN(oldnodesy);
4284 MEM_freeN(hedges);
4285 MEM_freeN(vedges);
4286
4287 /* create bsp */
4288 t = triangles = MEM_mallocN(sizeof(SmoothTriangle) * esize * 2, "PSmoothTris");
4289 trip = tri = MEM_mallocN(sizeof(SmoothTriangle *) * esize * 2, "PSmoothTriP");
4290
4291 if (!triangles || !tri) {
4292 MEM_freeN(nodes);
4293 MEM_freeN(nodesx);
4294 MEM_freeN(nodesy);
4295
4296 if (triangles) {
4297 MEM_freeN(triangles);
4298 }
4299 if (tri) {
4300 MEM_freeN(tri);
4301 }
4302
4303 // printf("Not enough memory for area smoothing grid");
4304 return;
4305 }
4306
4307 for (x = 0; x < edgesx; x++) {
4308 for (y = 0; y < edgesy; y++) {
4309 i = x + y * gridx;
4310
4311 t->co1[0] = nodesx[i];
4312 t->co1[1] = nodesy[i];
4313
4314 t->co2[0] = nodesx[i + 1];
4315 t->co2[1] = nodesy[i + 1];
4316
4317 t->co3[0] = nodesx[i + gridx];
4318 t->co3[1] = nodesy[i + gridx];
4319
4320 t->oco1[0] = minv[0] + x * median;
4321 t->oco1[1] = minv[1] + y * median;
4322
4323 t->oco2[0] = minv[0] + (x + 1) * median;
4324 t->oco2[1] = minv[1] + y * median;
4325
4326 t->oco3[0] = minv[0] + x * median;
4327 t->oco3[1] = minv[1] + (y + 1) * median;
4328
4329 t2 = t + 1;
4330
4331 t2->co1[0] = nodesx[i + gridx + 1];
4332 t2->co1[1] = nodesy[i + gridx + 1];
4333
4334 t2->oco1[0] = minv[0] + (x + 1) * median;
4335 t2->oco1[1] = minv[1] + (y + 1) * median;
4336
4337 t2->co2[0] = t->co2[0];
4338 t2->co2[1] = t->co2[1];
4339 t2->oco2[0] = t->oco2[0];
4340 t2->oco2[1] = t->oco2[1];
4341
4342 t2->co3[0] = t->co3[0];
4343 t2->co3[1] = t->co3[1];
4344 t2->oco3[0] = t->oco3[0];
4345 t2->oco3[1] = t->oco3[1];
4346
4347 *trip = t;
4348 trip++;
4349 t++;
4350 *trip = t;
4351 trip++;
4352 t++;
4353 }
4354 }
4355
4356 MEM_freeN(nodes);
4357 MEM_freeN(nodesx);
4358 MEM_freeN(nodesy);
4359
4360 arena = BLI_memarena_new(MEM_SIZE_OPTIMAL(1 << 16), "param smooth arena");
4361 root = p_node_new(arena, tri, esize * 2, minv, maxv, 0);
4362
4363 for (v = chart->verts; v; v = v->nextlink) {
4364 if (!p_node_intersect(root, v->uv)) {
4365 param_warning("area smoothing error: couldn't find mapping triangle\n");
4366 }
4367 }
4368
4369 p_node_delete(root);
4370 BLI_memarena_free(arena);
4371
4372 MEM_freeN(triangles);
4373 }
4374
4375 /* Exported */
4376
param_construct_begin(void)4377 ParamHandle *param_construct_begin(void)
4378 {
4379 PHandle *handle = MEM_callocN(sizeof(*handle), "PHandle");
4380 handle->construction_chart = p_chart_new(handle);
4381 handle->state = PHANDLE_STATE_ALLOCATED;
4382 handle->arena = BLI_memarena_new(MEM_SIZE_OPTIMAL(1 << 16), "param construct arena");
4383 handle->polyfill_arena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "param polyfill arena");
4384 handle->polyfill_heap = BLI_heap_new_ex(BLI_POLYFILL_ALLOC_NGON_RESERVE);
4385 handle->aspx = 1.0f;
4386 handle->aspy = 1.0f;
4387 handle->do_aspect = false;
4388
4389 handle->hash_verts = phash_new((PHashLink **)&handle->construction_chart->verts, 1);
4390 handle->hash_edges = phash_new((PHashLink **)&handle->construction_chart->edges, 1);
4391 handle->hash_faces = phash_new((PHashLink **)&handle->construction_chart->faces, 1);
4392
4393 return (ParamHandle *)handle;
4394 }
4395
param_aspect_ratio(ParamHandle * handle,float aspx,float aspy)4396 void param_aspect_ratio(ParamHandle *handle, float aspx, float aspy)
4397 {
4398 PHandle *phandle = (PHandle *)handle;
4399
4400 phandle->aspx = aspx;
4401 phandle->aspy = aspy;
4402 phandle->do_aspect = true;
4403 }
4404
param_delete(ParamHandle * handle)4405 void param_delete(ParamHandle *handle)
4406 {
4407 PHandle *phandle = (PHandle *)handle;
4408 int i;
4409
4410 param_assert((phandle->state == PHANDLE_STATE_ALLOCATED) ||
4411 (phandle->state == PHANDLE_STATE_CONSTRUCTED));
4412
4413 for (i = 0; i < phandle->ncharts; i++) {
4414 p_chart_delete(phandle->charts[i]);
4415 }
4416
4417 if (phandle->charts) {
4418 MEM_freeN(phandle->charts);
4419 }
4420
4421 if (phandle->construction_chart) {
4422 p_chart_delete(phandle->construction_chart);
4423
4424 phash_delete(phandle->hash_verts);
4425 phash_delete(phandle->hash_edges);
4426 phash_delete(phandle->hash_faces);
4427 }
4428
4429 BLI_memarena_free(phandle->arena);
4430 BLI_memarena_free(phandle->polyfill_arena);
4431 BLI_heap_free(phandle->polyfill_heap, NULL);
4432 MEM_freeN(phandle);
4433 }
4434
p_add_ngon(ParamHandle * handle,ParamKey key,int nverts,ParamKey * vkeys,float ** co,float ** uv,ParamBool * pin,ParamBool * select)4435 static void p_add_ngon(ParamHandle *handle,
4436 ParamKey key,
4437 int nverts,
4438 ParamKey *vkeys,
4439 float **co,
4440 float **uv,
4441 ParamBool *pin,
4442 ParamBool *select)
4443 {
4444 /* Allocate memory for polyfill. */
4445 PHandle *phandle = (PHandle *)handle;
4446 MemArena *arena = phandle->polyfill_arena;
4447 Heap *heap = phandle->polyfill_heap;
4448 uint nfilltri = nverts - 2;
4449 uint(*tris)[3] = BLI_memarena_alloc(arena, sizeof(*tris) * (size_t)nfilltri);
4450 float(*projverts)[2] = BLI_memarena_alloc(arena, sizeof(*projverts) * (size_t)nverts);
4451
4452 /* Calc normal, flipped: to get a positive 2d cross product. */
4453 float normal[3];
4454 zero_v3(normal);
4455
4456 const float *co_curr, *co_prev = co[nverts - 1];
4457 for (int j = 0; j < nverts; j++) {
4458 co_curr = co[j];
4459 add_newell_cross_v3_v3v3(normal, co_prev, co_curr);
4460 co_prev = co_curr;
4461 }
4462 if (UNLIKELY(normalize_v3(normal) == 0.0f)) {
4463 normal[2] = 1.0f;
4464 }
4465
4466 /* Project verts to 2d. */
4467 float axis_mat[3][3];
4468 axis_dominant_v3_to_m3_negate(axis_mat, normal);
4469 for (int j = 0; j < nverts; j++) {
4470 mul_v2_m3v3(projverts[j], axis_mat, co[j]);
4471 }
4472
4473 BLI_polyfill_calc_arena(projverts, nverts, 1, tris, arena);
4474
4475 /* Beautify helps avoid thin triangles that give numerical problems. */
4476 BLI_polyfill_beautify(projverts, nverts, tris, arena, heap);
4477
4478 /* Add triangles. */
4479 for (int j = 0; j < nfilltri; j++) {
4480 uint *tri = tris[j];
4481 uint v0 = tri[0];
4482 uint v1 = tri[1];
4483 uint v2 = tri[2];
4484
4485 ParamKey tri_vkeys[3] = {vkeys[v0], vkeys[v1], vkeys[v2]};
4486 float *tri_co[3] = {co[v0], co[v1], co[v2]};
4487 float *tri_uv[3] = {uv[v0], uv[v1], uv[v2]};
4488 ParamBool tri_pin[3] = {pin[v0], pin[v1], pin[v2]};
4489 ParamBool tri_select[3] = {select[v0], select[v1], select[v2]};
4490
4491 param_face_add(handle, key, 3, tri_vkeys, tri_co, tri_uv, tri_pin, tri_select);
4492 }
4493
4494 BLI_memarena_clear(arena);
4495 }
4496
param_face_add(ParamHandle * handle,ParamKey key,int nverts,ParamKey * vkeys,float * co[4],float * uv[4],ParamBool * pin,ParamBool * select)4497 void param_face_add(ParamHandle *handle,
4498 ParamKey key,
4499 int nverts,
4500 ParamKey *vkeys,
4501 float *co[4],
4502 float *uv[4],
4503 ParamBool *pin,
4504 ParamBool *select)
4505 {
4506 PHandle *phandle = (PHandle *)handle;
4507
4508 param_assert(phash_lookup(phandle->hash_faces, key) == NULL);
4509 param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
4510 param_assert((nverts == 3) || (nverts == 4));
4511
4512 if (nverts > 4) {
4513 /* ngon */
4514 p_add_ngon(handle, key, nverts, vkeys, co, uv, pin, select);
4515 }
4516 else if (nverts == 4) {
4517 /* quad */
4518 if (p_quad_split_direction(phandle, co, vkeys)) {
4519 p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 2, pin, select);
4520 p_face_add_construct(phandle, key, vkeys, co, uv, 0, 2, 3, pin, select);
4521 }
4522 else {
4523 p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 3, pin, select);
4524 p_face_add_construct(phandle, key, vkeys, co, uv, 1, 2, 3, pin, select);
4525 }
4526 }
4527 else if (!p_face_exists(phandle, vkeys, 0, 1, 2)) {
4528 /* triangle */
4529 p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 2, pin, select);
4530 }
4531 }
4532
param_edge_set_seam(ParamHandle * handle,ParamKey * vkeys)4533 void param_edge_set_seam(ParamHandle *handle, ParamKey *vkeys)
4534 {
4535 PHandle *phandle = (PHandle *)handle;
4536 PEdge *e;
4537
4538 param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
4539
4540 e = p_edge_lookup(phandle, vkeys);
4541 if (e) {
4542 e->flag |= PEDGE_SEAM;
4543 }
4544 }
4545
param_construct_end(ParamHandle * handle,ParamBool fill,ParamBool impl)4546 void param_construct_end(ParamHandle *handle, ParamBool fill, ParamBool impl)
4547 {
4548 PHandle *phandle = (PHandle *)handle;
4549 PChart *chart = phandle->construction_chart;
4550 int i, j, nboundaries = 0;
4551 PEdge *outer;
4552
4553 param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
4554
4555 phandle->ncharts = p_connect_pairs(phandle, (PBool)impl);
4556 phandle->charts = p_split_charts(phandle, chart, phandle->ncharts);
4557
4558 p_chart_delete(phandle->construction_chart);
4559 phandle->construction_chart = NULL;
4560
4561 phash_delete(phandle->hash_verts);
4562 phash_delete(phandle->hash_edges);
4563 phash_delete(phandle->hash_faces);
4564 phandle->hash_verts = phandle->hash_edges = phandle->hash_faces = NULL;
4565
4566 for (i = j = 0; i < phandle->ncharts; i++) {
4567 PVert *v;
4568 chart = phandle->charts[i];
4569
4570 p_chart_boundaries(chart, &nboundaries, &outer);
4571
4572 if (!impl && nboundaries == 0) {
4573 p_chart_delete(chart);
4574 continue;
4575 }
4576
4577 phandle->charts[j] = chart;
4578 j++;
4579
4580 if (fill && (nboundaries > 1)) {
4581 p_chart_fill_boundaries(chart, outer);
4582 }
4583
4584 for (v = chart->verts; v; v = v->nextlink) {
4585 p_vert_load_pin_select_uvs(handle, v);
4586 }
4587 }
4588
4589 phandle->ncharts = j;
4590
4591 phandle->state = PHANDLE_STATE_CONSTRUCTED;
4592 }
4593
param_lscm_begin(ParamHandle * handle,ParamBool live,ParamBool abf)4594 void param_lscm_begin(ParamHandle *handle, ParamBool live, ParamBool abf)
4595 {
4596 PHandle *phandle = (PHandle *)handle;
4597 PFace *f;
4598 int i;
4599
4600 param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
4601 phandle->state = PHANDLE_STATE_LSCM;
4602
4603 for (i = 0; i < phandle->ncharts; i++) {
4604 for (f = phandle->charts[i]->faces; f; f = f->nextlink) {
4605 p_face_backup_uvs(f);
4606 }
4607 p_chart_lscm_begin(phandle->charts[i], (PBool)live, (PBool)abf);
4608 }
4609 }
4610
param_lscm_solve(ParamHandle * handle)4611 void param_lscm_solve(ParamHandle *handle)
4612 {
4613 PHandle *phandle = (PHandle *)handle;
4614 PChart *chart;
4615 int i;
4616 PBool result;
4617
4618 param_assert(phandle->state == PHANDLE_STATE_LSCM);
4619
4620 for (i = 0; i < phandle->ncharts; i++) {
4621 chart = phandle->charts[i];
4622
4623 if (chart->u.lscm.context) {
4624 result = p_chart_lscm_solve(phandle, chart);
4625
4626 if (result && !(chart->flag & PCHART_HAS_PINS)) {
4627 p_chart_rotate_minimum_area(chart);
4628 }
4629 else if (result && chart->u.lscm.single_pin) {
4630 p_chart_rotate_fit_aabb(chart);
4631 p_chart_lscm_transform_single_pin(chart);
4632 }
4633
4634 if (!result || !(chart->flag & PCHART_HAS_PINS)) {
4635 p_chart_lscm_end(chart);
4636 }
4637 }
4638 }
4639 }
4640
param_lscm_end(ParamHandle * handle)4641 void param_lscm_end(ParamHandle *handle)
4642 {
4643 PHandle *phandle = (PHandle *)handle;
4644 int i;
4645
4646 param_assert(phandle->state == PHANDLE_STATE_LSCM);
4647
4648 for (i = 0; i < phandle->ncharts; i++) {
4649 p_chart_lscm_end(phandle->charts[i]);
4650 #if 0
4651 p_chart_complexify(phandle->charts[i]);
4652 #endif
4653 }
4654
4655 phandle->state = PHANDLE_STATE_CONSTRUCTED;
4656 }
4657
param_stretch_begin(ParamHandle * handle)4658 void param_stretch_begin(ParamHandle *handle)
4659 {
4660 PHandle *phandle = (PHandle *)handle;
4661 PChart *chart;
4662 PVert *v;
4663 PFace *f;
4664 int i;
4665
4666 param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
4667 phandle->state = PHANDLE_STATE_STRETCH;
4668
4669 phandle->rng = BLI_rng_new(31415926);
4670 phandle->blend = 0.0f;
4671
4672 for (i = 0; i < phandle->ncharts; i++) {
4673 chart = phandle->charts[i];
4674
4675 for (v = chart->verts; v; v = v->nextlink) {
4676 v->flag &= ~PVERT_PIN; /* don't use user-defined pins */
4677 }
4678
4679 p_stretch_pin_boundary(chart);
4680
4681 for (f = chart->faces; f; f = f->nextlink) {
4682 p_face_backup_uvs(f);
4683 f->u.area3d = p_face_area(f);
4684 }
4685 }
4686 }
4687
param_stretch_blend(ParamHandle * handle,float blend)4688 void param_stretch_blend(ParamHandle *handle, float blend)
4689 {
4690 PHandle *phandle = (PHandle *)handle;
4691
4692 param_assert(phandle->state == PHANDLE_STATE_STRETCH);
4693 phandle->blend = blend;
4694 }
4695
param_stretch_iter(ParamHandle * handle)4696 void param_stretch_iter(ParamHandle *handle)
4697 {
4698 PHandle *phandle = (PHandle *)handle;
4699 PChart *chart;
4700 int i;
4701
4702 param_assert(phandle->state == PHANDLE_STATE_STRETCH);
4703
4704 for (i = 0; i < phandle->ncharts; i++) {
4705 chart = phandle->charts[i];
4706 p_chart_stretch_minimize(chart, phandle->rng);
4707 }
4708 }
4709
param_stretch_end(ParamHandle * handle)4710 void param_stretch_end(ParamHandle *handle)
4711 {
4712 PHandle *phandle = (PHandle *)handle;
4713
4714 param_assert(phandle->state == PHANDLE_STATE_STRETCH);
4715 phandle->state = PHANDLE_STATE_CONSTRUCTED;
4716
4717 BLI_rng_free(phandle->rng);
4718 phandle->rng = NULL;
4719 }
4720
param_smooth_area(ParamHandle * handle)4721 void param_smooth_area(ParamHandle *handle)
4722 {
4723 PHandle *phandle = (PHandle *)handle;
4724 int i;
4725
4726 param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
4727
4728 for (i = 0; i < phandle->ncharts; i++) {
4729 PChart *chart = phandle->charts[i];
4730 PVert *v;
4731
4732 for (v = chart->verts; v; v = v->nextlink) {
4733 v->flag &= ~PVERT_PIN;
4734 }
4735
4736 p_smooth(chart);
4737 }
4738 }
4739
4740 /* don't pack, just rotate (used for better packing) */
param_pack_rotate(ParamHandle * handle,bool ignore_pinned)4741 static void param_pack_rotate(ParamHandle *handle, bool ignore_pinned)
4742 {
4743 PChart *chart;
4744 int i;
4745
4746 PHandle *phandle = (PHandle *)handle;
4747
4748 for (i = 0; i < phandle->ncharts; i++) {
4749 chart = phandle->charts[i];
4750
4751 if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
4752 continue;
4753 }
4754
4755 p_chart_rotate_fit_aabb(chart);
4756 }
4757 }
4758
param_pack(ParamHandle * handle,float margin,bool do_rotate,bool ignore_pinned)4759 void param_pack(ParamHandle *handle, float margin, bool do_rotate, bool ignore_pinned)
4760 {
4761 /* box packing variables */
4762 BoxPack *boxarray, *box;
4763 float tot_width, tot_height, scale;
4764
4765 PChart *chart;
4766 int i, unpacked = 0;
4767 float trans[2];
4768 double area = 0.0;
4769
4770 PHandle *phandle = (PHandle *)handle;
4771
4772 if (phandle->ncharts == 0) {
4773 return;
4774 }
4775
4776 /* this could be its own function */
4777 if (do_rotate) {
4778 param_pack_rotate(handle, ignore_pinned);
4779 }
4780
4781 if (phandle->aspx != phandle->aspy) {
4782 param_scale(handle, 1.0f / phandle->aspx, 1.0f / phandle->aspy);
4783 }
4784
4785 /* we may not use all these boxes */
4786 boxarray = MEM_mallocN(phandle->ncharts * sizeof(BoxPack), "BoxPack box");
4787
4788 for (i = 0; i < phandle->ncharts; i++) {
4789 chart = phandle->charts[i];
4790
4791 if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
4792 unpacked++;
4793 continue;
4794 }
4795
4796 box = boxarray + (i - unpacked);
4797
4798 p_chart_uv_bbox(chart, trans, chart->u.pack.size);
4799
4800 trans[0] = -trans[0];
4801 trans[1] = -trans[1];
4802
4803 p_chart_uv_translate(chart, trans);
4804
4805 box->w = chart->u.pack.size[0] + trans[0];
4806 box->h = chart->u.pack.size[1] + trans[1];
4807 box->index = i; /* warning this index skips PCHART_HAS_PINS boxes */
4808
4809 if (margin > 0.0f) {
4810 area += (double)sqrtf(box->w * box->h);
4811 }
4812 }
4813
4814 if (margin > 0.0f) {
4815 /* multiply the margin by the area to give predictable results not dependent on UV scale,
4816 * ...Without using the area running pack multiple times also gives a bad feedback loop.
4817 * multiply by 0.1 so the margin value from the UI can be from
4818 * 0.0 to 1.0 but not give a massive margin */
4819 margin = (margin * (float)area) * 0.1f;
4820 unpacked = 0;
4821 for (i = 0; i < phandle->ncharts; i++) {
4822 chart = phandle->charts[i];
4823
4824 if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
4825 unpacked++;
4826 continue;
4827 }
4828
4829 box = boxarray + (i - unpacked);
4830 trans[0] = margin;
4831 trans[1] = margin;
4832 p_chart_uv_translate(chart, trans);
4833 box->w += margin * 2;
4834 box->h += margin * 2;
4835 }
4836 }
4837
4838 BLI_box_pack_2d(boxarray, phandle->ncharts - unpacked, &tot_width, &tot_height);
4839
4840 if (tot_height > tot_width) {
4841 scale = 1.0f / tot_height;
4842 }
4843 else {
4844 scale = 1.0f / tot_width;
4845 }
4846
4847 for (i = 0; i < phandle->ncharts - unpacked; i++) {
4848 box = boxarray + i;
4849 trans[0] = box->x;
4850 trans[1] = box->y;
4851
4852 chart = phandle->charts[box->index];
4853 p_chart_uv_translate(chart, trans);
4854 p_chart_uv_scale(chart, scale);
4855 }
4856 MEM_freeN(boxarray);
4857
4858 if (phandle->aspx != phandle->aspy) {
4859 param_scale(handle, phandle->aspx, phandle->aspy);
4860 }
4861 }
4862
param_average(ParamHandle * handle,bool ignore_pinned)4863 void param_average(ParamHandle *handle, bool ignore_pinned)
4864 {
4865 PChart *chart;
4866 int i;
4867 float tot_uvarea = 0.0f, tot_facearea = 0.0f;
4868 float tot_fac, fac;
4869 float minv[2], maxv[2], trans[2];
4870 PHandle *phandle = (PHandle *)handle;
4871
4872 if (phandle->ncharts == 0) {
4873 return;
4874 }
4875
4876 for (i = 0; i < phandle->ncharts; i++) {
4877 PFace *f;
4878 chart = phandle->charts[i];
4879
4880 if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
4881 continue;
4882 }
4883
4884 chart->u.pack.area = 0.0f; /* 3d area */
4885 chart->u.pack.rescale = 0.0f; /* UV area, abusing rescale for tmp storage, oh well :/ */
4886
4887 for (f = chart->faces; f; f = f->nextlink) {
4888 chart->u.pack.area += p_face_area(f);
4889 chart->u.pack.rescale += fabsf(p_face_uv_area_signed(f));
4890 }
4891
4892 tot_facearea += chart->u.pack.area;
4893 tot_uvarea += chart->u.pack.rescale;
4894 }
4895
4896 if (tot_facearea == tot_uvarea || tot_facearea == 0.0f || tot_uvarea == 0.0f) {
4897 /* nothing to do */
4898 return;
4899 }
4900
4901 tot_fac = tot_facearea / tot_uvarea;
4902
4903 for (i = 0; i < phandle->ncharts; i++) {
4904 chart = phandle->charts[i];
4905
4906 if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
4907 continue;
4908 }
4909
4910 if (chart->u.pack.area != 0.0f && chart->u.pack.rescale != 0.0f) {
4911 fac = chart->u.pack.area / chart->u.pack.rescale;
4912
4913 /* Get the island center */
4914 p_chart_uv_bbox(chart, minv, maxv);
4915 trans[0] = (minv[0] + maxv[0]) / -2.0f;
4916 trans[1] = (minv[1] + maxv[1]) / -2.0f;
4917
4918 /* Move center to 0,0 */
4919 p_chart_uv_translate(chart, trans);
4920 p_chart_uv_scale(chart, sqrtf(fac / tot_fac));
4921
4922 /* Move to original center */
4923 trans[0] = -trans[0];
4924 trans[1] = -trans[1];
4925 p_chart_uv_translate(chart, trans);
4926 }
4927 }
4928 }
4929
param_scale(ParamHandle * handle,float x,float y)4930 void param_scale(ParamHandle *handle, float x, float y)
4931 {
4932 PHandle *phandle = (PHandle *)handle;
4933 PChart *chart;
4934 int i;
4935
4936 for (i = 0; i < phandle->ncharts; i++) {
4937 chart = phandle->charts[i];
4938 p_chart_uv_scale_xy(chart, x, y);
4939 }
4940 }
4941
param_flush(ParamHandle * handle)4942 void param_flush(ParamHandle *handle)
4943 {
4944 PHandle *phandle = (PHandle *)handle;
4945 PChart *chart;
4946 int i;
4947
4948 for (i = 0; i < phandle->ncharts; i++) {
4949 chart = phandle->charts[i];
4950
4951 if ((phandle->state == PHANDLE_STATE_LSCM) && !chart->u.lscm.context) {
4952 continue;
4953 }
4954
4955 if (phandle->blend == 0.0f) {
4956 p_flush_uvs(phandle, chart);
4957 }
4958 else {
4959 p_flush_uvs_blend(phandle, chart, phandle->blend);
4960 }
4961 }
4962 }
4963
param_flush_restore(ParamHandle * handle)4964 void param_flush_restore(ParamHandle *handle)
4965 {
4966 PHandle *phandle = (PHandle *)handle;
4967 PChart *chart;
4968 PFace *f;
4969 int i;
4970
4971 for (i = 0; i < phandle->ncharts; i++) {
4972 chart = phandle->charts[i];
4973
4974 for (f = chart->faces; f; f = f->nextlink) {
4975 p_face_restore_uvs(f);
4976 }
4977 }
4978 }
4979