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