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