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