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