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