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