412d68a5d4ddf16174572fd037ebd0f39c7ff552
[blender.git] / source / blender / editors / uvedit / uvedit_parametrizer.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * Contributor(s):
19  *
20  * ***** END GPL LICENSE BLOCK *****
21  */
22
23 /** \file blender/editors/uvedit/uvedit_parametrizer.c
24  *  \ingroup eduv
25  */
26
27 #include "MEM_guardedalloc.h"
28
29 #include "BLI_utildefines.h"
30 #include "BLI_alloca.h"
31 #include "BLI_memarena.h"
32 #include "BLI_math.h"
33 #include "BLI_rand.h"
34 #include "BLI_heap.h"
35 #include "BLI_boxpack_2d.h"
36 #include "BLI_convexhull_2d.h"
37 #include "BLI_polyfill_2d.h"
38 #include "BLI_polyfill_2d_beautify.h"
39
40 #include "uvedit_parametrizer.h"
41
42 #include <math.h>
43 #include <stdlib.h>
44 #include <stdio.h>
45 #include <string.h>
46
47 #include "BLI_sys_types.h"  /* for intptr_t support */
48
49 #include "eigen_capi.h"
50
51 /* Utils */
52
53 #if 0
54 #  define param_assert(condition)
55 #  define param_warning(message)
56 #  define param_test_equals_ptr(condition)
57 #  define param_test_equals_int(condition)
58 #else
59 #  define param_assert(condition) \
60                 if (!(condition)) \
61                         { /*printf("Assertion %s:%d\n", __FILE__, __LINE__); abort();*/ } (void)0
62 #  define param_warning(message) \
63                 { /*printf("Warning %s:%d: %s\n", __FILE__, __LINE__, message);*/ } (void)0
64 #  if 0
65 #    define param_test_equals_ptr(str, a, b) \
66                 if (a != b) \
67                         { /*printf("Equals %s => %p != %p\n", str, a, b);*/ } (void)0
68 #    define param_test_equals_int(str, a, b) \
69                 if (a != b) \
70                         { /*printf("Equals %s => %d != %d\n", str, a, b);*/ } (void)0
71 #  endif
72 #endif
73 typedef enum PBool {
74         P_TRUE = 1,
75         P_FALSE = 0
76 } PBool;
77
78 /* Special Purpose Hash */
79
80 typedef intptr_t PHashKey;
81
82 typedef struct PHashLink {
83         struct PHashLink *next;
84         PHashKey key;
85 } PHashLink;
86
87 typedef struct PHash {
88         PHashLink **list;
89         PHashLink **buckets;
90         int size, cursize, cursize_id;
91 } PHash;
92
93
94
95 struct PVert;
96 struct PEdge;
97 struct PFace;
98 struct PChart;
99 struct PHandle;
100
101 /* Simplices */
102
103 typedef struct PVert {
104         struct PVert *nextlink;
105
106         union PVertUnion {
107                 PHashKey key;           /* construct */
108                 int id;                 /* abf/lscm matrix index */
109                 float distortion;       /* area smoothing */
110                 HeapNode *heaplink;     /* edge collapsing */
111         } u;
112
113         struct PEdge *edge;
114         float co[3];
115         float uv[2];
116         unsigned char flag;
117
118 } PVert;
119
120 typedef struct PEdge {
121         struct PEdge *nextlink;
122
123         union PEdgeUnion {
124                 PHashKey key;                   /* construct */
125                 int id;                         /* abf matrix index */
126                 HeapNode *heaplink;             /* fill holes */
127                 struct PEdge *nextcollapse;     /* simplification */
128         } u;
129
130         struct PVert *vert;
131         struct PEdge *pair;
132         struct PEdge *next;
133         struct PFace *face;
134         float *orig_uv, old_uv[2];
135         unsigned short flag;
136
137 } PEdge;
138
139 typedef struct PFace {
140         struct PFace *nextlink;
141
142         union PFaceUnion {
143                 PHashKey key;           /* construct */
144                 int chart;              /* construct splitting*/
145                 float area3d;           /* stretch */
146                 int id;                 /* abf matrix index */
147         } u;
148
149         struct PEdge *edge;
150         unsigned char flag;
151 } PFace;
152
153 enum PVertFlag {
154         PVERT_PIN = 1,
155         PVERT_SELECT = 2,
156         PVERT_INTERIOR = 4,
157         PVERT_COLLAPSE = 8,
158         PVERT_SPLIT = 16
159 };
160
161 enum PEdgeFlag {
162         PEDGE_SEAM = 1,
163         PEDGE_VERTEX_SPLIT = 2,
164         PEDGE_PIN = 4,
165         PEDGE_SELECT = 8,
166         PEDGE_DONE = 16,
167         PEDGE_FILLED = 32,
168         PEDGE_COLLAPSE = 64,
169         PEDGE_COLLAPSE_EDGE = 128,
170         PEDGE_COLLAPSE_PAIR = 256
171 };
172
173 /* for flipping faces */
174 #define PEDGE_VERTEX_FLAGS (PEDGE_PIN)
175
176 enum PFaceFlag {
177         PFACE_CONNECTED = 1,
178         PFACE_FILLED = 2,
179         PFACE_COLLAPSE = 4
180 };
181
182 /* Chart */
183
184 typedef struct PChart {
185         PVert *verts;
186         PEdge *edges;
187         PFace *faces;
188         int nverts, nedges, nfaces;
189
190         PVert *collapsed_verts;
191         PEdge *collapsed_edges;
192         PFace *collapsed_faces;
193
194         union PChartUnion {
195                 struct PChartLscm {
196                         LinearSolver *context;
197                         float *abf_alpha;
198                         PVert *pin1, *pin2;
199                 } lscm;
200                 struct PChartPack {
201                         float rescale, area;
202                         float size[2] /* , trans[2] */;
203                 } pack;
204         } u;
205
206         unsigned char flag;
207         struct PHandle *handle;
208 } PChart;
209
210 enum PChartFlag {
211         PCHART_NOPACK = 1
212 };
213
214 enum PHandleState {
215         PHANDLE_STATE_ALLOCATED,
216         PHANDLE_STATE_CONSTRUCTED,
217         PHANDLE_STATE_LSCM,
218         PHANDLE_STATE_STRETCH
219 };
220
221 typedef struct PHandle {
222         enum PHandleState state;
223         MemArena *arena;
224         MemArena *polyfill_arena;
225         Heap *polyfill_heap;
226
227         PChart *construction_chart;
228         PHash *hash_verts;
229         PHash *hash_edges;
230         PHash *hash_faces;
231
232         PChart **charts;
233         int ncharts;
234
235         float aspx, aspy;
236
237         RNG *rng;
238         float blend;
239         char do_aspect;
240 } PHandle;
241
242 /* PHash
243  * - special purpose hash that keeps all its elements in a single linked list.
244  * - after construction, this hash is thrown away, and the list remains.
245  * - removing elements is not possible efficiently.
246  */
247
248 static int PHashSizes[] = {
249         1, 3, 5, 11, 17, 37, 67, 131, 257, 521, 1031, 2053, 4099, 8209,
250         16411, 32771, 65537, 131101, 262147, 524309, 1048583, 2097169,
251         4194319, 8388617, 16777259, 33554467, 67108879, 134217757, 268435459
252 };
253
254 #define PHASH_hash(ph, item) (((uintptr_t) (item)) % ((unsigned int) (ph)->cursize))
255 #define PHASH_edge(v1, v2)   (((v1) < (v2)) ? ((v1) * 39) ^ ((v2) * 31) : ((v1) * 31) ^ ((v2) * 39))
256
257 static PHash *phash_new(PHashLink **list, int sizehint)
258 {
259         PHash *ph = (PHash *)MEM_callocN(sizeof(PHash), "PHash");
260         ph->size = 0;
261         ph->cursize_id = 0;
262         ph->list = list;
263
264         while (PHashSizes[ph->cursize_id] < sizehint)
265                 ph->cursize_id++;
266
267         ph->cursize = PHashSizes[ph->cursize_id];
268         ph->buckets = (PHashLink **)MEM_callocN(ph->cursize * sizeof(*ph->buckets), "PHashBuckets");
269
270         return ph;
271 }
272
273 static void phash_delete(PHash *ph)
274 {
275         MEM_freeN(ph->buckets);
276         MEM_freeN(ph);
277 }
278
279 static int phash_size(PHash *ph)
280 {
281         return ph->size;
282 }
283
284 static void phash_insert(PHash *ph, PHashLink *link)
285 {
286         int size = ph->cursize;
287         uintptr_t hash = PHASH_hash(ph, link->key);
288         PHashLink *lookup = ph->buckets[hash];
289
290         if (lookup == NULL) {
291                 /* insert in front of the list */
292                 ph->buckets[hash] = link;
293                 link->next = *(ph->list);
294                 *(ph->list) = link;
295         }
296         else {
297                 /* insert after existing element */
298                 link->next = lookup->next;
299                 lookup->next = link;
300         }
301
302         ph->size++;
303
304         if (ph->size > (size * 3)) {
305                 PHashLink *next = NULL, *first = *(ph->list);
306
307                 ph->cursize = PHashSizes[++ph->cursize_id];
308                 MEM_freeN(ph->buckets);
309                 ph->buckets = (PHashLink **)MEM_callocN(ph->cursize * sizeof(*ph->buckets), "PHashBuckets");
310                 ph->size = 0;
311                 *(ph->list) = NULL;
312
313                 for (link = first; link; link = next) {
314                         next = link->next;
315                         phash_insert(ph, link);
316                 }
317         }
318 }
319
320 static PHashLink *phash_lookup(PHash *ph, PHashKey key)
321 {
322         PHashLink *link;
323         uintptr_t hash = PHASH_hash(ph, key);
324
325         for (link = ph->buckets[hash]; link; link = link->next)
326                 if (link->key == key)
327                         return link;
328                 else if (PHASH_hash(ph, link->key) != hash)
329                         return NULL;
330
331         return link;
332 }
333
334 static PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
335 {
336         uintptr_t hash = PHASH_hash(ph, key);
337
338         for (link = link->next; link; link = link->next)
339                 if (link->key == key)
340                         return link;
341                 else if (PHASH_hash(ph, link->key) != hash)
342                         return NULL;
343
344         return link;
345 }
346
347 /* Geometry */
348
349 static float p_vec_angle_cos(float *v1, float *v2, float *v3)
350 {
351         float d1[3], d2[3];
352
353         d1[0] = v1[0] - v2[0];
354         d1[1] = v1[1] - v2[1];
355         d1[2] = v1[2] - v2[2];
356
357         d2[0] = v3[0] - v2[0];
358         d2[1] = v3[1] - v2[1];
359         d2[2] = v3[2] - v2[2];
360
361         normalize_v3(d1);
362         normalize_v3(d2);
363
364         return d1[0] * d2[0] + d1[1] * d2[1] + d1[2] * d2[2];
365 }
366
367 static float p_vec_angle(float *v1, float *v2, float *v3)
368 {
369         float dot = p_vec_angle_cos(v1, v2, v3);
370
371         if (dot <= -1.0f)
372                 return (float)M_PI;
373         else if (dot >= 1.0f)
374                 return 0.0f;
375         else
376                 return acosf(dot);
377 }
378
379 static float p_vec2_angle(float *v1, float *v2, float *v3)
380 {
381         float u1[3], u2[3], u3[3];
382
383         u1[0] = v1[0]; u1[1] = v1[1]; u1[2] = 0.0f;
384         u2[0] = v2[0]; u2[1] = v2[1]; u2[2] = 0.0f;
385         u3[0] = v3[0]; u3[1] = v3[1]; u3[2] = 0.0f;
386
387         return p_vec_angle(u1, u2, u3);
388 }
389
390 static void p_triangle_angles(float *v1, float *v2, float *v3, float *a1, float *a2, float *a3)
391 {
392         *a1 = p_vec_angle(v3, v1, v2);
393         *a2 = p_vec_angle(v1, v2, v3);
394         *a3 = (float)M_PI - *a2 - *a1;
395 }
396
397 static void p_face_angles(PFace *f, float *a1, float *a2, float *a3)
398 {
399         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
400         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
401
402         p_triangle_angles(v1->co, v2->co, v3->co, a1, a2, a3);
403 }
404
405 static float p_face_area(PFace *f)
406 {
407         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
408         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
409
410         return area_tri_v3(v1->co, v2->co, v3->co);
411 }
412
413 static float p_area_signed(float *v1, float *v2, float *v3)
414 {
415         return 0.5f * (((v2[0] - v1[0]) * (v3[1] - v1[1])) -
416                        ((v3[0] - v1[0]) * (v2[1] - v1[1])));
417 }
418
419 static float p_face_uv_area_signed(PFace *f)
420 {
421         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
422         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
423
424         return 0.5f * (((v2->uv[0] - v1->uv[0]) * (v3->uv[1] - v1->uv[1])) -
425                        ((v3->uv[0] - v1->uv[0]) * (v2->uv[1] - v1->uv[1])));
426 }
427
428 static float p_edge_length(PEdge *e)
429 {
430         PVert *v1 = e->vert, *v2 = e->next->vert;
431         float d[3];
432
433         d[0] = v2->co[0] - v1->co[0];
434         d[1] = v2->co[1] - v1->co[1];
435         d[2] = v2->co[2] - v1->co[2];
436
437         return sqrtf(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
438 }
439
440 static float p_edge_uv_length(PEdge *e)
441 {
442         PVert *v1 = e->vert, *v2 = e->next->vert;
443         float d[3];
444
445         d[0] = v2->uv[0] - v1->uv[0];
446         d[1] = v2->uv[1] - v1->uv[1];
447
448         return sqrtf(d[0] * d[0] + d[1] * d[1]);
449 }
450
451 static void p_chart_uv_bbox(PChart *chart, float minv[2], float maxv[2])
452 {
453         PVert *v;
454
455         INIT_MINMAX2(minv, maxv);
456
457         for (v = chart->verts; v; v = v->nextlink) {
458                 minmax_v2v2_v2(minv, maxv, v->uv);
459         }
460 }
461
462 static void p_chart_uv_scale(PChart *chart, float scale)
463 {
464         PVert *v;
465
466         for (v = chart->verts; v; v = v->nextlink) {
467                 v->uv[0] *= scale;
468                 v->uv[1] *= scale;
469         }
470 }
471
472 static void p_chart_uv_scale_xy(PChart *chart, float x, float y)
473 {
474         PVert *v;
475
476         for (v = chart->verts; v; v = v->nextlink) {
477                 v->uv[0] *= x;
478                 v->uv[1] *= y;
479         }
480 }
481
482 static void p_chart_uv_translate(PChart *chart, float trans[2])
483 {
484         PVert *v;
485
486         for (v = chart->verts; v; v = v->nextlink) {
487                 v->uv[0] += trans[0];
488                 v->uv[1] += trans[1];
489         }
490 }
491
492 static void p_chart_uv_transform(PChart *chart, float mat[2][2])
493 {
494         PVert *v;
495
496         for (v = chart->verts; v; v = v->nextlink) {
497                 mul_m2v2(mat, v->uv);
498         }
499 }
500
501 static void p_chart_uv_to_array(PChart *chart, float (*points)[2])
502 {
503         PVert *v;
504         unsigned int i = 0;
505
506         for (v = chart->verts; v; v = v->nextlink) {
507                 copy_v2_v2(points[i++], v->uv);
508         }
509 }
510
511 static void UNUSED_FUNCTION(p_chart_uv_from_array)(PChart *chart, float (*points)[2])
512 {
513         PVert *v;
514         unsigned int i = 0;
515
516         for (v = chart->verts; v; v = v->nextlink) {
517                 copy_v2_v2(v->uv, points[i++]);
518         }
519 }
520
521
522 static PBool p_intersect_line_2d_dir(float *v1, float *dir1, float *v2, float *dir2, float *isect)
523 {
524         float lmbda, div;
525
526         div = dir2[0] * dir1[1] - dir2[1] * dir1[0];
527
528         if (div == 0.0f)
529                 return P_FALSE;
530
531         lmbda = ((v1[1] - v2[1]) * dir1[0] - (v1[0] - v2[0]) * dir1[1]) / div;
532         isect[0] = v1[0] + lmbda * dir2[0];
533         isect[1] = v1[1] + lmbda * dir2[1];
534
535         return P_TRUE;
536 }
537
538 #if 0
539 static PBool p_intersect_line_2d(float *v1, float *v2, float *v3, float *v4, float *isect)
540 {
541         float dir1[2], dir2[2];
542
543         dir1[0] = v4[0] - v3[0];
544         dir1[1] = v4[1] - v3[1];
545
546         dir2[0] = v2[0] - v1[0];
547         dir2[1] = v2[1] - v1[1];
548
549         if (!p_intersect_line_2d_dir(v1, dir1, v2, dir2, isect)) {
550                 /* parallel - should never happen in theory for polygon kernel, but
551                  * let's give a point nearby in case things go wrong */
552                 isect[0] = (v1[0] + v2[0]) * 0.5f;
553                 isect[1] = (v1[1] + v2[1]) * 0.5f;
554                 return P_FALSE;
555         }
556
557         return P_TRUE;
558 }
559 #endif
560
561 /* Topological Utilities */
562
563 static PEdge *p_wheel_edge_next(PEdge *e)
564 {
565         return e->next->next->pair;
566 }
567
568 static PEdge *p_wheel_edge_prev(PEdge *e)
569 {
570         return (e->pair) ? e->pair->next : NULL;
571 }
572
573 static PEdge *p_boundary_edge_next(PEdge *e)
574 {
575         return e->next->vert->edge;
576 }
577
578 static PEdge *p_boundary_edge_prev(PEdge *e)
579 {
580         PEdge *we = e, *last;
581
582         do {
583                 last = we;
584                 we = p_wheel_edge_next(we);
585         } while (we && (we != e));
586
587         return last->next->next;
588 }
589
590 static PBool p_vert_interior(PVert *v)
591 {
592         return (v->edge->pair != NULL);
593 }
594
595 static void p_face_flip(PFace *f)
596 {
597         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
598         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
599         int f1 = e1->flag, f2 = e2->flag, f3 = e3->flag;
600         float *orig_uv1 = e1->orig_uv, *orig_uv2 = e2->orig_uv, *orig_uv3 = e3->orig_uv;
601
602         e1->vert = v2;
603         e1->next = e3;
604         e1->orig_uv = orig_uv2;
605         e1->flag = (f1 & ~PEDGE_VERTEX_FLAGS) | (f2 & PEDGE_VERTEX_FLAGS);
606
607         e2->vert = v3;
608         e2->next = e1;
609         e2->orig_uv = orig_uv3;
610         e2->flag = (f2 & ~PEDGE_VERTEX_FLAGS) | (f3 & PEDGE_VERTEX_FLAGS);
611
612         e3->vert = v1;
613         e3->next = e2;
614         e3->orig_uv = orig_uv1;
615         e3->flag = (f3 & ~PEDGE_VERTEX_FLAGS) | (f1 & PEDGE_VERTEX_FLAGS);
616 }
617
618 #if 0
619 static void p_chart_topological_sanity_check(PChart *chart)
620 {
621         PVert *v;
622         PEdge *e;
623
624         for (v = chart->verts; v; v = v->nextlink)
625                 param_test_equals_ptr("v->edge->vert", v, v->edge->vert);
626
627         for (e = chart->edges; e; e = e->nextlink) {
628                 if (e->pair) {
629                         param_test_equals_ptr("e->pair->pair", e, e->pair->pair);
630                         param_test_equals_ptr("pair->vert", e->vert, e->pair->next->vert);
631                         param_test_equals_ptr("pair->next->vert", e->next->vert, e->pair->vert);
632                 }
633         }
634 }
635 #endif
636
637 /* Loading / Flushing */
638
639 static void p_vert_load_pin_select_uvs(PHandle *handle, PVert *v)
640 {
641         PEdge *e;
642         int nedges = 0, npins = 0;
643         float pinuv[2];
644
645         v->uv[0] = v->uv[1] = 0.0f;
646         pinuv[0] = pinuv[1] = 0.0f;
647         e = v->edge;
648         do {
649                 if (e->orig_uv) {
650                         if (e->flag & PEDGE_SELECT)
651                                 v->flag |= PVERT_SELECT;
652
653                         if (e->flag & PEDGE_PIN) {
654                                 pinuv[0] += e->orig_uv[0] * handle->aspx;
655                                 pinuv[1] += e->orig_uv[1] * handle->aspy;
656                                 npins++;
657                         }
658                         else {
659                                 v->uv[0] += e->orig_uv[0] * handle->aspx;
660                                 v->uv[1] += e->orig_uv[1] * handle->aspy;
661                         }
662
663                         nedges++;
664                 }
665
666                 e = p_wheel_edge_next(e);
667         } while (e && e != (v->edge));
668
669         if (npins > 0) {
670                 v->uv[0] = pinuv[0] / npins;
671                 v->uv[1] = pinuv[1] / npins;
672                 v->flag |= PVERT_PIN;
673         }
674         else if (nedges > 0) {
675                 v->uv[0] /= nedges;
676                 v->uv[1] /= nedges;
677         }
678 }
679
680 static void p_flush_uvs(PHandle *handle, PChart *chart)
681 {
682         PEdge *e;
683
684         for (e = chart->edges; e; e = e->nextlink) {
685                 if (e->orig_uv) {
686                         e->orig_uv[0] = e->vert->uv[0] / handle->aspx;
687                         e->orig_uv[1] = e->vert->uv[1] / handle->aspy;
688                 }
689         }
690 }
691
692 static void p_flush_uvs_blend(PHandle *handle, PChart *chart, float blend)
693 {
694         PEdge *e;
695         float invblend = 1.0f - blend;
696
697         for (e = chart->edges; e; e = e->nextlink) {
698                 if (e->orig_uv) {
699                         e->orig_uv[0] = blend * e->old_uv[0] + invblend * e->vert->uv[0] / handle->aspx;
700                         e->orig_uv[1] = blend * e->old_uv[1] + invblend * e->vert->uv[1] / handle->aspy;
701                 }
702         }
703 }
704
705 static void p_face_backup_uvs(PFace *f)
706 {
707         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
708
709         if (e1->orig_uv) {
710                 e1->old_uv[0] = e1->orig_uv[0];
711                 e1->old_uv[1] = e1->orig_uv[1];
712         }
713         if (e2->orig_uv) {
714                 e2->old_uv[0] = e2->orig_uv[0];
715                 e2->old_uv[1] = e2->orig_uv[1];
716         }
717         if (e3->orig_uv) {
718                 e3->old_uv[0] = e3->orig_uv[0];
719                 e3->old_uv[1] = e3->orig_uv[1];
720         }
721 }
722
723 static void p_face_restore_uvs(PFace *f)
724 {
725         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
726
727         if (e1->orig_uv) {
728                 e1->orig_uv[0] = e1->old_uv[0];
729                 e1->orig_uv[1] = e1->old_uv[1];
730         }
731         if (e2->orig_uv) {
732                 e2->orig_uv[0] = e2->old_uv[0];
733                 e2->orig_uv[1] = e2->old_uv[1];
734         }
735         if (e3->orig_uv) {
736                 e3->orig_uv[0] = e3->old_uv[0];
737                 e3->orig_uv[1] = e3->old_uv[1];
738         }
739 }
740
741 /* Construction (use only during construction, relies on u.key being set */
742
743 static PVert *p_vert_add(PHandle *handle, PHashKey key, const float co[3], PEdge *e)
744 {
745         PVert *v = (PVert *)BLI_memarena_alloc(handle->arena, sizeof(*v));
746         copy_v3_v3(v->co, co);
747
748         /* Sanity check, a single nan/inf point causes the entire result to be invalid.
749          * Note that values within the calculation may _become_ non-finite,
750          * so the rest of the code still needs to take this possibility into account. */
751         for (int i = 0; i < 3; i++) {
752                 if (UNLIKELY(!isfinite(v->co[i]))) {
753                         v->co[i] = 0.0f;
754                 }
755         }
756
757         v->u.key = key;
758         v->edge = e;
759         v->flag = 0;
760
761         phash_insert(handle->hash_verts, (PHashLink *)v);
762
763         return v;
764 }
765
766 static PVert *p_vert_lookup(PHandle *handle, PHashKey key, const float co[3], PEdge *e)
767 {
768         PVert *v = (PVert *)phash_lookup(handle->hash_verts, key);
769
770         if (v)
771                 return v;
772         else
773                 return p_vert_add(handle, key, co, e);
774 }
775
776 static PVert *p_vert_copy(PChart *chart, PVert *v)
777 {
778         PVert *nv = (PVert *)BLI_memarena_alloc(chart->handle->arena, sizeof(*nv));
779
780         copy_v3_v3(nv->co, v->co);
781         nv->uv[0] = v->uv[0];
782         nv->uv[1] = v->uv[1];
783         nv->u.key = v->u.key;
784         nv->edge = v->edge;
785         nv->flag = v->flag;
786
787         return nv;
788 }
789
790 static PEdge *p_edge_lookup(PHandle *handle, PHashKey *vkeys)
791 {
792         PHashKey key = PHASH_edge(vkeys[0], vkeys[1]);
793         PEdge *e = (PEdge *)phash_lookup(handle->hash_edges, key);
794
795         while (e) {
796                 if ((e->vert->u.key == vkeys[0]) && (e->next->vert->u.key == vkeys[1]))
797                         return e;
798                 else if ((e->vert->u.key == vkeys[1]) && (e->next->vert->u.key == vkeys[0]))
799                         return e;
800
801                 e = (PEdge *)phash_next(handle->hash_edges, key, (PHashLink *)e);
802         }
803
804         return NULL;
805 }
806
807 static int p_face_exists(ParamHandle *phandle, ParamKey *pvkeys, int i1, int i2, int i3)
808 {
809         PHandle *handle = (PHandle *)phandle;
810         PHashKey *vkeys = (PHashKey *)pvkeys;
811         PHashKey key = PHASH_edge(vkeys[i1], vkeys[i2]);
812         PEdge *e = (PEdge *)phash_lookup(handle->hash_edges, key);
813
814         while (e) {
815                 if ((e->vert->u.key == vkeys[i1]) && (e->next->vert->u.key == vkeys[i2])) {
816                         if (e->next->next->vert->u.key == vkeys[i3])
817                                 return P_TRUE;
818                 }
819                 else if ((e->vert->u.key == vkeys[i2]) && (e->next->vert->u.key == vkeys[i1])) {
820                         if (e->next->next->vert->u.key == vkeys[i3])
821                                 return P_TRUE;
822                 }
823
824                 e = (PEdge *)phash_next(handle->hash_edges, key, (PHashLink *)e);
825         }
826
827         return P_FALSE;
828 }
829
830 static PChart *p_chart_new(PHandle *handle)
831 {
832         PChart *chart = (PChart *)MEM_callocN(sizeof(*chart), "PChart");
833         chart->handle = handle;
834
835         return chart;
836 }
837
838 static void p_chart_delete(PChart *chart)
839 {
840         /* the actual links are free by memarena */
841         MEM_freeN(chart);
842 }
843
844 static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
845 {
846         float *uv1, *uv2, *uvp1, *uvp2;
847         float limit[2];
848
849         limit[0] = 0.00001;
850         limit[1] = 0.00001;
851
852         uv1 = e->orig_uv;
853         uv2 = e->next->orig_uv;
854
855         if (e->vert->u.key == ep->vert->u.key) {
856                 uvp1 = ep->orig_uv;
857                 uvp2 = ep->next->orig_uv;
858         }
859         else {
860                 uvp1 = ep->next->orig_uv;
861                 uvp2 = ep->orig_uv;
862         }
863
864         if ((fabsf(uv1[0] - uvp1[0]) > limit[0]) || (fabsf(uv1[1] - uvp1[1]) > limit[1])) {
865                 e->flag |= PEDGE_SEAM;
866                 ep->flag |= PEDGE_SEAM;
867                 return P_TRUE;
868         }
869         if ((fabsf(uv2[0] - uvp2[0]) > limit[0]) || (fabsf(uv2[1] - uvp2[1]) > limit[1])) {
870                 e->flag |= PEDGE_SEAM;
871                 ep->flag |= PEDGE_SEAM;
872                 return P_TRUE;
873         }
874
875         return P_FALSE;
876 }
877
878 static PBool p_edge_has_pair(PHandle *handle, PEdge *e, PEdge **pair, PBool impl)
879 {
880         PHashKey key;
881         PEdge *pe;
882         PVert *v1, *v2;
883         PHashKey key1 = e->vert->u.key;
884         PHashKey key2 = e->next->vert->u.key;
885
886         if (e->flag & PEDGE_SEAM)
887                 return P_FALSE;
888
889         key = PHASH_edge(key1, key2);
890         pe = (PEdge *)phash_lookup(handle->hash_edges, key);
891         *pair = NULL;
892
893         while (pe) {
894                 if (pe != e) {
895                         v1 = pe->vert;
896                         v2 = pe->next->vert;
897
898                         if (((v1->u.key == key1) && (v2->u.key == key2)) ||
899                             ((v1->u.key == key2) && (v2->u.key == key1)))
900                         {
901
902                                 /* don't connect seams and t-junctions */
903                                 if ((pe->flag & PEDGE_SEAM) || *pair ||
904                                     (impl && p_edge_implicit_seam(e, pe)))
905                                 {
906                                         *pair = NULL;
907                                         return P_FALSE;
908                                 }
909
910                                 *pair = pe;
911                         }
912                 }
913
914                 pe = (PEdge *)phash_next(handle->hash_edges, key, (PHashLink *)pe);
915         }
916
917         if (*pair && (e->vert == (*pair)->vert)) {
918                 if ((*pair)->next->pair || (*pair)->next->next->pair) {
919                         /* non unfoldable, maybe mobius ring or klein bottle */
920                         *pair = NULL;
921                         return P_FALSE;
922                 }
923         }
924
925         return (*pair != NULL);
926 }
927
928 static PBool p_edge_connect_pair(PHandle *handle, PEdge *e, PEdge ***stack, PBool impl)
929 {
930         PEdge *pair = NULL;
931
932         if (!e->pair && p_edge_has_pair(handle, e, &pair, impl)) {
933                 if (e->vert == pair->vert)
934                         p_face_flip(pair->face);
935
936                 e->pair = pair;
937                 pair->pair = e;
938
939                 if (!(pair->face->flag & PFACE_CONNECTED)) {
940                         **stack = pair;
941                         (*stack)++;
942                 }
943         }
944
945         return (e->pair != NULL);
946 }
947
948 static int p_connect_pairs(PHandle *handle, PBool impl)
949 {
950         PEdge **stackbase = MEM_mallocN(sizeof(*stackbase) * phash_size(handle->hash_faces), "Pstackbase");
951         PEdge **stack = stackbase;
952         PFace *f, *first;
953         PEdge *e, *e1, *e2;
954         PChart *chart = handle->construction_chart;
955         int ncharts = 0;
956
957         /* connect pairs, count edges, set vertex-edge pointer to a pairless edge */
958         for (first = chart->faces; first; first = first->nextlink) {
959                 if (first->flag & PFACE_CONNECTED)
960                         continue;
961
962                 *stack = first->edge;
963                 stack++;
964
965                 while (stack != stackbase) {
966                         stack--;
967                         e = *stack;
968                         e1 = e->next;
969                         e2 = e1->next;
970
971                         f = e->face;
972                         f->flag |= PFACE_CONNECTED;
973
974                         /* assign verts to charts so we can sort them later */
975                         f->u.chart = ncharts;
976
977                         if (!p_edge_connect_pair(handle, e, &stack, impl))
978                                 e->vert->edge = e;
979                         if (!p_edge_connect_pair(handle, e1, &stack, impl))
980                                 e1->vert->edge = e1;
981                         if (!p_edge_connect_pair(handle, e2, &stack, impl))
982                                 e2->vert->edge = e2;
983                 }
984
985                 ncharts++;
986         }
987
988         MEM_freeN(stackbase);
989
990         return ncharts;
991 }
992
993 static void p_split_vert(PChart *chart, PEdge *e)
994 {
995         PEdge *we, *lastwe = NULL;
996         PVert *v = e->vert;
997         PBool copy = P_TRUE;
998
999         if (e->flag & PEDGE_VERTEX_SPLIT)
1000                 return;
1001
1002         /* rewind to start */
1003         lastwe = e;
1004         for (we = p_wheel_edge_prev(e); we && (we != e); we = p_wheel_edge_prev(we))
1005                 lastwe = we;
1006
1007         /* go over all edges in wheel */
1008         for (we = lastwe; we; we = p_wheel_edge_next(we)) {
1009                 if (we->flag & PEDGE_VERTEX_SPLIT)
1010                         break;
1011
1012                 we->flag |= PEDGE_VERTEX_SPLIT;
1013
1014                 if (we == v->edge) {
1015                         /* found it, no need to copy */
1016                         copy = P_FALSE;
1017                         v->nextlink = chart->verts;
1018                         chart->verts = v;
1019                         chart->nverts++;
1020                 }
1021         }
1022
1023         if (copy) {
1024                 /* not found, copying */
1025                 v->flag |= PVERT_SPLIT;
1026                 v = p_vert_copy(chart, v);
1027                 v->flag |= PVERT_SPLIT;
1028
1029                 v->nextlink = chart->verts;
1030                 chart->verts = v;
1031                 chart->nverts++;
1032
1033                 v->edge = lastwe;
1034
1035                 we = lastwe;
1036                 do {
1037                         we->vert = v;
1038                         we = p_wheel_edge_next(we);
1039                 } while (we && (we != lastwe));
1040         }
1041 }
1042
1043 static PChart **p_split_charts(PHandle *handle, PChart *chart, int ncharts)
1044 {
1045         PChart **charts = MEM_mallocN(sizeof(*charts) * ncharts, "PCharts"), *nchart;
1046         PFace *f, *nextf;
1047         int i;
1048
1049         for (i = 0; i < ncharts; i++)
1050                 charts[i] = p_chart_new(handle);
1051
1052         f = chart->faces;
1053         while (f) {
1054                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1055                 nextf = f->nextlink;
1056
1057                 nchart = charts[f->u.chart];
1058
1059                 f->nextlink = nchart->faces;
1060                 nchart->faces = f;
1061                 e1->nextlink = nchart->edges;
1062                 nchart->edges = e1;
1063                 e2->nextlink = nchart->edges;
1064                 nchart->edges = e2;
1065                 e3->nextlink = nchart->edges;
1066                 nchart->edges = e3;
1067
1068                 nchart->nfaces++;
1069                 nchart->nedges += 3;
1070
1071                 p_split_vert(nchart, e1);
1072                 p_split_vert(nchart, e2);
1073                 p_split_vert(nchart, e3);
1074
1075                 f = nextf;
1076         }
1077
1078         return charts;
1079 }
1080
1081 static PFace *p_face_add(PHandle *handle)
1082 {
1083         PFace *f;
1084         PEdge *e1, *e2, *e3;
1085
1086         /* allocate */
1087         f = (PFace *)BLI_memarena_alloc(handle->arena, sizeof(*f));
1088         f->flag = 0;  /* init ! */
1089
1090         e1 = (PEdge *)BLI_memarena_alloc(handle->arena, sizeof(*e1));
1091         e2 = (PEdge *)BLI_memarena_alloc(handle->arena, sizeof(*e2));
1092         e3 = (PEdge *)BLI_memarena_alloc(handle->arena, sizeof(*e3));
1093
1094         /* set up edges */
1095         f->edge = e1;
1096         e1->face = e2->face = e3->face = f;
1097
1098         e1->next = e2;
1099         e2->next = e3;
1100         e3->next = e1;
1101
1102         e1->pair = NULL;
1103         e2->pair = NULL;
1104         e3->pair = NULL;
1105
1106         e1->flag = 0;
1107         e2->flag = 0;
1108         e3->flag = 0;
1109
1110         return f;
1111 }
1112
1113 static PFace *p_face_add_construct(PHandle *handle, ParamKey key, ParamKey *vkeys,
1114                                    float *co[4], float *uv[4], int i1, int i2, int i3,
1115                                    ParamBool *pin, ParamBool *select)
1116 {
1117         PFace *f = p_face_add(handle);
1118         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1119
1120         e1->vert = p_vert_lookup(handle, vkeys[i1], co[i1], e1);
1121         e2->vert = p_vert_lookup(handle, vkeys[i2], co[i2], e2);
1122         e3->vert = p_vert_lookup(handle, vkeys[i3], co[i3], e3);
1123
1124         e1->orig_uv = uv[i1];
1125         e2->orig_uv = uv[i2];
1126         e3->orig_uv = uv[i3];
1127
1128         if (pin) {
1129                 if (pin[i1]) e1->flag |= PEDGE_PIN;
1130                 if (pin[i2]) e2->flag |= PEDGE_PIN;
1131                 if (pin[i3]) e3->flag |= PEDGE_PIN;
1132         }
1133
1134         if (select) {
1135                 if (select[i1]) e1->flag |= PEDGE_SELECT;
1136                 if (select[i2]) e2->flag |= PEDGE_SELECT;
1137                 if (select[i3]) e3->flag |= PEDGE_SELECT;
1138         }
1139
1140         /* insert into hash */
1141         f->u.key = key;
1142         phash_insert(handle->hash_faces, (PHashLink *)f);
1143
1144         e1->u.key = PHASH_edge(vkeys[i1], vkeys[i2]);
1145         e2->u.key = PHASH_edge(vkeys[i2], vkeys[i3]);
1146         e3->u.key = PHASH_edge(vkeys[i3], vkeys[i1]);
1147
1148         phash_insert(handle->hash_edges, (PHashLink *)e1);
1149         phash_insert(handle->hash_edges, (PHashLink *)e2);
1150         phash_insert(handle->hash_edges, (PHashLink *)e3);
1151
1152         return f;
1153 }
1154
1155 static PFace *p_face_add_fill(PChart *chart, PVert *v1, PVert *v2, PVert *v3)
1156 {
1157         PFace *f = p_face_add(chart->handle);
1158         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1159
1160         e1->vert = v1;
1161         e2->vert = v2;
1162         e3->vert = v3;
1163
1164         e1->orig_uv = e2->orig_uv = e3->orig_uv = NULL;
1165
1166         f->nextlink = chart->faces;
1167         chart->faces = f;
1168         e1->nextlink = chart->edges;
1169         chart->edges = e1;
1170         e2->nextlink = chart->edges;
1171         chart->edges = e2;
1172         e3->nextlink = chart->edges;
1173         chart->edges = e3;
1174
1175         chart->nfaces++;
1176         chart->nedges += 3;
1177
1178         return f;
1179 }
1180
1181 static PBool p_quad_split_direction(PHandle *handle, float **co, PHashKey *vkeys)
1182 {
1183         /* slight bias to prefer one edge over the other in case they are equal, so
1184          * that in symmetric models we choose the same split direction instead of
1185          * depending on floating point errors to decide */
1186         float bias = 1.0f + 1e-6f;
1187         float fac = len_v3v3(co[0], co[2]) * bias - len_v3v3(co[1], co[3]);
1188         PBool dir = (fac <= 0.0f);
1189
1190         /* the face exists check is there because of a special case: when
1191          * two quads share three vertices, they can each be split into two
1192          * triangles, resulting in two identical triangles. for example in
1193          * suzanne's nose. */
1194         if (dir) {
1195                 if (p_face_exists(handle, vkeys, 0, 1, 2) || p_face_exists(handle, vkeys, 0, 2, 3))
1196                         return !dir;
1197         }
1198         else {
1199                 if (p_face_exists(handle, vkeys, 0, 1, 3) || p_face_exists(handle, vkeys, 1, 2, 3))
1200                         return !dir;
1201         }
1202
1203         return dir;
1204 }
1205
1206 /* Construction: boundary filling */
1207
1208 static void p_chart_boundaries(PChart *chart, int *nboundaries, PEdge **outer)
1209 {
1210         PEdge *e, *be;
1211         float len, maxlen = -1.0;
1212
1213         if (nboundaries)
1214                 *nboundaries = 0;
1215         if (outer)
1216                 *outer = NULL;
1217
1218         for (e = chart->edges; e; e = e->nextlink) {
1219                 if (e->pair || (e->flag & PEDGE_DONE))
1220                         continue;
1221
1222                 if (nboundaries)
1223                         (*nboundaries)++;
1224
1225                 len = 0.0f;
1226
1227                 be = e;
1228                 do {
1229                         be->flag |= PEDGE_DONE;
1230                         len += p_edge_length(be);
1231                         be = be->next->vert->edge;
1232                 } while (be != e);
1233
1234                 if (outer && (len > maxlen)) {
1235                         *outer = e;
1236                         maxlen = len;
1237                 }
1238         }
1239
1240         for (e = chart->edges; e; e = e->nextlink)
1241                 e->flag &= ~PEDGE_DONE;
1242 }
1243
1244 static float p_edge_boundary_angle(PEdge *e)
1245 {
1246         PEdge *we;
1247         PVert *v, *v1, *v2;
1248         float angle;
1249         int n = 0;
1250
1251         v = e->vert;
1252
1253         /* concave angle check -- could be better */
1254         angle = M_PI;
1255
1256         we = v->edge;
1257         do {
1258                 v1 = we->next->vert;
1259                 v2 = we->next->next->vert;
1260                 angle -= p_vec_angle(v1->co, v->co, v2->co);
1261
1262                 we = we->next->next->pair;
1263                 n++;
1264         } while (we && (we != v->edge));
1265
1266         return angle;
1267 }
1268
1269 static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
1270 {
1271         PEdge *e, *e1, *e2;
1272
1273         PFace *f;
1274         struct Heap *heap = BLI_heap_new();
1275         float angle;
1276
1277         e = be;
1278         do {
1279                 angle = p_edge_boundary_angle(e);
1280                 e->u.heaplink = BLI_heap_insert(heap, angle, e);
1281
1282                 e = p_boundary_edge_next(e);
1283         } while (e != be);
1284
1285         if (nedges == 2) {
1286                 /* no real boundary, but an isolated seam */
1287                 e = be->next->vert->edge;
1288                 e->pair = be;
1289                 be->pair = e;
1290
1291                 BLI_heap_remove(heap, e->u.heaplink);
1292                 BLI_heap_remove(heap, be->u.heaplink);
1293         }
1294         else {
1295                 while (nedges > 2) {
1296                         PEdge *ne, *ne1, *ne2;
1297
1298                         e = (PEdge *)BLI_heap_pop_min(heap);
1299
1300                         e1 = p_boundary_edge_prev(e);
1301                         e2 = p_boundary_edge_next(e);
1302
1303                         BLI_heap_remove(heap, e1->u.heaplink);
1304                         BLI_heap_remove(heap, e2->u.heaplink);
1305                         e->u.heaplink = e1->u.heaplink = e2->u.heaplink = NULL;
1306
1307                         e->flag |= PEDGE_FILLED;
1308                         e1->flag |= PEDGE_FILLED;
1309
1310                         f = p_face_add_fill(chart, e->vert, e1->vert, e2->vert);
1311                         f->flag |= PFACE_FILLED;
1312
1313                         ne = f->edge->next->next;
1314                         ne1 = f->edge;
1315                         ne2 = f->edge->next;
1316
1317                         ne->flag = ne1->flag = ne2->flag = PEDGE_FILLED;
1318
1319                         e->pair = ne;
1320                         ne->pair = e;
1321                         e1->pair = ne1;
1322                         ne1->pair = e1;
1323
1324                         ne->vert = e2->vert;
1325                         ne1->vert = e->vert;
1326                         ne2->vert = e1->vert;
1327
1328                         if (nedges == 3) {
1329                                 e2->pair = ne2;
1330                                 ne2->pair = e2;
1331                         }
1332                         else {
1333                                 ne2->vert->edge = ne2;
1334
1335                                 ne2->u.heaplink = BLI_heap_insert(heap, p_edge_boundary_angle(ne2), ne2);
1336                                 e2->u.heaplink = BLI_heap_insert(heap, p_edge_boundary_angle(e2), e2);
1337                         }
1338
1339                         nedges--;
1340                 }
1341         }
1342
1343         BLI_heap_free(heap, NULL);
1344 }
1345
1346 static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
1347 {
1348         PEdge *e, *be; /* *enext - as yet unused */
1349         int nedges;
1350
1351         for (e = chart->edges; e; e = e->nextlink) {
1352                 /* enext = e->nextlink; - as yet unused */
1353
1354                 if (e->pair || (e->flag & PEDGE_FILLED))
1355                         continue;
1356
1357                 nedges = 0;
1358                 be = e;
1359                 do {
1360                         be->flag |= PEDGE_FILLED;
1361                         be = be->next->vert->edge;
1362                         nedges++;
1363                 } while (be != e);
1364
1365                 if (e != outer)
1366                         p_chart_fill_boundary(chart, e, nedges);
1367         }
1368 }
1369
1370 #if 0
1371 /* Polygon kernel for inserting uv's non overlapping */
1372
1373 static int p_polygon_point_in(float *cp1, float *cp2, float *p)
1374 {
1375         if ((cp1[0] == p[0]) && (cp1[1] == p[1]))
1376                 return 2;
1377         else if ((cp2[0] == p[0]) && (cp2[1] == p[1]))
1378                 return 3;
1379         else
1380                 return (p_area_signed(cp1, cp2, p) >= 0.0f);
1381 }
1382
1383 static void p_polygon_kernel_clip(float (*oldpoints)[2], int noldpoints, float (*newpoints)[2], int *nnewpoints, float *cp1, float *cp2)
1384 {
1385         float *p2, *p1, isect[2];
1386         int i, p2in, p1in;
1387
1388         p1 = oldpoints[noldpoints - 1];
1389         p1in = p_polygon_point_in(cp1, cp2, p1);
1390         *nnewpoints = 0;
1391
1392         for (i = 0; i < noldpoints; i++) {
1393                 p2 = oldpoints[i];
1394                 p2in = p_polygon_point_in(cp1, cp2, p2);
1395
1396                 if ((p2in >= 2) || (p1in && p2in)) {
1397                         newpoints[*nnewpoints][0] = p2[0];
1398                         newpoints[*nnewpoints][1] = p2[1];
1399                         (*nnewpoints)++;
1400                 }
1401                 else if (p1in && !p2in) {
1402                         if (p1in != 3) {
1403                                 p_intersect_line_2d(p1, p2, cp1, cp2, isect);
1404                                 newpoints[*nnewpoints][0] = isect[0];
1405                                 newpoints[*nnewpoints][1] = isect[1];
1406                                 (*nnewpoints)++;
1407                         }
1408                 }
1409                 else if (!p1in && p2in) {
1410                         p_intersect_line_2d(p1, p2, cp1, cp2, isect);
1411                         newpoints[*nnewpoints][0] = isect[0];
1412                         newpoints[*nnewpoints][1] = isect[1];
1413                         (*nnewpoints)++;
1414
1415                         newpoints[*nnewpoints][0] = p2[0];
1416                         newpoints[*nnewpoints][1] = p2[1];
1417                         (*nnewpoints)++;
1418                 }
1419
1420                 p1in = p2in;
1421                 p1 = p2;
1422         }
1423 }
1424
1425 static void p_polygon_kernel_center(float (*points)[2], int npoints, float *center)
1426 {
1427         int i, size, nnewpoints = npoints;
1428         float (*oldpoints)[2], (*newpoints)[2], *p1, *p2;
1429
1430         size = npoints * 3;
1431         oldpoints = MEM_mallocN(sizeof(float) * 2 * size, "PPolygonOldPoints");
1432         newpoints = MEM_mallocN(sizeof(float) * 2 * size, "PPolygonNewPoints");
1433
1434         memcpy(oldpoints, points, sizeof(float) * 2 * npoints);
1435
1436         for (i = 0; i < npoints; i++) {
1437                 p1 = points[i];
1438                 p2 = points[(i + 1) % npoints];
1439                 p_polygon_kernel_clip(oldpoints, nnewpoints, newpoints, &nnewpoints, p1, p2);
1440
1441                 if (nnewpoints == 0) {
1442                         /* degenerate case, use center of original polygon */
1443                         memcpy(oldpoints, points, sizeof(float) * 2 * npoints);
1444                         nnewpoints = npoints;
1445                         break;
1446                 }
1447                 else if (nnewpoints == 1) {
1448                         /* degenerate case, use remaining point */
1449                         center[0] = newpoints[0][0];
1450                         center[1] = newpoints[0][1];
1451
1452                         MEM_freeN(oldpoints);
1453                         MEM_freeN(newpoints);
1454
1455                         return;
1456                 }
1457
1458                 if (nnewpoints * 2 > size) {
1459                         size *= 2;
1460                         MEM_freeN(oldpoints);
1461                         oldpoints = MEM_mallocN(sizeof(float) * 2 * size, "oldpoints");
1462                         memcpy(oldpoints, newpoints, sizeof(float) * 2 * nnewpoints);
1463                         MEM_freeN(newpoints);
1464                         newpoints = MEM_mallocN(sizeof(float) * 2 * size, "newpoints");
1465                 }
1466                 else {
1467                         float (*sw_points)[2] = oldpoints;
1468                         oldpoints = newpoints;
1469                         newpoints = sw_points;
1470                 }
1471         }
1472
1473         center[0] = center[1] = 0.0f;
1474
1475         for (i = 0; i < nnewpoints; i++) {
1476                 center[0] += oldpoints[i][0];
1477                 center[1] += oldpoints[i][1];
1478         }
1479
1480         center[0] /= nnewpoints;
1481         center[1] /= nnewpoints;
1482
1483         MEM_freeN(oldpoints);
1484         MEM_freeN(newpoints);
1485 }
1486 #endif
1487
1488 #if 0
1489 /* Edge Collapser */
1490
1491 int NCOLLAPSE = 1;
1492 int NCOLLAPSEX = 0;
1493
1494 static float p_vert_cotan(float *v1, float *v2, float *v3)
1495 {
1496         float a[3], b[3], c[3], clen;
1497
1498         sub_v3_v3v3(a, v2, v1);
1499         sub_v3_v3v3(b, v3, v1);
1500         cross_v3_v3v3(c, a, b);
1501
1502         clen = len_v3(c);
1503
1504         if (clen == 0.0f)
1505                 return 0.0f;
1506
1507         return dot_v3v3(a, b) / clen;
1508 }
1509
1510 static PBool p_vert_flipped_wheel_triangle(PVert *v)
1511 {
1512         PEdge *e = v->edge;
1513
1514         do {
1515                 if (p_face_uv_area_signed(e->face) < 0.0f)
1516                         return P_TRUE;
1517
1518                 e = p_wheel_edge_next(e);
1519         } while (e && (e != v->edge));
1520
1521         return P_FALSE;
1522 }
1523
1524 static PBool p_vert_map_harmonic_weights(PVert *v)
1525 {
1526         float weightsum, positionsum[2], olduv[2];
1527
1528         weightsum = 0.0f;
1529         positionsum[0] = positionsum[1] = 0.0f;
1530
1531         if (p_vert_interior(v)) {
1532                 PEdge *e = v->edge;
1533
1534                 do {
1535                         float t1, t2, weight;
1536                         PVert *v1, *v2;
1537
1538                         v1 = e->next->vert;
1539                         v2 = e->next->next->vert;
1540                         t1 = p_vert_cotan(v2->co, e->vert->co, v1->co);
1541
1542                         v1 = e->pair->next->vert;
1543                         v2 = e->pair->next->next->vert;
1544                         t2 = p_vert_cotan(v2->co, e->pair->vert->co, v1->co);
1545
1546                         weight = 0.5f * (t1 + t2);
1547                         weightsum += weight;
1548                         positionsum[0] += weight * e->pair->vert->uv[0];
1549                         positionsum[1] += weight * e->pair->vert->uv[1];
1550
1551                         e = p_wheel_edge_next(e);
1552                 } while (e && (e != v->edge));
1553         }
1554         else {
1555                 PEdge *e = v->edge;
1556
1557                 do {
1558                         float t1, t2;
1559                         PVert *v1, *v2;
1560
1561                         v2 = e->next->vert;
1562                         v1 = e->next->next->vert;
1563
1564                         t1 = p_vert_cotan(v1->co, v->co, v2->co);
1565                         t2 = p_vert_cotan(v2->co, v->co, v1->co);
1566
1567                         weightsum += t1 + t2;
1568                         positionsum[0] += (v2->uv[1] - v1->uv[1]) + (t1 * v2->uv[0] + t2 * v1->uv[0]);
1569                         positionsum[1] += (v1->uv[0] - v2->uv[0]) + (t1 * v2->uv[1] + t2 * v1->uv[1]);
1570
1571                         e = p_wheel_edge_next(e);
1572                 } while (e && (e != v->edge));
1573         }
1574
1575         if (weightsum != 0.0f) {
1576                 weightsum = 1.0f / weightsum;
1577                 positionsum[0] *= weightsum;
1578                 positionsum[1] *= weightsum;
1579         }
1580
1581         olduv[0] = v->uv[0];
1582         olduv[1] = v->uv[1];
1583         v->uv[0] = positionsum[0];
1584         v->uv[1] = positionsum[1];
1585
1586         if (p_vert_flipped_wheel_triangle(v)) {
1587                 v->uv[0] = olduv[0];
1588                 v->uv[1] = olduv[1];
1589
1590                 return P_FALSE;
1591         }
1592
1593         return P_TRUE;
1594 }
1595
1596 static void p_vert_harmonic_insert(PVert *v)
1597 {
1598         PEdge *e;
1599
1600         if (!p_vert_map_harmonic_weights(v)) {
1601                 /* do polygon kernel center insertion: this is quite slow, but should
1602                  * only be needed for 0.01 % of verts or so, when insert with harmonic
1603                  * weights fails */
1604
1605                 int npoints = 0, i;
1606                 float (*points)[2];
1607
1608                 e = v->edge;
1609                 do {
1610                         npoints++;
1611                         e = p_wheel_edge_next(e);
1612                 } while (e && (e != v->edge));
1613
1614                 if (e == NULL)
1615                         npoints++;
1616
1617                 points = MEM_mallocN(sizeof(float) * 2 * npoints, "PHarmonicPoints");
1618
1619                 e = v->edge;
1620                 i = 0;
1621                 do {
1622                         PEdge *nexte = p_wheel_edge_next(e);
1623
1624                         points[i][0] = e->next->vert->uv[0];
1625                         points[i][1] = e->next->vert->uv[1];
1626
1627                         if (nexte == NULL) {
1628                                 i++;
1629                                 points[i][0] = e->next->next->vert->uv[0];
1630                                 points[i][1] = e->next->next->vert->uv[1];
1631                                 break;
1632                         }
1633
1634                         e = nexte;
1635                         i++;
1636                 } while (e != v->edge);
1637
1638                 p_polygon_kernel_center(points, npoints, v->uv);
1639
1640                 MEM_freeN(points);
1641         }
1642
1643         e = v->edge;
1644         do {
1645                 if (!(e->next->vert->flag & PVERT_PIN))
1646                         p_vert_map_harmonic_weights(e->next->vert);
1647                 e = p_wheel_edge_next(e);
1648         } while (e && (e != v->edge));
1649
1650         p_vert_map_harmonic_weights(v);
1651 }
1652
1653 static void p_vert_fix_edge_pointer(PVert *v)
1654 {
1655         PEdge *start = v->edge;
1656
1657         /* set v->edge pointer to the edge with no pair, if there is one */
1658         while (v->edge->pair) {
1659                 v->edge = p_wheel_edge_prev(v->edge);
1660
1661                 if (v->edge == start)
1662                         break;
1663         }
1664 }
1665
1666 static void p_collapsing_verts(PEdge *edge, PEdge *pair, PVert **newv, PVert **keepv)
1667 {
1668         /* the two vertices that are involved in the collapse */
1669         if (edge) {
1670                 *newv = edge->vert;
1671                 *keepv = edge->next->vert;
1672         }
1673         else {
1674                 *newv = pair->next->vert;
1675                 *keepv = pair->vert;
1676         }
1677 }
1678
1679 static void p_collapse_edge(PEdge *edge, PEdge *pair)
1680 {
1681         PVert *oldv, *keepv;
1682         PEdge *e;
1683
1684         p_collapsing_verts(edge, pair, &oldv, &keepv);
1685
1686         /* change e->vert pointers from old vertex to the target vertex */
1687         e = oldv->edge;
1688         do {
1689                 if ((e != edge) && !(pair && pair->next == e))
1690                         e->vert = keepv;
1691
1692                 e = p_wheel_edge_next(e);
1693         } while (e && (e != oldv->edge));
1694
1695         /* set keepv->edge pointer */
1696         if ((edge && (keepv->edge == edge->next)) || (keepv->edge == pair)) {
1697                 if (edge && edge->next->pair)
1698                         keepv->edge = edge->next->pair->next;
1699                 else if (pair && pair->next->next->pair)
1700                         keepv->edge = pair->next->next->pair;
1701                 else if (edge && edge->next->next->pair)
1702                         keepv->edge = edge->next->next->pair;
1703                 else
1704                         keepv->edge = pair->next->pair->next;
1705         }
1706
1707         /* update pairs and v->edge pointers */
1708         if (edge) {
1709                 PEdge *e1 = edge->next, *e2 = e1->next;
1710
1711                 if (e1->pair)
1712                         e1->pair->pair = e2->pair;
1713
1714                 if (e2->pair) {
1715                         e2->pair->pair = e1->pair;
1716                         e2->vert->edge = p_wheel_edge_prev(e2);
1717                 }
1718                 else
1719                         e2->vert->edge = p_wheel_edge_next(e2);
1720
1721                 p_vert_fix_edge_pointer(e2->vert);
1722         }
1723
1724         if (pair) {
1725                 PEdge *e1 = pair->next, *e2 = e1->next;
1726
1727                 if (e1->pair)
1728                         e1->pair->pair = e2->pair;
1729
1730                 if (e2->pair) {
1731                         e2->pair->pair = e1->pair;
1732                         e2->vert->edge = p_wheel_edge_prev(e2);
1733                 }
1734                 else
1735                         e2->vert->edge = p_wheel_edge_next(e2);
1736
1737                 p_vert_fix_edge_pointer(e2->vert);
1738         }
1739
1740         p_vert_fix_edge_pointer(keepv);
1741
1742         /* mark for move to collapsed list later */
1743         oldv->flag |= PVERT_COLLAPSE;
1744
1745         if (edge) {
1746                 PFace *f = edge->face;
1747                 PEdge *e1 = edge->next, *e2 = e1->next;
1748
1749                 f->flag |= PFACE_COLLAPSE;
1750                 edge->flag |= PEDGE_COLLAPSE;
1751                 e1->flag |= PEDGE_COLLAPSE;
1752                 e2->flag |= PEDGE_COLLAPSE;
1753         }
1754
1755         if (pair) {
1756                 PFace *f = pair->face;
1757                 PEdge *e1 = pair->next, *e2 = e1->next;
1758
1759                 f->flag |= PFACE_COLLAPSE;
1760                 pair->flag |= PEDGE_COLLAPSE;
1761                 e1->flag |= PEDGE_COLLAPSE;
1762                 e2->flag |= PEDGE_COLLAPSE;
1763         }
1764 }
1765
1766 static void p_split_vertex(PEdge *edge, PEdge *pair)
1767 {
1768         PVert *newv, *keepv;
1769         PEdge *e;
1770
1771         p_collapsing_verts(edge, pair, &newv, &keepv);
1772
1773         /* update edge pairs */
1774         if (edge) {
1775                 PEdge *e1 = edge->next, *e2 = e1->next;
1776
1777                 if (e1->pair)
1778                         e1->pair->pair = e1;
1779                 if (e2->pair)
1780                         e2->pair->pair = e2;
1781
1782                 e2->vert->edge = e2;
1783                 p_vert_fix_edge_pointer(e2->vert);
1784                 keepv->edge = e1;
1785         }
1786
1787         if (pair) {
1788                 PEdge *e1 = pair->next, *e2 = e1->next;
1789
1790                 if (e1->pair)
1791                         e1->pair->pair = e1;
1792                 if (e2->pair)
1793                         e2->pair->pair = e2;
1794
1795                 e2->vert->edge = e2;
1796                 p_vert_fix_edge_pointer(e2->vert);
1797                 keepv->edge = pair;
1798         }
1799
1800         p_vert_fix_edge_pointer(keepv);
1801
1802         /* set e->vert pointers to restored vertex */
1803         e = newv->edge;
1804         do {
1805                 e->vert = newv;
1806                 e = p_wheel_edge_next(e);
1807         } while (e && (e != newv->edge));
1808 }
1809
1810 static PBool p_collapse_allowed_topologic(PEdge *edge, PEdge *pair)
1811 {
1812         PVert *oldv, *keepv;
1813
1814         p_collapsing_verts(edge, pair, &oldv, &keepv);
1815
1816         /* boundary edges */
1817         if (!edge || !pair) {
1818                 /* avoid collapsing chart into an edge */
1819                 if (edge && !edge->next->pair && !edge->next->next->pair)
1820                         return P_FALSE;
1821                 else if (pair && !pair->next->pair && !pair->next->next->pair)
1822                         return P_FALSE;
1823         }
1824         /* avoid merging two boundaries (oldv and keepv are on the 'other side' of
1825          * the chart) */
1826         else if (!p_vert_interior(oldv) && !p_vert_interior(keepv))
1827                 return P_FALSE;
1828
1829         return P_TRUE;
1830 }
1831
1832 static PBool p_collapse_normal_flipped(float *v1, float *v2, float *vold, float *vnew)
1833 {
1834         float nold[3], nnew[3], sub1[3], sub2[3];
1835
1836         sub_v3_v3v3(sub1, vold, v1);
1837         sub_v3_v3v3(sub2, vold, v2);
1838         cross_v3_v3v3(nold, sub1, sub2);
1839
1840         sub_v3_v3v3(sub1, vnew, v1);
1841         sub_v3_v3v3(sub2, vnew, v2);
1842         cross_v3_v3v3(nnew, sub1, sub2);
1843
1844         return (dot_v3v3(nold, nnew) <= 0.0f);
1845 }
1846
1847 static PBool p_collapse_allowed_geometric(PEdge *edge, PEdge *pair)
1848 {
1849         PVert *oldv, *keepv;
1850         PEdge *e;
1851         float angulardefect, angle;
1852
1853         p_collapsing_verts(edge, pair, &oldv, &keepv);
1854
1855         angulardefect = 2 * M_PI;
1856
1857         e = oldv->edge;
1858         do {
1859                 float a[3], b[3], minangle, maxangle;
1860                 PEdge *e1 = e->next, *e2 = e1->next;
1861                 PVert *v1 = e1->vert, *v2 = e2->vert;
1862                 int i;
1863
1864                 angle = p_vec_angle(v1->co, oldv->co, v2->co);
1865                 angulardefect -= angle;
1866
1867                 /* skip collapsing faces */
1868                 if (v1 == keepv || v2 == keepv) {
1869                         e = p_wheel_edge_next(e);
1870                         continue;
1871                 }
1872
1873                 if (p_collapse_normal_flipped(v1->co, v2->co, oldv->co, keepv->co))
1874                         return P_FALSE;
1875
1876                 a[0] = angle;
1877                 a[1] = p_vec_angle(v2->co, v1->co, oldv->co);
1878                 a[2] = M_PI - a[0] - a[1];
1879
1880                 b[0] = p_vec_angle(v1->co, keepv->co, v2->co);
1881                 b[1] = p_vec_angle(v2->co, v1->co, keepv->co);
1882                 b[2] = M_PI - b[0] - b[1];
1883
1884                 /* abf criterion 1: avoid sharp and obtuse angles */
1885                 minangle = 15.0f * M_PI / 180.0f;
1886                 maxangle = M_PI - minangle;
1887
1888                 for (i = 0; i < 3; i++) {
1889                         if ((b[i] < a[i]) && (b[i] < minangle))
1890                                 return P_FALSE;
1891                         else if ((b[i] > a[i]) && (b[i] > maxangle))
1892                                 return P_FALSE;
1893                 }
1894
1895                 e = p_wheel_edge_next(e);
1896         } while (e && (e != oldv->edge));
1897
1898         if (p_vert_interior(oldv)) {
1899                 /* hlscm criterion: angular defect smaller than threshold */
1900                 if (fabsf(angulardefect) > (float)(M_PI * 30.0 / 180.0))
1901                         return P_FALSE;
1902         }
1903         else {
1904                 PVert *v1 = p_boundary_edge_next(oldv->edge)->vert;
1905                 PVert *v2 = p_boundary_edge_prev(oldv->edge)->vert;
1906
1907                 /* abf++ criterion 2: avoid collapsing verts inwards */
1908                 if (p_vert_interior(keepv))
1909                         return P_FALSE;
1910
1911                 /* don't collapse significant boundary changes */
1912                 angle = p_vec_angle(v1->co, oldv->co, v2->co);
1913                 if (angle < (M_PI * 160.0 / 180.0))
1914                         return P_FALSE;
1915         }
1916
1917         return P_TRUE;
1918 }
1919
1920 static PBool p_collapse_allowed(PEdge *edge, PEdge *pair)
1921 {
1922         PVert *oldv, *keepv;
1923
1924         p_collapsing_verts(edge, pair, &oldv, &keepv);
1925
1926         if (oldv->flag & PVERT_PIN)
1927                 return P_FALSE;
1928
1929         return (p_collapse_allowed_topologic(edge, pair) &&
1930                 p_collapse_allowed_geometric(edge, pair));
1931 }
1932
1933 static float p_collapse_cost(PEdge *edge, PEdge *pair)
1934 {
1935         /* based on volume and boundary optimization from:
1936          * "Fast and Memory Efficient Polygonal Simplification" P. Lindstrom, G. Turk */
1937
1938         PVert *oldv, *keepv;
1939         PEdge *e;
1940         PFace *oldf1, *oldf2;
1941         float volumecost = 0.0f, areacost = 0.0f, edgevec[3], cost, weight, elen;
1942         float shapecost = 0.0f;
1943         float shapeold = 0.0f, shapenew = 0.0f;
1944         int nshapeold = 0, nshapenew = 0;
1945
1946         p_collapsing_verts(edge, pair, &oldv, &keepv);
1947         oldf1 = (edge) ? edge->face : NULL;
1948         oldf2 = (pair) ? pair->face : NULL;
1949
1950         sub_v3_v3v3(edgevec, keepv->co, oldv->co);
1951
1952         e = oldv->edge;
1953         do {
1954                 float a1, a2, a3;
1955                 float *co1 = e->next->vert->co;
1956                 float *co2 = e->next->next->vert->co;
1957
1958                 if ((e->face != oldf1) && (e->face != oldf2)) {
1959                         float tetrav2[3], tetrav3[3], c[3];
1960
1961                         /* tetrahedron volume = (1/3!)*|a.(b x c)| */
1962                         sub_v3_v3v3(tetrav2, co1, oldv->co);
1963                         sub_v3_v3v3(tetrav3, co2, oldv->co);
1964                         cross_v3_v3v3(c, tetrav2, tetrav3);
1965
1966                         volumecost += fabsf(dot_v3v3(edgevec, c) / 6.0f);
1967 #if 0
1968                         shapecost += dot_v3v3(co1, keepv->co);
1969
1970                         if (p_wheel_edge_next(e) == NULL)
1971                                 shapecost += dot_v3v3(co2, keepv->co);
1972 #endif
1973
1974                         p_triangle_angles(oldv->co, co1, co2, &a1, &a2, &a3);
1975                         a1 = a1 - M_PI / 3.0;
1976                         a2 = a2 - M_PI / 3.0;
1977                         a3 = a3 - M_PI / 3.0;
1978                         shapeold = (a1 * a1 + a2 * a2 + a3 * a3) / ((M_PI / 2) * (M_PI / 2));
1979
1980                         nshapeold++;
1981                 }
1982                 else {
1983                         p_triangle_angles(keepv->co, co1, co2, &a1, &a2, &a3);
1984                         a1 = a1 - M_PI / 3.0;
1985                         a2 = a2 - M_PI / 3.0;
1986                         a3 = a3 - M_PI / 3.0;
1987                         shapenew = (a1 * a1 + a2 * a2 + a3 * a3) / ((M_PI / 2) * (M_PI / 2));
1988
1989                         nshapenew++;
1990                 }
1991
1992                 e = p_wheel_edge_next(e);
1993         } while (e && (e != oldv->edge));
1994
1995         if (!p_vert_interior(oldv)) {
1996                 PVert *v1 = p_boundary_edge_prev(oldv->edge)->vert;
1997                 PVert *v2 = p_boundary_edge_next(oldv->edge)->vert;
1998
1999                 areacost = area_tri_v3(oldv->co, v1->co, v2->co);
2000         }
2001
2002         elen = len_v3(edgevec);
2003         weight = 1.0f; /* 0.2f */
2004         cost = weight * volumecost * volumecost + elen * elen * areacost * areacost;
2005 #if 0
2006         cost += shapecost;
2007 #else
2008         shapeold /= nshapeold;
2009         shapenew /= nshapenew;
2010         shapecost = (shapeold + 0.00001) / (shapenew + 0.00001);
2011
2012         cost *= shapecost;
2013 #endif
2014
2015         return cost;
2016 }
2017
2018 static void p_collapse_cost_vertex(PVert *vert, float *mincost, PEdge **mine)
2019 {
2020         PEdge *e, *enext, *pair;
2021
2022         *mine = NULL;
2023         *mincost = 0.0f;
2024         e = vert->edge;
2025         do {
2026                 if (p_collapse_allowed(e, e->pair)) {
2027                         float cost = p_collapse_cost(e, e->pair);
2028
2029                         if ((*mine == NULL) || (cost < *mincost)) {
2030                                 *mincost = cost;
2031                                 *mine = e;
2032                         }
2033                 }
2034
2035                 enext = p_wheel_edge_next(e);
2036
2037                 if (enext == NULL) {
2038                         /* the other boundary edge, where we only have the pair halfedge */
2039                         pair = e->next->next;
2040
2041                         if (p_collapse_allowed(NULL, pair)) {
2042                                 float cost = p_collapse_cost(NULL, pair);
2043
2044                                 if ((*mine == NULL) || (cost < *mincost)) {
2045                                         *mincost = cost;
2046                                         *mine = pair;
2047                                 }
2048                         }
2049
2050                         break;
2051                 }
2052
2053                 e = enext;
2054         } while (e != vert->edge);
2055 }
2056
2057 static void p_chart_post_collapse_flush(PChart *chart, PEdge *collapsed)
2058 {
2059         /* move to collapsed_ */
2060
2061         PVert *v, *nextv = NULL, *verts = chart->verts;
2062         PEdge *e, *nexte = NULL, *edges = chart->edges, *laste = NULL;
2063         PFace *f, *nextf = NULL, *faces = chart->faces;
2064
2065         chart->verts = chart->collapsed_verts = NULL;
2066         chart->edges = chart->collapsed_edges = NULL;
2067         chart->faces = chart->collapsed_faces = NULL;
2068
2069         chart->nverts = chart->nedges = chart->nfaces = 0;
2070
2071         for (v = verts; v; v = nextv) {
2072                 nextv = v->nextlink;
2073
2074                 if (v->flag & PVERT_COLLAPSE) {
2075                         v->nextlink = chart->collapsed_verts;
2076                         chart->collapsed_verts = v;
2077                 }
2078                 else {
2079                         v->nextlink = chart->verts;
2080                         chart->verts = v;
2081                         chart->nverts++;
2082                 }
2083         }
2084
2085         for (e = edges; e; e = nexte) {
2086                 nexte = e->nextlink;
2087
2088                 if (!collapsed || !(e->flag & PEDGE_COLLAPSE_EDGE)) {
2089                         if (e->flag & PEDGE_COLLAPSE) {
2090                                 e->nextlink = chart->collapsed_edges;
2091                                 chart->collapsed_edges = e;
2092                         }
2093                         else {
2094                                 e->nextlink = chart->edges;
2095                                 chart->edges = e;
2096                                 chart->nedges++;
2097                         }
2098                 }
2099         }
2100
2101         /* these are added last so they can be popped of in the right order
2102          * for splitting */
2103         for (e = collapsed; e; e = e->nextlink) {
2104                 e->nextlink = e->u.nextcollapse;
2105                 laste = e;
2106         }
2107         if (laste) {
2108                 laste->nextlink = chart->collapsed_edges;
2109                 chart->collapsed_edges = collapsed;
2110         }
2111
2112         for (f = faces; f; f = nextf) {
2113                 nextf = f->nextlink;
2114
2115                 if (f->flag & PFACE_COLLAPSE) {
2116                         f->nextlink = chart->collapsed_faces;
2117                         chart->collapsed_faces = f;
2118                 }
2119                 else {
2120                         f->nextlink = chart->faces;
2121                         chart->faces = f;
2122                         chart->nfaces++;
2123                 }
2124         }
2125 }
2126
2127 static void p_chart_post_split_flush(PChart *chart)
2128 {
2129         /* move from collapsed_ */
2130
2131         PVert *v, *nextv = NULL;
2132         PEdge *e, *nexte = NULL;
2133         PFace *f, *nextf = NULL;
2134
2135         for (v = chart->collapsed_verts; v; v = nextv) {
2136                 nextv = v->nextlink;
2137                 v->nextlink = chart->verts;
2138                 chart->verts = v;
2139                 chart->nverts++;
2140         }
2141
2142         for (e = chart->collapsed_edges; e; e = nexte) {
2143                 nexte = e->nextlink;
2144                 e->nextlink = chart->edges;
2145                 chart->edges = e;
2146                 chart->nedges++;
2147         }
2148
2149         for (f = chart->collapsed_faces; f; f = nextf) {
2150                 nextf = f->nextlink;
2151                 f->nextlink = chart->faces;
2152                 chart->faces = f;
2153                 chart->nfaces++;
2154         }
2155
2156         chart->collapsed_verts = NULL;
2157         chart->collapsed_edges = NULL;
2158         chart->collapsed_faces = NULL;
2159 }
2160
2161 static void p_chart_simplify_compute(PChart *chart)
2162 {
2163         /* Computes a list of edge collapses / vertex splits. The collapsed
2164          * simplices go in the chart->collapsed_* lists, The original and
2165          * collapsed may then be view as stacks, where the next collapse/split
2166          * is at the top of the respective lists. */
2167
2168         Heap *heap = BLI_heap_new();
2169         PVert *v, **wheelverts;
2170         PEdge *collapsededges = NULL, *e;
2171         int nwheelverts, i, ncollapsed = 0;
2172
2173         wheelverts = MEM_mallocN(sizeof(PVert *) * chart->nverts, "PChartWheelVerts");
2174
2175         /* insert all potential collapses into heap */
2176         for (v = chart->verts; v; v = v->nextlink) {
2177                 float cost;
2178                 PEdge *e = NULL;
2179
2180                 p_collapse_cost_vertex(v, &cost, &e);
2181
2182                 if (e)
2183                         v->u.heaplink = BLI_heap_insert(heap, cost, e);
2184                 else
2185                         v->u.heaplink = NULL;
2186         }
2187
2188         for (e = chart->edges; e; e = e->nextlink)
2189                 e->u.nextcollapse = NULL;
2190
2191         /* pop edge collapse out of heap one by one */
2192         while (!BLI_heap_is_empty(heap)) {
2193                 if (ncollapsed == NCOLLAPSE)
2194                         break;
2195
2196                 HeapNode *link = BLI_heap_top(heap);
2197                 PEdge *edge = (PEdge *)BLI_heap_pop_min(heap), *pair = edge->pair;
2198                 PVert *oldv, *keepv;
2199                 PEdge *wheele, *nexte;
2200
2201                 /* remember the edges we collapsed */
2202                 edge->u.nextcollapse = collapsededges;
2203                 collapsededges = edge;
2204
2205                 if (edge->vert->u.heaplink != link) {
2206                         edge->flag |= (PEDGE_COLLAPSE_EDGE | PEDGE_COLLAPSE_PAIR);
2207                         edge->next->vert->u.heaplink = NULL;
2208                         SWAP(PEdge *, edge, pair);
2209                 }
2210                 else {
2211                         edge->flag |= PEDGE_COLLAPSE_EDGE;
2212                         edge->vert->u.heaplink = NULL;
2213                 }
2214
2215                 p_collapsing_verts(edge, pair, &oldv, &keepv);
2216
2217                 /* gather all wheel verts and remember them before collapse */
2218                 nwheelverts = 0;
2219                 wheele = oldv->edge;
2220
2221                 do {
2222                         wheelverts[nwheelverts++] = wheele->next->vert;
2223                         nexte = p_wheel_edge_next(wheele);
2224
2225                         if (nexte == NULL)
2226                                 wheelverts[nwheelverts++] = wheele->next->next->vert;
2227
2228                         wheele = nexte;
2229                 } while (wheele && (wheele != oldv->edge));
2230
2231                 /* collapse */
2232                 p_collapse_edge(edge, pair);
2233
2234                 for (i = 0; i < nwheelverts; i++) {
2235                         float cost;
2236                         PEdge *collapse = NULL;
2237
2238                         v = wheelverts[i];
2239
2240                         if (v->u.heaplink) {
2241                                 BLI_heap_remove(heap, v->u.heaplink);
2242                                 v->u.heaplink = NULL;
2243                         }
2244
2245                         p_collapse_cost_vertex(v, &cost, &collapse);
2246
2247                         if (collapse)
2248                                 v->u.heaplink = BLI_heap_insert(heap, cost, collapse);
2249                 }
2250
2251                 ncollapsed++;
2252         }
2253
2254         MEM_freeN(wheelverts);
2255         BLI_heap_free(heap, NULL);
2256
2257         p_chart_post_collapse_flush(chart, collapsededges);
2258 }
2259
2260 static void p_chart_complexify(PChart *chart)
2261 {
2262         PEdge *e, *pair, *edge;
2263         PVert *newv, *keepv;
2264         int x = 0;
2265
2266         for (e = chart->collapsed_edges; e; e = e->nextlink) {
2267                 if (!(e->flag & PEDGE_COLLAPSE_EDGE))
2268                         break;
2269
2270                 edge = e;
2271                 pair = e->pair;
2272
2273                 if (edge->flag & PEDGE_COLLAPSE_PAIR) {
2274                         SWAP(PEdge *, edge, pair);
2275                 }
2276
2277                 p_split_vertex(edge, pair);
2278                 p_collapsing_verts(edge, pair, &newv, &keepv);
2279
2280                 if (x >= NCOLLAPSEX) {
2281                         newv->uv[0] = keepv->uv[0];
2282                         newv->uv[1] = keepv->uv[1];
2283                 }
2284                 else {
2285                         p_vert_harmonic_insert(newv);
2286                         x++;
2287                 }
2288         }
2289
2290         p_chart_post_split_flush(chart);
2291 }
2292
2293 #if 0
2294 static void p_chart_simplify(PChart *chart)
2295 {
2296         /* Not implemented, needs proper reordering in split_flush. */
2297 }
2298 #endif
2299 #endif
2300
2301 /* ABF */
2302
2303 #define ABF_MAX_ITER 20
2304
2305 typedef struct PAbfSystem {
2306         int ninterior, nfaces, nangles;
2307         float *alpha, *beta, *sine, *cosine, *weight;
2308         float *bAlpha, *bTriangle, *bInterior;
2309         float *lambdaTriangle, *lambdaPlanar, *lambdaLength;
2310         float (*J2dt)[3], *bstar, *dstar;
2311         float minangle, maxangle;
2312 } PAbfSystem;
2313
2314 static void p_abf_setup_system(PAbfSystem *sys)
2315 {
2316         int i;
2317
2318         sys->alpha = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFalpha");
2319         sys->beta = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFbeta");
2320         sys->sine = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFsine");
2321         sys->cosine = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFcosine");
2322         sys->weight = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFweight");
2323
2324         sys->bAlpha = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFbalpha");
2325         sys->bTriangle = (float *)MEM_mallocN(sizeof(float) * sys->nfaces, "ABFbtriangle");
2326         sys->bInterior = (float *)MEM_mallocN(sizeof(float) * 2 * sys->ninterior, "ABFbinterior");
2327
2328         sys->lambdaTriangle = (float *)MEM_callocN(sizeof(float) * sys->nfaces, "ABFlambdatri");
2329         sys->lambdaPlanar = (float *)MEM_callocN(sizeof(float) * sys->ninterior, "ABFlamdaplane");
2330         sys->lambdaLength = (float *)MEM_mallocN(sizeof(float) * sys->ninterior, "ABFlambdalen");
2331
2332         sys->J2dt = MEM_mallocN(sizeof(float) * sys->nangles * 3, "ABFj2dt");
2333         sys->bstar = (float *)MEM_mallocN(sizeof(float) * sys->nfaces, "ABFbstar");
2334         sys->dstar = (float *)MEM_mallocN(sizeof(float) * sys->nfaces, "ABFdstar");
2335
2336         for (i = 0; i < sys->ninterior; i++)
2337                 sys->lambdaLength[i] = 1.0;
2338
2339         sys->minangle = 1.0 * M_PI / 180.0;
2340         sys->maxangle = (float)M_PI - sys->minangle;
2341 }
2342
2343 static void p_abf_free_system(PAbfSystem *sys)
2344 {
2345         MEM_freeN(sys->alpha);
2346         MEM_freeN(sys->beta);
2347         MEM_freeN(sys->sine);
2348         MEM_freeN(sys->cosine);
2349         MEM_freeN(sys->weight);
2350         MEM_freeN(sys->bAlpha);
2351         MEM_freeN(sys->bTriangle);
2352         MEM_freeN(sys->bInterior);
2353         MEM_freeN(sys->lambdaTriangle);
2354         MEM_freeN(sys->lambdaPlanar);
2355         MEM_freeN(sys->lambdaLength);
2356         MEM_freeN(sys->J2dt);
2357         MEM_freeN(sys->bstar);
2358         MEM_freeN(sys->dstar);
2359 }
2360
2361 static void p_abf_compute_sines(PAbfSystem *sys)
2362 {
2363         int i;
2364         float *sine = sys->sine, *cosine = sys->cosine, *alpha = sys->alpha;
2365
2366         for (i = 0; i < sys->nangles; i++, sine++, cosine++, alpha++) {
2367                 *sine = sinf(*alpha);
2368                 *cosine = cosf(*alpha);
2369         }
2370 }
2371
2372 static float p_abf_compute_sin_product(PAbfSystem *sys, PVert *v, int aid)
2373 {
2374         PEdge *e, *e1, *e2;
2375         float sin1, sin2;
2376
2377         sin1 = sin2 = 1.0;
2378
2379         e = v->edge;
2380         do {
2381                 e1 = e->next;
2382                 e2 = e->next->next;
2383
2384                 if (aid == e1->u.id) {
2385                         /* we are computing a derivative for this angle,
2386                          * so we use cos and drop the other part */
2387                         sin1 *= sys->cosine[e1->u.id];
2388                         sin2 = 0.0;
2389                 }
2390                 else
2391                         sin1 *= sys->sine[e1->u.id];
2392
2393                 if (aid == e2->u.id) {
2394                         /* see above */
2395                         sin1 = 0.0;
2396                         sin2 *= sys->cosine[e2->u.id];
2397                 }
2398                 else
2399                         sin2 *= sys->sine[e2->u.id];
2400
2401                 e = e->next->next->pair;
2402         } while (e && (e != v->edge));
2403
2404         return (sin1 - sin2);
2405 }
2406
2407 static float p_abf_compute_grad_alpha(PAbfSystem *sys, PFace *f, PEdge *e)
2408 {
2409         PVert *v = e->vert, *v1 = e->next->vert, *v2 = e->next->next->vert;
2410         float deriv;
2411
2412         deriv = (sys->alpha[e->u.id] - sys->beta[e->u.id]) * sys->weight[e->u.id];
2413         deriv += sys->lambdaTriangle[f->u.id];
2414
2415         if (v->flag & PVERT_INTERIOR) {
2416                 deriv += sys->lambdaPlanar[v->u.id];
2417         }
2418
2419         if (v1->flag & PVERT_INTERIOR) {
2420                 float product = p_abf_compute_sin_product(sys, v1, e->u.id);
2421                 deriv += sys->lambdaLength[v1->u.id] * product;
2422         }
2423
2424         if (v2->flag & PVERT_INTERIOR) {
2425                 float product = p_abf_compute_sin_product(sys, v2, e->u.id);
2426                 deriv += sys->lambdaLength[v2->u.id] * product;
2427         }
2428
2429         return deriv;
2430 }
2431
2432 static float p_abf_compute_gradient(PAbfSystem *sys, PChart *chart)
2433 {
2434         PFace *f;
2435         PEdge *e;
2436         PVert *v;
2437         float norm = 0.0;
2438
2439         for (f = chart->faces; f; f = f->nextlink) {
2440                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
2441                 float gtriangle, galpha1, galpha2, galpha3;
2442
2443                 galpha1 = p_abf_compute_grad_alpha(sys, f, e1);
2444                 galpha2 = p_abf_compute_grad_alpha(sys, f, e2);
2445                 galpha3 = p_abf_compute_grad_alpha(sys, f, e3);
2446
2447                 sys->bAlpha[e1->u.id] = -galpha1;
2448                 sys->bAlpha[e2->u.id] = -galpha2;
2449                 sys->bAlpha[e3->u.id] = -galpha3;
2450
2451                 norm += galpha1 * galpha1 + galpha2 * galpha2 + galpha3 * galpha3;
2452
2453                 gtriangle = sys->alpha[e1->u.id] + sys->alpha[e2->u.id] + sys->alpha[e3->u.id] - (float)M_PI;
2454                 sys->bTriangle[f->u.id] = -gtriangle;
2455                 norm += gtriangle * gtriangle;
2456         }
2457
2458         for (v = chart->verts; v; v = v->nextlink) {
2459                 if (v->flag & PVERT_INTERIOR) {
2460                         float gplanar = -2 * M_PI, glength;
2461
2462                         e = v->edge;
2463                         do {
2464                                 gplanar += sys->alpha[e->u.id];
2465                                 e = e->next->next->pair;
2466                         } while (e && (e != v->edge));
2467
2468                         sys->bInterior[v->u.id] = -gplanar;
2469                         norm += gplanar * gplanar;
2470
2471                         glength = p_abf_compute_sin_product(sys, v, -1);
2472                         sys->bInterior[sys->ninterior + v->u.id] = -glength;
2473                         norm += glength * glength;
2474                 }
2475         }
2476
2477         return norm;
2478 }
2479
2480 static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
2481 {
2482         PFace *f;
2483         PEdge *e;
2484         int i, j, ninterior = sys->ninterior, nvar = 2 * sys->ninterior;
2485         PBool success;
2486         LinearSolver *context;
2487
2488         context = EIG_linear_solver_new(0, nvar, 1);
2489
2490         for (i = 0; i < nvar; i++)
2491                 EIG_linear_solver_right_hand_side_add(context, 0, i, sys->bInterior[i]);
2492
2493         for (f = chart->faces; f; f = f->nextlink) {
2494                 float wi1, wi2, wi3, b, si, beta[3], j2[3][3], W[3][3];
2495                 float row1[6], row2[6], row3[6];
2496                 int vid[6];
2497                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
2498                 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
2499
2500                 wi1 = 1.0f / sys->weight[e1->u.id];
2501                 wi2 = 1.0f / sys->weight[e2->u.id];
2502                 wi3 = 1.0f / sys->weight[e3->u.id];
2503
2504                 /* bstar1 = (J1*dInv*bAlpha - bTriangle) */
2505                 b = sys->bAlpha[e1->u.id] * wi1;
2506                 b += sys->bAlpha[e2->u.id] * wi2;
2507                 b += sys->bAlpha[e3->u.id] * wi3;
2508                 b -= sys->bTriangle[f->u.id];
2509
2510                 /* si = J1*d*J1t */
2511                 si = 1.0f / (wi1 + wi2 + wi3);
2512
2513                 /* J1t*si*bstar1 - bAlpha */
2514                 beta[0] = b * si - sys->bAlpha[e1->u.id];
2515                 beta[1] = b * si - sys->bAlpha[e2->u.id];
2516                 beta[2] = b * si - sys->bAlpha[e3->u.id];
2517
2518                 /* use this later for computing other lambda's */
2519                 sys->bstar[f->u.id] = b;
2520                 sys->dstar[f->u.id] = si;
2521
2522                 /* set matrix */
2523                 W[0][0] = si - sys->weight[e1->u.id]; W[0][1] = si; W[0][2] = si;
2524                 W[1][0] = si; W[1][1] = si - sys->weight[e2->u.id]; W[1][2] = si;
2525                 W[2][0] = si; W[2][1] = si; W[2][2] = si - sys->weight[e3->u.id];
2526
2527                 vid[0] = vid[1] = vid[2] = vid[3] = vid[4] = vid[5] = -1;
2528
2529                 if (v1->flag & PVERT_INTERIOR) {
2530                         vid[0] = v1->u.id;
2531                         vid[3] = ninterior + v1->u.id;
2532
2533                         sys->J2dt[e1->u.id][0] = j2[0][0] = 1.0f * wi1;
2534                         sys->J2dt[e2->u.id][0] = j2[1][0] = p_abf_compute_sin_product(sys, v1, e2->u.id) * wi2;
2535                         sys->J2dt[e3->u.id][0] = j2[2][0] = p_abf_compute_sin_product(sys, v1, e3->u.id) * wi3;
2536
2537                         EIG_linear_solver_right_hand_side_add(context, 0, v1->u.id, j2[0][0] * beta[0]);
2538                         EIG_linear_solver_right_hand_side_add(context, 0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
2539
2540                         row1[0] = j2[0][0] * W[0][0];
2541                         row2[0] = j2[0][0] * W[1][0];
2542                         row3[0] = j2[0][0] * W[2][0];
2543
2544                         row1[3] = j2[1][0] * W[0][1] + j2[2][0] * W[0][2];
2545                         row2[3] = j2[1][0] * W[1][1] + j2[2][0] * W[1][2];
2546                         row3[3] = j2[1][0] * W[2][1] + j2[2][0] * W[2][2];
2547                 }
2548
2549                 if (v2->flag & PVERT_INTERIOR) {
2550                         vid[1] = v2->u.id;
2551                         vid[4] = ninterior + v2->u.id;
2552
2553                         sys->J2dt[e1->u.id][1] = j2[0][1] = p_abf_compute_sin_product(sys, v2, e1->u.id) * wi1;
2554                         sys->J2dt[e2->u.id][1] = j2[1][1] = 1.0f * wi2;
2555                         sys->J2dt[e3->u.id][1] = j2[2][1] = p_abf_compute_sin_product(sys, v2, e3->u.id) * wi3;
2556
2557                         EIG_linear_solver_right_hand_side_add(context, 0, v2->u.id, j2[1][1] * beta[1]);
2558                         EIG_linear_solver_right_hand_side_add(context, 0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
2559
2560                         row1[1] = j2[1][1] * W[0][1];
2561                         row2[1] = j2[1][1] * W[1][1];
2562                         row3[1] = j2[1][1] * W[2][1];
2563
2564                         row1[4] = j2[0][1] * W[0][0] + j2[2][1] * W[0][2];
2565                         row2[4] = j2[0][1] * W[1][0] + j2[2][1] * W[1][2];
2566                         row3[4] = j2[0][1] * W[2][0] + j2[2][1] * W[2][2];
2567                 }
2568
2569                 if (v3->flag & PVERT_INTERIOR) {
2570                         vid[2] = v3->u.id;
2571                         vid[5] = ninterior + v3->u.id;
2572
2573                         sys->J2dt[e1->u.id][2] = j2[0][2] = p_abf_compute_sin_product(sys, v3, e1->u.id) * wi1;
2574                         sys->J2dt[e2->u.id][2] = j2[1][2] = p_abf_compute_sin_product(sys, v3, e2->u.id) * wi2;
2575                         sys->J2dt[e3->u.id][2] = j2[2][2] = 1.0f * wi3;
2576
2577                         EIG_linear_solver_right_hand_side_add(context, 0, v3->u.id, j2[2][2] * beta[2]);
2578                         EIG_linear_solver_right_hand_side_add(context, 0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
2579
2580                         row1[2] = j2[2][2] * W[0][2];
2581                         row2[2] = j2[2][2] * W[1][2];
2582                         row3[2] = j2[2][2] * W[2][2];
2583
2584                         row1[5] = j2[0][2] * W[0][0] + j2[1][2] * W[0][1];
2585                         row2[5] = j2[0][2] * W[1][0] + j2[1][2] * W[1][1];
2586                         row3[5] = j2[0][2] * W[2][0] + j2[1][2] * W[2][1];
2587                 }
2588
2589                 for (i = 0; i < 3; i++) {
2590                         int r = vid[i];
2591
2592                         if (r == -1)
2593                                 continue;
2594
2595                         for (j = 0; j < 6; j++) {
2596                                 int c = vid[j];
2597
2598                                 if (c == -1)
2599                                         continue;
2600
2601                                 if (i == 0)
2602                                         EIG_linear_solver_matrix_add(context, r, c, j2[0][i] * row1[j]);
2603                                 else
2604                                         EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[0][i] * row1[j]);
2605
2606                                 if (i == 1)
2607                                         EIG_linear_solver_matrix_add(context, r, c, j2[1][i] * row2[j]);
2608                                 else
2609                                         EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[1][i] * row2[j]);
2610
2611
2612                                 if (i == 2)
2613                                         EIG_linear_solver_matrix_add(context, r, c, j2[2][i] * row3[j]);
2614                                 else
2615                                         EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[2][i] * row3[j]);
2616                         }
2617                 }
2618         }
2619
2620         success = EIG_linear_solver_solve(context);
2621
2622         if (success) {
2623                 for (f = chart->faces; f; f = f->nextlink) {
2624                         float dlambda1, pre[3], dalpha;
2625                         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
2626                         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
2627
2628                         pre[0] = pre[1] = pre[2] = 0.0;
2629
2630                         if (v1->flag & PVERT_INTERIOR) {
2631                                 float x = EIG_linear_solver_variable_get(context, 0, v1->u.id);
2632                                 float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v1->u.id);
2633                                 pre[0] += sys->J2dt[e1->u.id][0] * x;
2634                                 pre[1] += sys->J2dt[e2->u.id][0] * x2;
2635                                 pre[2] += sys->J2dt[e3->u.id][0] * x2;
2636                         }
2637
2638                         if (v2->flag & PVERT_INTERIOR) {
2639                                 float x = EIG_linear_solver_variable_get(context, 0, v2->u.id);
2640                                 float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v2->u.id);
2641                                 pre[0] += sys->J2dt[e1->u.id][1] * x2;
2642                                 pre[1] += sys->J2dt[e2->u.id][1] * x;
2643                                 pre[2] += sys->J2dt[e3->u.id][1] * x2;
2644                         }
2645
2646                         if (v3->flag & PVERT_INTERIOR) {
2647                                 float x = EIG_linear_solver_variable_get(context, 0, v3->u.id);
2648                                 float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v3->u.id);
2649                                 pre[0] += sys->J2dt[e1->u.id][2] * x2;
2650                                 pre[1] += sys->J2dt[e2->u.id][2] * x2;
2651                                 pre[2] += sys->J2dt[e3->u.id][2] * x;
2652                         }
2653
2654                         dlambda1 = pre[0] + pre[1] + pre[2];
2655                         dlambda1 = sys->dstar[f->u.id] * (sys->bstar[f->u.id] - dlambda1);
2656
2657                         sys->lambdaTriangle[f->u.id] += dlambda1;
2658
2659                         dalpha = (sys->bAlpha[e1->u.id] - dlambda1);
2660                         sys->alpha[e1->u.id] += dalpha / sys->weight[e1->u.id] - pre[0];
2661
2662                         dalpha = (sys->bAlpha[e2->u.id] - dlambda1);
2663                         sys->alpha[e2->u.id] += dalpha / sys->weight[e2->u.id] - pre[1];
2664
2665                         dalpha = (sys->bAlpha[e3->u.id] - dlambda1);
2666                         sys->alpha[e3->u.id] += dalpha / sys->weight[e3->u.id] - pre[2];
2667
2668                         /* clamp */
2669                         e = f->edge;
2670                         do {
2671                                 if (sys->alpha[e->u.id] > (float)M_PI)
2672                                         sys->alpha[e->u.id] = (float)M_PI;
2673                                 else if (sys->alpha[e->u.id] < 0.0f)
2674                                         sys->alpha[e->u.id] = 0.0f;
2675                         } while (e != f->edge);
2676                 }
2677
2678                 for (i = 0; i < ninterior; i++) {
2679                         sys->lambdaPlanar[i] += (float)EIG_linear_solver_variable_get(context, 0, i);
2680                         sys->lambdaLength[i] += (float)EIG_linear_solver_variable_get(context, 0, ninterior + i);
2681                 }
2682         }
2683
2684         EIG_linear_solver_delete(context);
2685
2686         return success;
2687 }
2688
2689 static PBool p_chart_abf_solve(PChart *chart)
2690 {
2691         PVert *v;
2692         PFace *f;
2693         PEdge *e, *e1, *e2, *e3;
2694         PAbfSystem sys;
2695         int i;
2696         float /* lastnorm, */ /* UNUSED */ limit = (chart->nfaces > 100) ? 1.0f : 0.001f;
2697
2698         /* setup id's */
2699         sys.ninterior = sys.nfaces = sys.nangles = 0;
2700
2701         for (v = chart->verts; v; v = v->nextlink) {
2702                 if (p_vert_interior(v)) {
2703                         v->flag |= PVERT_INTERIOR;
2704                         v->u.id = sys.ninterior++;
2705                 }
2706                 else
2707                         v->flag &= ~PVERT_INTERIOR;
2708         }
2709
2710         for (f = chart->faces; f; f = f->nextlink) {
2711                 e1 = f->edge; e2 = e1->next; e3 = e2->next;
2712                 f->u.id = sys.nfaces++;
2713
2714                 /* angle id's are conveniently stored in half edges */
2715                 e1->u.id = sys.nangles++;
2716                 e2->u.id = sys.nangles++;
2717                 e3->u.id = sys.nangles++;
2718         }
2719
2720         p_abf_setup_system(&sys);
2721
2722         /* compute initial angles */
2723         for (f = chart->faces; f; f = f->nextlink) {
2724                 float a1, a2, a3;
2725
2726                 e1 = f->edge; e2 = e1->next; e3 = e2->next;
2727                 p_face_angles(f, &a1, &a2, &a3);
2728
2729                 if (a1 < sys.minangle)
2730                         a1 = sys.minangle;
2731                 else if (a1 > sys.maxangle)
2732                         a1 = sys.maxangle;
2733                 if (a2 < sys.minangle)
2734                         a2 = sys.minangle;
2735                 else if (a2 > sys.maxangle)
2736                         a2 = sys.maxangle;
2737                 if (a3 < sys.minangle)
2738                         a3 = sys.minangle;
2739                 else if (a3 > sys.maxangle)
2740                         a3 = sys.maxangle;
2741
2742                 sys.alpha[e1->u.id] = sys.beta[e1->u.id] = a1;
2743                 sys.alpha[e2->u.id] = sys.beta[e2->u.id] = a2;
2744                 sys.alpha[e3->u.id] = sys.beta[e3->u.id] = a3;
2745
2746                 sys.weight[e1->u.id] = 2.0f / (a1 * a1);
2747                 sys.weight[e2->u.id] = 2.0f / (a2 * a2);
2748                 sys.weight[e3->u.id] = 2.0f / (a3 * a3);
2749         }
2750
2751         for (v = chart->verts; v; v = v->nextlink) {
2752                 if (v->flag & PVERT_INTERIOR) {
2753                         float anglesum = 0.0, scale;
2754
2755                         e = v->edge;
2756                         do {
2757                                 anglesum += sys.beta[e->u.id];
2758                                 e = e->next->next->pair;
2759                         } while (e && (e != v->edge));
2760
2761                         scale = (anglesum == 0.0f) ? 0.0f : 2.0f * (float)M_PI / anglesum;
2762
2763                         e = v->edge;
2764                         do {
2765                                 sys.beta[e->u.id] = sys.alpha[e->u.id] = sys.beta[e->u.id] * scale;
2766                                 e = e->next->next->pair;
2767                         } while (e && (e != v->edge));
2768                 }
2769         }
2770
2771         if (sys.ninterior > 0) {
2772                 p_abf_compute_sines(&sys);
2773
2774                 /* iteration */
2775                 /* lastnorm = 1e10; */ /* UNUSED */
2776
2777                 for (i = 0; i < ABF_MAX_ITER; i++) {
2778                         float norm = p_abf_compute_gradient(&sys, chart);
2779
2780                         /* lastnorm = norm; */ /* UNUSED */
2781
2782                         if (norm < limit)
2783                                 break;
2784
2785                         if (!p_abf_matrix_invert(&sys, chart)) {
2786                                 param_warning("ABF failed to invert matrix");
2787                                 p_abf_free_system(&sys);
2788                                 return P_FALSE;
2789                         }
2790
2791                         p_abf_compute_sines(&sys);
2792                 }
2793
2794                 if (i == ABF_MAX_ITER) {
2795                         param_warning("ABF maximum iterations reached");
2796                         p_abf_free_system(&sys);
2797                         return P_FALSE;
2798                 }
2799         }
2800
2801         chart->u.lscm.abf_alpha = MEM_dupallocN(sys.alpha);
2802         p_abf_free_system(&sys);
2803
2804         return P_TRUE;
2805 }
2806
2807 /* Least Squares Conformal Maps */
2808
2809 static void p_chart_pin_positions(PChart *chart, PVert **pin1, PVert **pin2)
2810 {
2811         if (!*pin1 || !*pin2 || *pin1 == *pin2) {
2812                 /* degenerate case */
2813                 PFace *f = chart->faces;
2814                 *pin1 = f->edge->vert;
2815                 *pin2 = f->edge->next->vert;
2816
2817                 (*pin1)->uv[0] = 0.0f;
2818                 (*pin1)->uv[1] = 0.5f;
2819                 (*pin2)->uv[0] = 1.0f;
2820                 (*pin2)->uv[1] = 0.5f;
2821         }
2822         else {
2823                 int diru, dirv, dirx, diry;
2824                 float sub[3];
2825
2826                 sub_v3_v3v3(sub, (*pin1)->co, (*pin2)->co);
2827                 sub[0] = fabsf(sub[0]);
2828                 sub[1] = fabsf(sub[1]);
2829                 sub[2] = fabsf(sub[2]);
2830
2831                 if ((sub[0] > sub[1]) && (sub[0] > sub[2])) {
2832                         dirx = 0;
2833                         diry = (sub[1] > sub[2]) ? 1 : 2;
2834                 }
2835                 else if ((sub[1] > sub[0]) && (sub[1] > sub[2])) {
2836                         dirx = 1;
2837                         diry = (sub[0] > sub[2]) ? 0 : 2;
2838                 }
2839                 else {
2840                         dirx = 2;
2841                         diry = (sub[0] > sub[1]) ? 0 : 1;
2842                 }
2843
2844                 if (dirx == 2) {
2845                         diru = 1;
2846                         dirv = 0;
2847                 }
2848                 else {
2849                         diru = 0;
2850                         dirv = 1;
2851                 }
2852
2853                 (*pin1)->uv[diru] = (*pin1)->co[dirx];
2854                 (*pin1)->uv[dirv] = (*pin1)->co[diry];
2855                 (*pin2)->uv[diru] = (*pin2)->co[dirx];
2856                 (*pin2)->uv[dirv] = (*pin2)->co[diry];
2857         }
2858 }
2859
2860 static PBool p_chart_symmetry_pins(PChart *chart, PEdge *outer, PVert **pin1, PVert **pin2)
2861 {
2862         PEdge *be, *lastbe = NULL, *maxe1 = NULL, *maxe2 = NULL, *be1, *be2;
2863         PEdge *cure = NULL, *firste1 = NULL, *firste2 = NULL, *nextbe;
2864         float maxlen = 0.0f, curlen = 0.0f, totlen = 0.0f, firstlen = 0.0f;
2865         float len1, len2;
2866
2867         /* find longest series of verts split in the chart itself, these are
2868          * marked during construction */
2869         be = outer;
2870         lastbe = p_boundary_edge_prev(be);
2871         do {
2872                 float len = p_edge_length(be);
2873                 totlen += len;
2874
2875                 nextbe = p_boundary_edge_next(be);
2876
2877                 if ((be->vert->flag & PVERT_SPLIT) ||
2878                     (lastbe->vert->flag & nextbe->vert->flag & PVERT_SPLIT))
2879                 {
2880                         if (!cure) {
2881                                 if (be == outer)
2882                                         firste1 = be;
2883                                 cure = be;
2884                         }
2885                         else
2886                                 curlen += p_edge_length(lastbe);
2887                 }
2888                 else if (cure) {
2889                         if (curlen > maxlen) {
2890                                 maxlen = curlen;
2891                                 maxe1 = cure;
2892                                 maxe2 = lastbe;
2893                         }
2894
2895                         if (firste1 == cure) {
2896                                 firstlen = curlen;
2897                                 firste2 = lastbe;
2898                         }
2899
2900                         curlen = 0.0f;
2901                         cure = NULL;
2902                 }
2903
2904                 lastbe = be;
2905                 be = nextbe;
2906         } while (be != outer);
2907
2908         /* make sure we also count a series of splits over the starting point */
2909         if (cure && (cure != outer)) {
2910                 firstlen += curlen + p_edge_length(be);
2911
2912                 if (firstlen > maxlen) {
2913                         maxlen = firstlen;
2914                         maxe1 = cure;
2915                         maxe2 = firste2;
2916                 }
2917         }
2918
2919         if (!maxe1 || !maxe2 || (maxlen < 0.5f * totlen))
2920                 return P_FALSE;
2921
2922         /* find pin1 in the split vertices */
2923         be1 = maxe1;
2924         be2 = maxe2;
2925         len1 = 0.0f;
2926         len2 = 0.0f;
2927
2928         do {
2929                 if (len1 < len2) {
2930                         len1 += p_edge_length(be1);
2931                         be1 = p_boundary_edge_next(be1);
2932                 }
2933                 else {
2934                         be2 = p_boundary_edge_prev(be2);
2935                         len2 += p_edge_length(be2);
2936                 }
2937         } while (be1 != be2);
2938
2939         *pin1 = be1->vert;
2940
2941         /* find pin2 outside the split vertices */
2942         be1 = maxe1;
2943         be2 = maxe2;
2944         len1 = 0.0f;
2945         len2 = 0.0f;
2946
2947         do {
2948                 if (len1 < len2) {
2949                         be1 = p_boundary_edge_prev(be1);
2950                         len1 += p_edge_length(be1);
2951                 }
2952                 else {
2953                         len2 += p_edge_length(be2);
2954                         be2 = p_boundary_edge_next(be2);
2955                 }
2956         } while (be1 != be2);
2957
2958         *pin2 = be1->vert;
2959
2960         p_chart_pin_positions(chart, pin1, pin2);
2961
2962         return !equals_v3v3((*pin1)->co, (*pin2)->co);
2963 }
2964
2965 static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **pin2)
2966 {
2967         float minv[3], maxv[3], dirlen;
2968         PVert *v, *minvert[3], *maxvert[3];
2969         int i, dir;
2970
2971         /* find minimum and maximum verts over x/y/z axes */
2972         minv[0] = minv[1] = minv[2] = 1e20;
2973         maxv[0] = maxv[1] = maxv[2] = -1e20;
2974
2975         minvert[0] = minvert[1] = minvert[2] = NULL;
2976         maxvert[0] = maxvert[1] = maxvert[2] = NULL;
2977
2978         for (v = chart->verts; v; v = v->nextlink) {
2979                 for (i = 0; i < 3; i++) {
2980                         if (v->co[i] < minv[i]) {
2981                                 minv[i] = v->co[i];
2982                                 minvert[i] = v;
2983                         }
2984                         if (v->co[i] > maxv[i]) {
2985                                 maxv[i] = v->co[i];
2986                                 maxvert[i] = v;
2987                         }
2988                 }
2989         }
2990
2991         /* find axes with longest distance */
2992         dir = 0;
2993         dirlen = -1.0;
2994
2995         for (i = 0; i < 3; i++) {
2996                 if (maxv[i] - minv[i] > dirlen) {
2997                         dir = i;
2998                         dirlen = maxv[i] - minv[i];
2999                 }
3000         }
3001
3002         *pin1 = minvert[dir];
3003         *pin2 = maxvert[dir];
3004
3005         p_chart_pin_positions(chart, pin1, pin2);
3006 }
3007
3008 static void p_chart_lscm_load_solution(PChart *chart)
3009 {
3010         LinearSolver *context = chart->u.lscm.context;
3011         PVert *v;
3012
3013         for (v = chart->verts; v; v = v->nextlink) {
3014                 v->uv[0] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id);
3015                 v->uv[1] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id + 1);
3016         }
3017 }
3018
3019 static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
3020 {
3021         PVert *v, *pin1, *pin2;
3022         PBool select = P_FALSE, deselect = P_FALSE;
3023         int npins = 0, id = 0;
3024
3025         /* give vertices matrix indices and count pins */
3026         for (v = chart->verts; v; v = v->nextlink) {
3027                 if (v->flag & PVERT_PIN) {
3028                         npins++;
3029                         if (v->flag & PVERT_SELECT)
3030                                 select = P_TRUE;
3031                 }
3032
3033                 if (!(v->flag & PVERT_SELECT))
3034                         deselect = P_TRUE;
3035         }
3036
3037         if ((live && (!select || !deselect)) || (npins == 1)) {
3038                 chart->u.lscm.context = NULL;
3039         }
3040         else {
3041 #if 0
3042                 p_chart_simplify_compute(chart);
3043                 p_chart_topological_sanity_check(chart);
3044 #endif
3045
3046                 if (abf) {
3047                         if (!p_chart_abf_solve(chart))
3048                                 param_warning("ABF solving failed: falling back to LSCM.\n");
3049                 }
3050
3051                 if (npins <= 1) {
3052                         /* not enough pins, lets find some ourself */
3053                         PEdge *outer;
3054
3055                         p_chart_boundaries(chart, NULL, &outer);
3056
3057                         /* outer can be NULL with non-finite coords. */
3058                         if (!(outer && p_chart_symmetry_pins(chart, outer, &pin1, &pin2))) {
3059                                 p_chart_extrema_verts(chart, &pin1, &pin2);
3060                         }
3061
3062                         chart->u.lscm.pin1 = pin1;
3063                         chart->u.lscm.pin2 = pin2;
3064                 }
3065                 else {
3066                         chart->flag |= PCHART_NOPACK;
3067                 }
3068
3069                 for (v = chart->verts; v; v = v->nextlink)
3070                         v->u.id = id++;
3071
3072                 chart->u.lscm.context = EIG_linear_least_squares_solver_new(2 * chart->nfaces, 2 * chart->nverts, 1);
3073         }
3074 }
3075
3076 static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
3077 {
3078         LinearSolver *context = chart->u.lscm.context;
3079         PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
3080         PFace *f;
3081         float *alpha = chart->u.lscm.abf_alpha;
3082         float area_pinned_up, area_pinned_down;
3083         bool flip_faces;
3084         int row;
3085
3086 #if 0
3087         /* TODO: make loading pins work for simplify/complexify. */
3088 #endif
3089
3090         for (v = chart->verts; v; v = v->nextlink)
3091                 if (v->flag & PVERT_PIN)
3092                         p_vert_load_pin_select_uvs(handle, v);  /* reload for live */
3093
3094         if (chart->u.lscm.pin1) {
3095                 EIG_linear_solver_variable_lock(context, 2 * pin1->u.id);
3096                 EIG_linear_solver_variable_lock(context, 2 * pin1->u.id + 1);
3097                 EIG_linear_solver_variable_lock(context, 2 * pin2->u.id);
3098                 EIG_linear_solver_variable_lock(context, 2 * pin2->u.id + 1);
3099
3100                 EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id, pin1->uv[0]);
3101                 EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id + 1, pin1->uv[1]);
3102                 EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id, pin2->uv[0]);
3103                 EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id + 1, pin2->uv[1]);
3104         }
3105         else {
3106                 /* set and lock the pins */
3107                 for (v = chart->verts; v; v = v->nextlink) {
3108                         if (v->flag & PVERT_PIN) {
3109                                 EIG_linear_solver_variable_lock(context, 2 * v->u.id);
3110                                 EIG_linear_solver_variable_lock(context, 2 * v->u.id + 1);
3111
3112                                 EIG_linear_solver_variable_set(context, 0, 2 * v->u.id, v->uv[0]);
3113                                 EIG_linear_solver_variable_set(context, 0, 2 * v->u.id + 1, v->uv[1]);
3114                         }
3115                 }
3116         }
3117
3118         /* detect up direction based on pinned vertices */
3119         area_pinned_up = 0.0f;
3120         area_pinned_down = 0.0f;
3121
3122         for (f = chart->faces; f; f = f->nextlink) {
3123                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
3124                 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
3125
3126                 if ((v1->flag & PVERT_PIN) && (v2->flag & PVERT_PIN) && (v3->flag & PVERT_PIN)) {
3127                         float area = p_face_uv_area_signed(f);
3128
3129                         if (area > 0.0f)
3130                                 area_pinned_up += area;
3131                         else
3132                                 area_pinned_down -= area;
3133                 }
3134         }
3135
3136         flip_faces = (area_pinned_down > area_pinned_up);
3137
3138         /* construct matrix */
3139
3140         row = 0;
3141         for (f = chart->faces; f; f = f->nextlink) {
3142                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
3143                 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
3144                 float a1, a2, a3, ratio, cosine, sine;
3145                 float sina1, sina2, sina3, sinmax;
3146
3147                 if (alpha) {
3148                         /* use abf angles if passed on */
3149                         a1 = *(alpha++);
3150                         a2 = *(alpha++);
3151                         a3 = *(alpha++);
3152                 }
3153                 else
3154                         p_face_angles(f, &a1, &a2, &a3);
3155
3156                 if (flip_faces) {
3157                         SWAP(float, a2, a3);
3158                         SWAP(PEdge *, e2, e3);
3159                         SWAP(PVert *, v2, v3);
3160                 }
3161
3162                 sina1 = sinf(a1);
3163                 sina2 = sinf(a2);
3164                 sina3 = sinf(a3);
3165
3166                 sinmax = max_fff(sina1, sina2, sina3);
3167
3168                 /* shift vertices to find most stable order */
3169                 if (sina3 != sinmax) {
3170                         SHIFT3(PVert *, v1, v2, v3);
3171                         SHIFT3(float, a1, a2, a3);
3172                         SHIFT3(float, sina1, sina2, sina3);
3173
3174                         if (sina2 == sinmax) {
3175                                 SHIFT3(PVert *, v1, v2, v3);
3176                                 SHIFT3(float, a1, a2, a3);
3177                                 SHIFT3(float, sina1, sina2, sina3);
3178                         }
3179                 }
3180
3181                 /* angle based lscm formulation */
3182                 ratio = (sina3 == 0.0f) ? 1.0f : sina2 / sina3;
3183                 cosine = cosf(a1) * ratio;
3184                 sine = sina1 * ratio;
3185
3186                 EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id,   cosine - 1.0f);
3187                 EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, -sine);
3188                 EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id,   -cosine);
3189                 EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, sine);
3190                 EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id,   1.0);
3191                 row++;
3192
3193                 EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id,   sine);
3194                 EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, cosine - 1.0f);
3195                 EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id,   -sine);
3196                 EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, -cosine);
3197                 EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id + 1, 1.0);
3198                 row++;
3199         }
3200
3201         if (EIG_linear_solver_solve(context)) {
3202                 p_chart_lscm_load_solution(chart);
3203                 return P_TRUE;
3204         }
3205         else {
3206                 for (v = chart->verts; v; v = v->nextlink) {
3207                         v->uv[0] = 0.0f;
3208                         v->uv[1] = 0.0f;
3209                 }
3210         }
3211
3212         return P_FALSE;
3213 }
3214
3215 static void p_chart_lscm_end(PChart *chart)
3216 {
3217         if (chart->u.lscm.context)
3218                 EIG_linear_solver_delete(chart->u.lscm.context);
3219
3220         if (chart->u.lscm.abf_alpha) {
3221                 MEM_freeN(chart->u.lscm.abf_alpha);
3222                 chart->u.lscm.abf_alpha = NULL;
3223         }
3224