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