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