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