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