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