Wednesday merger of Orange branch and bf-blender
[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 /* Hash */
28
29 static int PHashSizes[] = {
30         1, 3, 5, 11, 17, 37, 67, 131, 257, 521, 1031, 2053, 4099, 8209, 
31         16411, 32771, 65537, 131101, 262147, 524309, 1048583, 2097169, 
32         4194319, 8388617, 16777259, 33554467, 67108879, 134217757, 268435459
33 };
34
35 #define PHASH_hash(ph, item) (((unsigned long) (item))%((unsigned int) (ph)->cursize))
36
37 PHash *phash_new(int sizehint)
38 {
39         PHash *ph = (PHash*)MEM_callocN(sizeof(PHash), "PHash");
40         ph->size = 0;
41         ph->cursize_id = 0;
42         ph->first = NULL;
43
44         while (PHashSizes[ph->cursize_id] < sizehint)
45                 ph->cursize_id++;
46
47         ph->cursize = PHashSizes[ph->cursize_id];
48         ph->buckets = (PHashLink**)MEM_callocN(ph->cursize*sizeof(*ph->buckets), "PHashBuckets");
49
50         return ph;
51 }
52
53 void phash_delete(PHash *ph)
54 {
55         MEM_freeN(ph->buckets);
56         MEM_freeN(ph);
57 }
58
59 void phash_delete_with_links(PHash *ph)
60 {
61         PHashLink *link, *next=NULL;
62
63         for (link = ph->first; link; link = next) {
64                 next = link->next;
65                 MEM_freeN(link);
66         }
67
68         phash_delete(ph);
69 }
70
71 int phash_size(PHash *ph)
72 {
73         return ph->size;
74 }
75
76 void phash_insert(PHash *ph, PHashLink *link)
77 {
78         int size = ph->cursize;
79         int hash = PHASH_hash(ph, link->key);
80         PHashLink *lookup = ph->buckets[hash];
81
82         if (lookup == NULL) {
83                 /* insert in front of the list */
84                 ph->buckets[hash] = link;
85                 link->next = ph->first;
86                 ph->first = link;
87         }
88         else {
89                 /* insert after existing element */
90                 link->next = lookup->next;
91                 lookup->next = link;
92         }
93                 
94         ph->size++;
95
96         if (ph->size > (size*3)) {
97                 PHashLink *next = NULL, *first = ph->first;
98
99                 ph->cursize = PHashSizes[++ph->cursize_id];
100                 MEM_freeN(ph->buckets);
101                 ph->buckets = (PHashLink**)MEM_callocN(ph->cursize*sizeof(*ph->buckets), "PHashBuckets");
102                 ph->size = 0;
103                 ph->first = NULL;
104
105                 for (link = first; link; link = next) {
106                         next = link->next;
107                         phash_insert(ph, link);
108                 }
109         }
110 }
111
112 PHashLink *phash_lookup(PHash *ph, PHashKey key)
113 {
114         PHashLink *link;
115         int hash = PHASH_hash(ph, key);
116
117         for (link = ph->buckets[hash]; link; link = link->next)
118                 if (link->key == key)
119                         return link;
120                 else if (PHASH_hash(ph, link->key) != hash)
121                         return NULL;
122         
123         return link;
124 }
125
126 PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
127 {
128         int hash = PHASH_hash(ph, key);
129
130         for (link = link->next; link; link = link->next)
131                 if (link->key == key)
132                         return link;
133                 else if (PHASH_hash(ph, link->key) != hash)
134                         return NULL;
135         
136         return link;
137 }
138
139 /* Heap */
140
141 #define PHEAP_PARENT(i) ((i-1)>>1)
142 #define PHEAP_LEFT(i)   ((i<<1)+1)
143 #define PHEAP_RIGHT(i)  ((i<<1)+2)
144 #define PHEAP_COMPARE(a, b) (a->value < b->value)
145 #define PHEAP_EQUALS(a, b) (a->value == b->value)
146 #define PHEAP_SWAP(heap, i, j) \
147         { SWAP(int, heap->tree[i]->index, heap->tree[j]->index); \
148           SWAP(PHeapLink*, heap->tree[i], heap->tree[j]);  }
149
150 static void pheap_down(PHeap *heap, int i)
151 {
152         while (P_TRUE) {
153                 int size = heap->size, smallest;
154                 int l = PHEAP_LEFT(i);
155                 int r = PHEAP_RIGHT(i);
156
157                 smallest = ((l < size) && PHEAP_COMPARE(heap->tree[l], heap->tree[i]))? l: i;
158
159                 if ((r < size) && PHEAP_COMPARE(heap->tree[r], heap->tree[smallest]))
160                         smallest = r;
161                 
162                 if (smallest == i)
163                         break;
164
165                 PHEAP_SWAP(heap, i, smallest);
166                 i = smallest;
167         }
168 }
169
170 static void pheap_up(PHeap *heap, int i)
171 {
172         while (i > 0) {
173                 int p = PHEAP_PARENT(i);
174
175                 if (PHEAP_COMPARE(heap->tree[p], heap->tree[i]))
176                         break;
177
178                 PHEAP_SWAP(heap, p, i);
179                 i = p;
180         }
181 }
182
183 PHeap *pheap_new()
184 {
185         /* TODO: replace mallocN with something faster */
186
187         PHeap *heap = (PHeap*)MEM_callocN(sizeof(PHeap), "PHeap");
188         heap->bufsize = 1;
189         heap->tree = (PHeapLink**)MEM_mallocN(sizeof(PHeapLink*), "PHeapTree");
190
191         return heap;
192 }
193
194 void pheap_delete(PHeap *heap)
195 {
196         MEM_freeN(heap->tree);
197         MEM_freeN(heap);
198 }
199
200 PHeapLink *pheap_insert(PHeap *heap, float value, void *ptr)
201 {
202         PHeapLink *link;
203
204         if ((heap->size + 1) > heap->bufsize) {
205                 int newsize = heap->bufsize*2;
206
207                 PHeapLink **ntree = (PHeapLink**)MEM_mallocN(newsize*sizeof(PHeapLink*), "PHeapTree");
208                 memcpy(ntree, heap->tree, sizeof(PHeapLink*)*heap->size);
209                 MEM_freeN(heap->tree);
210
211                 heap->tree = ntree;
212                 heap->bufsize = newsize;
213         }
214
215         param_assert(heap->size < heap->bufsize);
216
217         link = MEM_mallocN(sizeof *link, "PHeapLink");
218         link->value = value;
219         link->ptr = ptr;
220         link->index = heap->size;
221
222         heap->tree[link->index] = link;
223
224         heap->size++;
225
226         pheap_up(heap, heap->size-1);
227
228         return link;
229 }
230
231 int pheap_empty(PHeap *heap)
232 {
233         return (heap->size == 0);
234 }
235
236 int pheap_size(PHeap *heap)
237 {
238         return heap->size;
239 }
240
241 void *pheap_min(PHeap *heap)
242 {
243         return heap->tree[0]->ptr;
244 }
245
246 void *pheap_popmin(PHeap *heap)
247 {
248         void *ptr = heap->tree[0]->ptr;
249
250         MEM_freeN(heap->tree[0]);
251
252         if (heap->size == 1)
253                 heap->size--;
254         else {
255                 PHEAP_SWAP(heap, 0, heap->size-1);
256                 heap->size--;
257
258                 pheap_down(heap, 0);
259         }
260
261         return ptr;
262 }
263
264 static void pheap_remove(PHeap *heap, PHeapLink *link)
265 {
266         int i = link->index;
267
268         while (i > 0) {
269                 int p = PHEAP_PARENT(i);
270
271                 PHEAP_SWAP(heap, p, i);
272                 i = p;
273         }
274
275         pheap_popmin(heap);
276 }
277
278 /* Construction */
279
280 PEdge *p_wheel_edge_next(PEdge *e)
281 {
282         return e->next->next->pair;
283 }
284
285 PEdge *p_wheel_edge_prev(PEdge *e)
286 {   
287         return (e->pair)? e->pair->next: NULL;
288 }
289
290 static PVert *p_vert_add(PChart *chart, PHashKey key, float *co, PEdge *e)
291 {
292         PVert *v = (PVert*)BLI_memarena_alloc(chart->handle->arena, sizeof *v);
293         v->co = co;
294         v->link.key = key;
295         v->edge = e;
296         v->flag = 0;
297
298         phash_insert(chart->verts, (PHashLink*)v);
299
300         return v;
301 }
302
303 static PVert *p_vert_lookup(PChart *chart, PHashKey key, float *co, PEdge *e)
304 {
305         PVert *v = (PVert*)phash_lookup(chart->verts, key);
306
307         if (v)
308                 return v;
309         else
310                 return p_vert_add(chart, key, co, e);
311 }
312
313 static PVert *p_vert_copy(PChart *chart, PVert *v)
314 {
315         PVert *nv = (PVert*)BLI_memarena_alloc(chart->handle->arena, sizeof *nv);
316         nv->co = v->co;
317         nv->uv[0] = v->uv[0];
318         nv->uv[1] = v->uv[1];
319         nv->link.key = v->link.key;
320         nv->edge = v->edge;
321         nv->flag = v->flag;
322
323         phash_insert(chart->verts, (PHashLink*)nv);
324
325         return nv;
326 }
327
328 static PEdge *p_edge_lookup(PChart *chart, PHashKey *vkeys)
329 {
330         PHashKey key = vkeys[0]^vkeys[1];
331         PEdge *e = (PEdge*)phash_lookup(chart->edges, key);
332
333         while (e) {
334                 if ((e->vert->link.key == vkeys[0]) && (e->next->vert->link.key == vkeys[1]))
335                         return e;
336                 else if ((e->vert->link.key == vkeys[1]) && (e->next->vert->link.key == vkeys[0]))
337                         return e;
338
339                 e = (PEdge*)phash_next(chart->edges, key, (PHashLink*)e);
340         }
341
342         return NULL;
343 }
344
345 static void p_face_flip(PFace *f)
346 {
347         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
348         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
349         int f1 = e1->flag, f2 = e2->flag, f3 = e3->flag;
350
351         e1->vert = v2;
352         e1->next = e3;
353         e1->flag = (f1 & ~PEDGE_VERTEX_FLAGS) | (f2 & PEDGE_VERTEX_FLAGS);
354
355         e2->vert = v3;
356         e2->next = e1;
357         e2->flag = (f2 & ~PEDGE_VERTEX_FLAGS) | (f3 & PEDGE_VERTEX_FLAGS);
358
359         e3->vert = v1;
360         e3->next = e2;
361         e3->flag = (f3 & ~PEDGE_VERTEX_FLAGS) | (f1 & PEDGE_VERTEX_FLAGS);
362 }
363
364 static void p_vert_load_pin_select_uvs(PVert *v)
365 {
366         PEdge *e;
367         int nedges = 0;
368
369         v->uv[0] = v->uv[1] = 0.0f;
370         nedges = 0;
371         e = v->edge;
372         do {
373                 if (e->orig_uv && (e->flag & PEDGE_PIN)) {
374                         if (e->flag & PEDGE_SELECT)
375                                 v->flag |= PVERT_SELECT;
376
377                         v->flag |= PVERT_PIN;
378                         v->uv[0] += e->orig_uv[0];
379                         v->uv[1] += e->orig_uv[1];
380                         nedges++;
381                 }
382
383                 e = p_wheel_edge_next(e);
384         } while (e && e != (v->edge));
385
386         if (nedges > 0) {
387                 v->uv[0] /= nedges;
388                 v->uv[1] /= nedges;
389         }
390 }
391
392 static void p_vert_load_select_uvs(PVert *v)
393 {
394         PEdge *e;
395         int nedges = 0;
396
397         v->uv[0] = v->uv[1] = 0.0f;
398         nedges = 0;
399         e = v->edge;
400         do {
401                 if (e->orig_uv && (e->flag & PEDGE_SELECT))
402                         v->flag |= PVERT_SELECT;
403
404                 v->uv[0] += e->orig_uv[0];
405                 v->uv[1] += e->orig_uv[1];
406                 nedges++;
407
408                 e = p_wheel_edge_next(e);
409         } while (e && e != (v->edge));
410
411         if (nedges > 0) {
412                 v->uv[0] /= nedges;
413                 v->uv[1] /= nedges;
414         }
415 }
416
417 static void p_extrema_verts(PChart *chart, PVert **v1, PVert **v2)
418 {
419         float minv[3], maxv[3], dirlen;
420         PVert *v, *minvert[3], *maxvert[3];
421         int i, dir;
422
423         /* find minimum and maximum verts over x/y/z axes */
424         minv[0] = minv[1] = minv[2] = 1e20;
425         maxv[0] = maxv[1] = maxv[2] = -1e20;
426
427         minvert[0] = minvert[1] = minvert[2] = NULL;
428         maxvert[0] = maxvert[1] = maxvert[2] = NULL;
429
430         for (v = (PVert*)chart->verts->first; v; v=v->link.next) {
431                 for (i = 0; i < 3; i++) {
432                         if (v->co[i] < minv[i]) {
433                                 minv[i] = v->co[i];
434                                 minvert[i] = v;
435                         }
436                         if (v->co[i] > maxv[i]) {
437                                 maxv[i] = v->co[i];
438                                 maxvert[i] = v;
439                         }
440                 }
441         }
442
443         /* find axes with longest distance */
444         dir = 0;
445         dirlen = -1.0;
446
447         for (i = 0; i < 3; i++) {
448                 if (maxv[i] - minv[i] > dirlen) {
449                         dir = i;
450                         dirlen = maxv[i] - minv[i];
451                 }
452         }
453
454         if (minvert[dir] == maxvert[dir]) {
455                 /* degenerate case */
456                 PFace *f = (PFace*)chart->faces->first;
457                 *v1 = f->edge->vert;
458                 *v2 = f->edge->next->vert;
459
460                 (*v1)->uv[0] = 0.0f;
461                 (*v1)->uv[1] = 0.5f;
462                 (*v2)->uv[0] = 1.0f;
463                 (*v2)->uv[1] = 0.5f;
464         }
465         else {
466                 *v1 = minvert[dir];
467                 *v2 = maxvert[dir];
468
469                 (*v1)->uv[0] = (*v1)->co[dir];
470                 (*v1)->uv[1] = (*v1)->co[(dir+1)%3];
471                 (*v2)->uv[0] = (*v2)->co[dir];
472                 (*v2)->uv[1] = (*v2)->co[(dir+1)%3];
473         }
474 }
475
476 static float p_vec_normalise(float *v)
477 {
478    float d;
479    
480     d = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
481
482     if(d != 0.0f) {
483                 d = 1.0f/d;
484
485                 v[0] *= d;
486                 v[1] *= d;
487                 v[2] *= d;
488         }
489
490         return d;
491 }
492
493 static float p_vec_angle_cos(float *v1, float *v2, float *v3)
494 {
495         float d1[3], d2[3];
496
497         d1[0] = v1[0] - v2[0];
498         d1[1] = v1[1] - v2[1];
499         d1[2] = v1[2] - v2[2];
500
501         d2[0] = v3[0] - v2[0];
502         d2[1] = v3[1] - v2[1];
503         d2[2] = v3[2] - v2[2];
504
505         p_vec_normalise(d1);
506         p_vec_normalise(d2);
507
508         return d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2];
509 }
510
511 static float p_vec_angle(float *v1, float *v2, float *v3)
512 {
513         float dot = p_vec_angle_cos(v1, v2, v3);
514
515         if (dot <= -1.0f)
516                 return (float)M_PI;
517         else if (dot >= 1.0f)
518                 return 0.0f;
519         else
520                 return (float)acos(dot);
521 }
522
523 static void p_face_angles(PFace *f, float *a1, float *a2, float *a3)
524 {
525         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
526         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
527
528         *a1 = p_vec_angle(v3->co, v1->co, v2->co);
529         *a2 = p_vec_angle(v1->co, v2->co, v3->co);
530         *a3 = M_PI - *a2 - *a1;
531 }
532
533 static float p_face_area(PFace *f)
534 {
535         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
536         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
537
538         return AreaT3Dfl(v1->co, v2->co, v3->co);
539 }
540
541 static float p_face_uv_area_signed(PFace *f)
542 {
543         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
544         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
545
546         return 0.5f*(((v2->uv[0]-v1->uv[0]) * (v3->uv[1]-v1->uv[1])) - 
547                     ((v3->uv[0]-v1->uv[0]) * (v2->uv[1]-v1->uv[1])));
548 }
549
550 static float p_face_uv_area(PFace *f)
551 {
552         return fabs(p_face_uv_area_signed(f));
553 }
554
555 static void p_chart_area(PChart *chart, float *uv_area, float *area)
556 {
557         PFace *f;
558
559         *uv_area = *area = 0.0f;
560
561         for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
562                 *uv_area += p_face_uv_area(f);
563                 *area += p_face_area(f);
564         }
565 }
566
567 static PChart *p_chart_new(PHandle *handle)
568 {
569         PChart *chart = (PChart*)MEM_callocN(sizeof*chart, "PChart");
570         chart->verts = phash_new(1);
571         chart->edges = phash_new(1);
572         chart->faces = phash_new(1);
573         chart->handle = handle;
574
575         return chart;
576 }
577
578 static void p_chart_delete(PChart *chart)
579 {
580         /* the actual links are free by memarena */
581         phash_delete(chart->verts);
582         phash_delete(chart->edges);
583         phash_delete(chart->faces);
584
585         MEM_freeN(chart);
586 }
587
588 static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
589 {
590         float *uv1, *uv2, *uvp1, *uvp2;
591         float limit[2];
592
593         uv1 = e->orig_uv;
594         uv2 = e->next->orig_uv;
595
596         if (e->vert->link.key == ep->vert->link.key) {
597                 uvp1 = ep->orig_uv;
598                 uvp2 = ep->next->orig_uv;
599         }
600         else {
601                 uvp1 = ep->next->orig_uv;
602                 uvp2 = ep->orig_uv;
603         }
604
605         get_connected_limit_tface_uv(limit);
606
607         if((fabs(uv1[0]-uvp1[0]) > limit[0]) && (fabs(uv1[1]-uvp1[1]) > limit[1])) {
608                 e->flag |= PEDGE_SEAM;
609                 ep->flag |= PEDGE_SEAM;
610                 return P_TRUE;
611         }
612         if((fabs(uv2[0]-uvp2[0]) > limit[0]) && (fabs(uv2[1]-uvp2[1]) > limit[1])) {
613                 e->flag |= PEDGE_SEAM;
614                 ep->flag |= PEDGE_SEAM;
615                 return P_TRUE;
616         }
617         
618         return P_FALSE;
619 }
620
621 static PBool p_edge_has_pair(PChart *chart, PEdge *e, PEdge **pair, PBool impl)
622 {
623         PHashKey key;
624         PEdge *pe;
625         PVert *v1, *v2;
626         PHashKey key1 = e->vert->link.key;
627         PHashKey key2 = e->next->vert->link.key;
628
629         if (e->flag & PEDGE_SEAM)
630                 return P_FALSE;
631         
632         key = key1 ^ key2;
633         pe = (PEdge*)phash_lookup(chart->edges, key);
634         *pair = NULL;
635
636         while (pe) {
637                 if (pe != e) {
638                         v1 = pe->vert;
639                         v2 = pe->next->vert;
640
641                         if (((v1->link.key == key1) && (v2->link.key == key2)) ||
642                                 ((v1->link.key == key2) && (v2->link.key == key1))) {
643
644                                 /* don't connect seams and t-junctions */
645                                 if ((pe->flag & PEDGE_SEAM) || *pair ||
646                                     (impl && p_edge_implicit_seam(e, pe))) {
647                                         *pair = NULL;
648                                         return P_FALSE;
649                                 }
650
651                                 *pair = pe;
652                         }
653                 }
654
655                 pe = (PEdge*)phash_next(chart->edges, key, (PHashLink*)pe);
656         }
657
658         if (*pair && (e->vert == (*pair)->vert)) {
659                 if ((*pair)->next->pair || (*pair)->next->next->pair) {
660                         /* non unfoldable, maybe mobius ring or klein bottle */
661                         *pair = NULL;
662                         return P_FALSE;
663                 }
664         }
665
666         return (*pair != NULL);
667 }
668
669 static PBool p_edge_connect_pair(PChart *chart, PEdge *e, PEdge ***stack, PBool impl)
670 {
671         PEdge *pair = NULL;
672
673         if(!e->pair && p_edge_has_pair(chart, e, &pair, impl)) {
674                 if (e->vert == pair->vert)
675                         p_face_flip(pair->face);
676
677                 e->pair = pair;
678                 pair->pair = e;
679
680                 if (!(pair->face->flag & PFACE_CONNECTED)) {
681                         **stack = pair;
682                         (*stack)++;
683                 }
684         }
685
686         return (e->pair != NULL);
687 }
688
689 static int p_connect_pairs(PChart *chart, PBool impl)
690 {
691         PEdge **stackbase = MEM_mallocN(sizeof*stackbase * phash_size(chart->faces), "Pstackbase");
692         PEdge **stack = stackbase;
693         PFace *f, *first;
694         PEdge *e, *e1, *e2;
695         int ncharts = 0;
696
697         /* connect pairs, count edges, set vertex-edge pointer to a pairless edge */
698         for (first=(PFace*)chart->faces->first; first; first=first->link.next) {
699                 if (first->flag & PFACE_CONNECTED)
700                         continue;
701
702                 *stack = first->edge;
703                 stack++;
704
705                 while (stack != stackbase) {
706                         stack--;
707                         e = *stack;
708                         e1 = e->next;
709                         e2 = e1->next;
710
711                         f = e->face;
712                         f->flag |= PFACE_CONNECTED;
713
714                         /* assign verts to charts so we can sort them later */
715                         f->u.chart = ncharts;
716
717                         if (!p_edge_connect_pair(chart, e, &stack, impl))
718                                 e->vert->edge = e;
719                         if (!p_edge_connect_pair(chart, e1, &stack, impl))
720                                 e1->vert->edge = e1;
721                         if (!p_edge_connect_pair(chart, e2, &stack, impl))
722                                 e2->vert->edge = e2;
723                 }
724
725                 ncharts++;
726         }
727
728         MEM_freeN(stackbase);
729
730         return ncharts;
731 }
732
733 static void p_split_vert(PChart *chart, PEdge *e)
734 {
735         PEdge *we, *lastwe = NULL;
736         PVert *v = e->vert;
737         PBool copy = P_TRUE;
738
739         if (e->flag & PEDGE_VERTEX_SPLIT)
740                 return;
741
742         /* rewind to start */
743         lastwe = e;
744         for (we = p_wheel_edge_prev(e); we && (we != e); we = p_wheel_edge_prev(we))
745                 lastwe = we;
746         
747         /* go over all edges in wheel */
748         for (we = lastwe; we; we = p_wheel_edge_next(we)) {
749                 if (we->flag & PEDGE_VERTEX_SPLIT)
750                         break;
751
752                 we->flag |= PEDGE_VERTEX_SPLIT;
753
754                 if (we == v->edge) {
755                         /* found it, no need to copy */
756                         copy = P_FALSE;
757                         phash_insert(chart->verts, (PHashLink*)v);
758                 }
759         }
760
761         if (copy) {
762                 /* not found, copying */
763                 v = p_vert_copy(chart, v);
764                 v->edge = lastwe;
765
766                 we = lastwe;
767                 do {
768                         we->vert = v;
769                         we = p_wheel_edge_next(we);
770                 } while (we && (we != lastwe));
771         }
772 }
773
774 static PChart **p_split_charts(PHandle *handle, PChart *chart, int ncharts)
775 {
776         PChart **charts = MEM_mallocN(sizeof*charts * ncharts, "PCharts"), *nchart;
777         PFace *f, *nextf;
778         int i;
779
780         for (i = 0; i < ncharts; i++)
781                 charts[i] = p_chart_new(handle);
782
783         f = (PFace*)chart->faces->first;
784         while (f) {
785                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
786                 nextf = f->link.next;
787
788                 nchart = charts[f->u.chart];
789
790                 phash_insert(nchart->faces, (PHashLink*)f);
791                 phash_insert(nchart->edges, (PHashLink*)e1);
792                 phash_insert(nchart->edges, (PHashLink*)e2);
793                 phash_insert(nchart->edges, (PHashLink*)e3);
794
795                 p_split_vert(nchart, e1);
796                 p_split_vert(nchart, e2);
797                 p_split_vert(nchart, e3);
798
799                 f = nextf;
800         }
801
802         return charts;
803 }
804
805 static void p_face_backup_uvs(PFace *f)
806 {
807         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
808
809         e1->old_uv[0] = e1->orig_uv[0];
810         e1->old_uv[1] = e1->orig_uv[1];
811         e2->old_uv[0] = e2->orig_uv[0];
812         e2->old_uv[1] = e2->orig_uv[1];
813         e3->old_uv[0] = e3->orig_uv[0];
814         e3->old_uv[1] = e3->orig_uv[1];
815 }
816
817 static void p_face_restore_uvs(PFace *f)
818 {
819         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
820
821         e1->orig_uv[0] = e1->old_uv[0];
822         e1->orig_uv[1] = e1->old_uv[1];
823         e2->orig_uv[0] = e2->old_uv[0];
824         e2->orig_uv[1] = e2->old_uv[1];
825         e3->orig_uv[0] = e3->old_uv[0];
826         e3->orig_uv[1] = e3->old_uv[1];
827 }
828
829 static PFace *p_face_add(PChart *chart, ParamKey key, ParamKey *vkeys,
830                          float *co[3], float *uv[3], int i1, int i2, int i3,
831                          ParamBool *pin, ParamBool *select)
832 {
833         PFace *f;
834         PEdge *e1, *e2, *e3;
835
836         /* allocate */
837         f = (PFace*)BLI_memarena_alloc(chart->handle->arena, sizeof *f);
838         f->flag=0; // init !
839
840         e1 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e1);
841         e2 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e2);
842         e3 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e3);
843
844         
845
846
847         /* set up edges */
848         f->edge = e1;
849         e1->face = e2->face = e3->face = f;
850
851         e1->next = e2;
852         e2->next = e3;
853         e3->next = e1;
854
855         e1->pair = NULL;
856         e2->pair = NULL;
857         e3->pair = NULL;
858    
859         e1->flag =0;
860         e2->flag =0;
861         e3->flag =0;
862
863
864         if (co && uv) {
865                 e1->vert = p_vert_lookup(chart, vkeys[i1], co[i1], e1);
866                 e2->vert = p_vert_lookup(chart, vkeys[i2], co[i2], e2);
867                 e3->vert = p_vert_lookup(chart, vkeys[i3], co[i3], e3);
868
869                 e1->orig_uv = uv[i1];
870                 e2->orig_uv = uv[i2];
871                 e3->orig_uv = uv[i3];
872
873         }
874         else {
875                 /* internal call to add face */
876                 e1->vert = e2->vert = e3->vert = NULL;
877                 e1->orig_uv = e2->orig_uv = e3->orig_uv = NULL;
878         }
879
880         if (pin) {
881                 if (pin[i1]) e1->flag |= PEDGE_PIN;
882                 if (pin[i2]) e2->flag |= PEDGE_PIN;
883                 if (pin[i3]) e3->flag |= PEDGE_PIN;
884         }
885
886         if (select) {
887                 if (select[i1]) e1->flag |= PEDGE_SELECT;
888                 if (select[i2]) e2->flag |= PEDGE_SELECT;
889                 if (select[i3]) e3->flag |= PEDGE_SELECT;
890         }
891
892         /* insert into hash */
893         f->link.key = key;
894         phash_insert(chart->faces, (PHashLink*)f);
895
896         e1->link.key = vkeys[i1]^vkeys[i2];
897         e2->link.key = vkeys[i2]^vkeys[i3];
898         e3->link.key = vkeys[i3]^vkeys[i1];
899
900         phash_insert(chart->edges, (PHashLink*)e1);
901         phash_insert(chart->edges, (PHashLink*)e2);
902         phash_insert(chart->edges, (PHashLink*)e3);
903
904         return f;
905 }
906
907 static PBool p_quad_split_direction(float **co)
908 {
909     float a1, a2;
910         
911         a1 = p_vec_angle_cos(co[0], co[1], co[2]);
912         a1 += p_vec_angle_cos(co[1], co[0], co[2]);
913         a1 += p_vec_angle_cos(co[2], co[0], co[1]);
914
915         a2 = p_vec_angle_cos(co[0], co[1], co[3]);
916         a2 += p_vec_angle_cos(co[1], co[0], co[3]);
917         a2 += p_vec_angle_cos(co[3], co[0], co[1]);
918
919         return (a1 > a2);
920 }
921
922 static float p_edge_length(PEdge *e)
923 {
924     PVert *v1 = e->vert, *v2 = e->next->vert;
925     float d[3];
926
927     d[0] = v2->co[0] - v1->co[0];
928     d[1] = v2->co[1] - v1->co[1];
929     d[2] = v2->co[2] - v1->co[2];
930
931     return sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
932 }
933
934 static float p_edge_uv_length(PEdge *e)
935 {
936     PVert *v1 = e->vert, *v2 = e->next->vert;
937     float d[3];
938
939     d[0] = v2->uv[0] - v1->uv[0];
940     d[1] = v2->uv[1] - v1->uv[1];
941
942     return sqrt(d[0]*d[0] + d[1]*d[1]);
943 }
944
945 void p_chart_uv_bbox(PChart *chart, float *minv, float *maxv)
946 {
947         PVert *v;
948
949         INIT_MINMAX2(minv, maxv);
950
951         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
952                 DO_MINMAX2(v->uv, minv, maxv);
953         }
954 }
955
956 static void p_chart_uv_scale(PChart *chart, float scale)
957 {
958         PVert *v;
959
960         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
961                 v->uv[0] *= scale;
962                 v->uv[1] *= scale;
963         }
964 }
965
966 static void p_chart_uv_translate(PChart *chart, float trans[2])
967 {
968         PVert *v;
969
970         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
971                 v->uv[0] += trans[0];
972                 v->uv[1] += trans[1];
973         }
974 }
975
976 static void p_chart_boundaries(PChart *chart, int *nboundaries, PEdge **outer)
977 {   
978     PEdge *e, *be;
979     float len, maxlen = -1.0;
980
981         *nboundaries = 0;
982         *outer = NULL;
983
984         for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
985         if (e->pair || (e->flag & PEDGE_DONE))
986             continue;
987
988                 (*nboundaries)++;
989         len = 0.0f;
990     
991                 be = e;
992                 do {
993             be->flag |= PEDGE_DONE;
994             len += p_edge_length(be);
995                         be = be->next->vert->edge;
996         } while(be != e);
997
998         if (len > maxlen) {
999                         *outer = e;
1000             maxlen = len;
1001         }
1002     }
1003
1004         for (e=(PEdge*)chart->edges->first; e; e=e->link.next)
1005         e->flag &= ~PEDGE_DONE;
1006 }
1007
1008 static float p_edge_boundary_angle(PEdge *e)
1009 {
1010         PEdge *we;
1011         PVert *v, *v1, *v2;
1012         float angle;
1013         int n = 0;
1014
1015         v = e->vert;
1016
1017         /* concave angle check -- could be better */
1018         angle = M_PI;
1019
1020         we = v->edge;
1021         do {
1022                 v1 = we->next->vert;
1023                 v2 = we->next->next->vert;
1024                 angle -= p_vec_angle(v1->co, v->co, v2->co);
1025
1026                 we = we->next->next->pair;
1027                 n++;
1028         } while (we && (we != v->edge));
1029
1030         return angle;
1031 }
1032
1033 static PEdge *p_boundary_edge_next(PEdge *e)
1034 {
1035         return e->next->vert->edge;
1036 }
1037
1038 static PEdge *p_boundary_edge_prev(PEdge *e)
1039 {
1040     PEdge *we = e, *last;
1041
1042         do {
1043                 last = we;
1044                 we = p_wheel_edge_next(we);
1045         } while (we && (we != e));
1046
1047         return last->next->next;
1048 }
1049
1050 static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
1051 {
1052         PEdge *e, *e1, *e2;
1053         PHashKey vkeys[3];
1054         PFace *f;
1055         struct PHeap *heap = pheap_new(nedges);
1056         float angle;
1057
1058         e = be;
1059         do {
1060                 angle = p_edge_boundary_angle(e);
1061                 e->u.heaplink = pheap_insert(heap, angle, e);
1062
1063                 e = e->next->vert->edge;
1064         } while(e != be);
1065
1066         if (nedges == 2) {
1067                 /* no real boundary, but an isolated seam */
1068                 e = be->next->vert->edge;
1069                 e->pair = be;
1070                 be->pair = e;
1071
1072                 pheap_remove(heap, e->u.heaplink);
1073                 pheap_remove(heap, be->u.heaplink);
1074         }
1075         else {
1076                 while (nedges > 2) {
1077                         PEdge *ne, *ne1, *ne2;
1078
1079                         e = pheap_popmin(heap);
1080
1081                         e1 = p_boundary_edge_prev(e);
1082                         e2 = p_boundary_edge_next(e);
1083
1084                         pheap_remove(heap, e1->u.heaplink);
1085                         pheap_remove(heap, e2->u.heaplink);
1086                         e->u.heaplink = e1->u.heaplink = e2->u.heaplink = NULL;
1087
1088                         e->flag |= PEDGE_FILLED;
1089                         e1->flag |= PEDGE_FILLED;
1090
1091                         vkeys[0] = e->vert->link.key;
1092                         vkeys[1] = e1->vert->link.key;
1093                         vkeys[2] = e2->vert->link.key;
1094
1095                         f = p_face_add(chart, -1, vkeys, NULL, NULL, 0, 1, 2, NULL, NULL);
1096                         f->flag |= PFACE_FILLED;
1097
1098                         ne = f->edge->next->next;
1099                         ne1 = f->edge;
1100                         ne2 = f->edge->next;
1101
1102                         ne->flag = ne1->flag = ne2->flag = PEDGE_FILLED;
1103
1104                         e->pair = ne;
1105                         ne->pair = e;
1106                         e1->pair = ne1;
1107                         ne1->pair = e1;
1108
1109                         ne->vert = e2->vert;
1110                         ne1->vert = e->vert;
1111                         ne2->vert = e1->vert;
1112
1113                         if (nedges == 3) {
1114                                 e2->pair = ne2;
1115                                 ne2->pair = e2;
1116                         }
1117                         else {
1118                                 ne2->vert->edge = ne2;
1119                                 
1120                                 ne2->u.heaplink = pheap_insert(heap, p_edge_boundary_angle(ne2), ne2);
1121                                 e2->u.heaplink = pheap_insert(heap, p_edge_boundary_angle(e2), e2);
1122                         }
1123
1124                         nedges--;
1125                 }
1126         }
1127
1128         pheap_delete(heap);
1129 }
1130
1131 static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
1132 {
1133         PEdge *e, *enext, *be;
1134         int nedges;
1135
1136         for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
1137                 enext = e->link.next;
1138
1139         if (e->pair || (e->flag & PEDGE_FILLED))
1140             continue;
1141
1142                 nedges = 0;
1143                 be = e;
1144                 do {
1145                         be->flag |= PEDGE_FILLED;
1146                         be = be->next->vert->edge;
1147                         nedges++;
1148                 } while(be != e);
1149
1150                 if (e != outer)
1151                         p_chart_fill_boundary(chart, e, nedges);
1152     }
1153 }
1154
1155 static void p_flush_uvs(PChart *chart)
1156 {
1157         PEdge *e;
1158
1159         for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
1160                 if (e->orig_uv) {
1161                         e->orig_uv[0] = e->vert->uv[0];
1162                         e->orig_uv[1] = e->vert->uv[1];
1163                 }
1164         }
1165 }
1166
1167 static void p_flush_uvs_blend(PChart *chart, float blend)
1168 {
1169         PEdge *e;
1170         float invblend = 1.0f - blend;
1171
1172         for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
1173                 if (e->orig_uv) {
1174                         e->orig_uv[0] = blend*e->old_uv[0] + invblend*e->vert->uv[0];
1175                         e->orig_uv[1] = blend*e->old_uv[1] + invblend*e->vert->uv[1];
1176                 }
1177         }
1178 }
1179
1180 /* Exported */
1181
1182 ParamHandle *param_construct_begin()
1183 {
1184         PHandle *handle = MEM_callocN(sizeof*handle, "PHandle");
1185         handle->construction_chart = p_chart_new(handle);
1186         handle->state = PHANDLE_STATE_ALLOCATED;
1187         handle->arena = BLI_memarena_new((1<<16));
1188         
1189         return (ParamHandle*)handle;
1190 }
1191
1192 void param_delete(ParamHandle *handle)
1193 {
1194         PHandle *phandle = (PHandle*)handle;
1195         int i;
1196
1197         param_assert((phandle->state == PHANDLE_STATE_ALLOCATED) ||
1198                      (phandle->state == PHANDLE_STATE_CONSTRUCTED));
1199
1200         for (i = 0; i < phandle->ncharts; i++)
1201                 p_chart_delete(phandle->charts[i]);
1202         
1203         if (phandle->charts)
1204                 MEM_freeN(phandle->charts);
1205
1206         if (phandle->construction_chart)
1207                 p_chart_delete(phandle->construction_chart);
1208
1209         BLI_memarena_free(phandle->arena);
1210         MEM_freeN(phandle);
1211 }
1212
1213 void param_face_add(ParamHandle *handle, ParamKey key, int nverts,
1214                     ParamKey *vkeys, float **co, float **uv,
1215                     ParamBool *pin, ParamBool *select)
1216 {
1217         PHandle *phandle = (PHandle*)handle;
1218         PChart *chart = phandle->construction_chart;
1219
1220         param_assert(phash_lookup(chart->faces, key) == NULL);
1221         param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
1222         param_assert((nverts == 3) || (nverts == 4));
1223
1224         if (nverts == 4) {
1225                 if (!p_quad_split_direction(co)) {
1226                         p_face_add(chart, key, vkeys, co, uv, 0, 1, 2, pin, select);
1227                         p_face_add(chart, key, vkeys, co, uv, 0, 2, 3, pin, select);
1228                 }
1229                 else {
1230                         p_face_add(chart, key, vkeys, co, uv, 0, 1, 3, pin, select);
1231                         p_face_add(chart, key, vkeys, co, uv, 1, 2, 3, pin, select);
1232                 }
1233         }
1234         else
1235                 p_face_add(chart, key, vkeys, co, uv, 0, 1, 2, pin, select);
1236 }
1237
1238 void param_edge_set_seam(ParamHandle *handle, ParamKey *vkeys)
1239 {
1240         PHandle *phandle = (PHandle*)handle;
1241         PChart *chart = phandle->construction_chart;
1242         PEdge *e;
1243
1244         param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
1245
1246         e = p_edge_lookup(chart, vkeys);
1247         if (e)
1248                 e->flag |= PEDGE_SEAM;
1249 }
1250
1251 void param_construct_end(ParamHandle *handle, ParamBool fill, ParamBool impl)
1252 {
1253         PHandle *phandle = (PHandle*)handle;
1254         PChart *chart = phandle->construction_chart;
1255         int i, j, nboundaries = 0;
1256         PEdge *outer;
1257
1258         param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
1259
1260         phandle->ncharts = p_connect_pairs(chart, impl);
1261         phandle->charts = p_split_charts(phandle, chart, phandle->ncharts);
1262
1263         p_chart_delete(chart);
1264         phandle->construction_chart = NULL;
1265
1266         for (i = j = 0; i < phandle->ncharts; i++) {
1267                 p_chart_boundaries(phandle->charts[i], &nboundaries, &outer);
1268
1269                 if (nboundaries == 0) {
1270                         p_chart_delete(phandle->charts[i]);
1271                         continue;
1272                 }
1273
1274                 phandle->charts[j] = phandle->charts[i];
1275                 j++;
1276
1277 #if 0
1278                 if (fill && (nboundaries > 1))
1279                         p_chart_fill_boundaries(phandle->charts[i], outer);
1280 #endif
1281         }
1282
1283         phandle->ncharts = j;
1284
1285         phandle->state = PHANDLE_STATE_CONSTRUCTED;
1286 }
1287
1288 /* Least Squares Conformal Maps */
1289
1290 static void p_chart_lscm_load_solution(PChart *chart)
1291 {
1292         PVert *v;
1293
1294         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
1295                 v->uv[0] = nlGetVariable(2*v->u.index);
1296                 v->uv[1] = nlGetVariable(2*v->u.index + 1);
1297         }
1298 }
1299
1300 static void p_chart_lscm_begin(PChart *chart, PBool live)
1301 {
1302         PVert *v, *pin1, *pin2;
1303         PBool select = P_FALSE;
1304         int npins = 0, id = 0;
1305
1306         /* give vertices matrix indices and count pins */
1307         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
1308                 p_vert_load_pin_select_uvs(v);
1309
1310                 if (v->flag & PVERT_PIN)
1311                         npins++;
1312
1313                 if (v->flag & PVERT_SELECT)
1314                         select = P_TRUE;
1315
1316                 v->u.index = id++;
1317         }
1318
1319         if ((live && !select) || (npins == 1)) {
1320                 chart->u.lscm.context = NULL;
1321         }
1322         else {
1323                 if (npins <= 1) {
1324                         /* not enough pins, lets find some ourself */
1325                         p_extrema_verts(chart, &pin1, &pin2);
1326
1327                         chart->u.lscm.pin1 = pin1;
1328                         chart->u.lscm.pin2 = pin2;
1329                 }
1330                 else {
1331                         chart->flag |= PCHART_NOPACK;
1332                 }
1333
1334                 nlNewContext();
1335                 nlSolverParameteri(NL_NB_VARIABLES, 2*phash_size(chart->verts));
1336                 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
1337
1338                 chart->u.lscm.context = nlGetCurrent();
1339         }
1340 }
1341
1342 static PBool p_chart_lscm_solve(PChart *chart)
1343 {
1344         PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
1345         PFace *f;
1346
1347         nlMakeCurrent(chart->u.lscm.context);
1348
1349         nlBegin(NL_SYSTEM);
1350
1351         for (v=(PVert*)chart->verts->first; v; v=v->link.next)
1352                 if (v->flag & PVERT_PIN)
1353                         p_vert_load_pin_select_uvs(v);
1354
1355         if (chart->u.lscm.pin1) {
1356                 nlLockVariable(2*pin1->u.index);
1357                 nlLockVariable(2*pin1->u.index + 1);
1358                 nlLockVariable(2*pin2->u.index);
1359                 nlLockVariable(2*pin2->u.index + 1);
1360         
1361                 nlSetVariable(2*pin1->u.index, pin1->uv[0]);
1362                 nlSetVariable(2*pin1->u.index + 1, pin1->uv[1]);
1363                 nlSetVariable(2*pin2->u.index, pin2->uv[0]);
1364                 nlSetVariable(2*pin2->u.index + 1, pin2->uv[1]);
1365         }
1366         else {
1367                 /* set and lock the pins */
1368                 for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
1369                         if (v->flag & PVERT_PIN) {
1370                                 nlLockVariable(2*v->u.index);
1371                                 nlLockVariable(2*v->u.index + 1);
1372
1373                                 nlSetVariable(2*v->u.index, v->uv[0]);
1374                                 nlSetVariable(2*v->u.index + 1, v->uv[1]);
1375                         }
1376                 }
1377         }
1378
1379         /* construct matrix */
1380
1381         nlBegin(NL_MATRIX);
1382
1383         for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
1384                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1385                 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
1386                 float a1, a2, a3, ratio, cosine, sine;
1387                 float sina1, sina2, sina3, sinmax;
1388
1389                 if (chart->u.lscm.abf_alpha) {
1390                         /* use abf angles if passed on */
1391                         a1 = *(chart->u.lscm.abf_alpha++);
1392                         a2 = *(chart->u.lscm.abf_alpha++);
1393                         a3 = *(chart->u.lscm.abf_alpha++);
1394                 }
1395                 else
1396                         p_face_angles(f, &a1, &a2, &a3);
1397
1398                 sina1 = sin(a1);
1399                 sina2 = sin(a2);
1400                 sina3 = sin(a3);
1401
1402                 sinmax = MAX3(sina1, sina2, sina3);
1403
1404                 /* shift vertices to find most stable order */
1405                 #define SHIFT3(type, a, b, c) \
1406                         { type tmp; tmp = a; a = c; c = b; b = tmp; }
1407
1408                 if (sina3 != sinmax) {
1409                         SHIFT3(PVert*, v1, v2, v3);
1410                         SHIFT3(float, a1, a2, a3);
1411                         SHIFT3(float, sina1, sina2, sina3);
1412
1413                         if (sina2 == sinmax) {
1414                                 SHIFT3(PVert*, v1, v2, v3);
1415                                 SHIFT3(float, a1, a2, a3);
1416                                 SHIFT3(float, sina1, sina2, sina3);
1417                         }
1418                 }
1419
1420                 /* angle based lscm formulation */
1421                 ratio = (sina3 == 0.0f)? 0.0f: sina2/sina3;
1422                 cosine = cos(a1)*ratio;
1423                 sine = sina1*ratio;
1424
1425                 nlBegin(NL_ROW);
1426                 nlCoefficient(2*v1->u.index,   cosine - 1.0);
1427                 nlCoefficient(2*v1->u.index+1, -sine);
1428                 nlCoefficient(2*v2->u.index,   -cosine);
1429                 nlCoefficient(2*v2->u.index+1, sine);
1430                 nlCoefficient(2*v3->u.index,   1.0);
1431                 nlEnd(NL_ROW);
1432
1433                 nlBegin(NL_ROW);
1434                 nlCoefficient(2*v1->u.index,   sine);
1435                 nlCoefficient(2*v1->u.index+1, cosine - 1.0);
1436                 nlCoefficient(2*v2->u.index,   -sine);
1437                 nlCoefficient(2*v2->u.index+1, -cosine);
1438                 nlCoefficient(2*v3->u.index+1, 1.0);
1439                 nlEnd(NL_ROW);
1440         }
1441
1442         nlEnd(NL_MATRIX);
1443
1444         nlEnd(NL_SYSTEM);
1445
1446         if (nlSolveAdvanced(NULL, NL_TRUE)) {
1447                 p_chart_lscm_load_solution(chart);
1448                 return P_TRUE;
1449         }
1450
1451         return P_FALSE;
1452 }
1453
1454 static void p_chart_lscm_end(PChart *chart)
1455 {
1456         if (chart->u.lscm.context)
1457                 nlDeleteContext(chart->u.lscm.context);
1458
1459         chart->u.lscm.context = NULL;
1460         chart->u.lscm.pin1 = NULL;
1461         chart->u.lscm.pin2 = NULL;
1462 }
1463
1464 void param_lscm_begin(ParamHandle *handle, ParamBool live)
1465 {
1466         PHandle *phandle = (PHandle*)handle;
1467         int i;
1468
1469         param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
1470         phandle->state = PHANDLE_STATE_LSCM;
1471
1472         for (i = 0; i < phandle->ncharts; i++)
1473                 p_chart_lscm_begin(phandle->charts[i], live);
1474 }
1475
1476 void param_lscm_solve(ParamHandle *handle)
1477 {
1478         PHandle *phandle = (PHandle*)handle;
1479         PChart *chart;
1480         int i;
1481         PBool result;
1482
1483         param_assert(phandle->state == PHANDLE_STATE_LSCM);
1484
1485         for (i = 0; i < phandle->ncharts; i++) {
1486                 chart = phandle->charts[i];
1487
1488                 if (chart->u.lscm.context) {
1489                         result = p_chart_lscm_solve(chart);
1490
1491                         if (!result || (chart->u.lscm.pin1))
1492                                 p_chart_lscm_end(chart);
1493                 }
1494         }
1495 }
1496
1497 void param_lscm_end(ParamHandle *handle)
1498 {
1499         PHandle *phandle = (PHandle*)handle;
1500         int i;
1501
1502         param_assert(phandle->state == PHANDLE_STATE_LSCM);
1503
1504         for (i = 0; i < phandle->ncharts; i++)
1505                 p_chart_lscm_end(phandle->charts[i]);
1506
1507         phandle->state = PHANDLE_STATE_CONSTRUCTED;
1508 }
1509
1510 /* Stretch */
1511
1512 #define P_STRETCH_ITER 20
1513
1514 static void p_stretch_pin_boundary(PChart *chart)
1515 {
1516         PVert *v;
1517
1518         for(v=(PVert*)chart->verts->first; v; v=v->link.next)
1519                 if (v->edge->pair == NULL)
1520                         v->flag |= PVERT_PIN;
1521                 else
1522                         v->flag &= ~PVERT_PIN;
1523 }
1524
1525 static float p_face_stretch(PFace *f)
1526 {
1527         float T, w, tmp[3];
1528         float Ps[3], Pt[3];
1529         float a, c, area;
1530         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1531         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
1532
1533         area = p_face_uv_area_signed(f);
1534
1535         if (area <= 0.0f) /* flipped face -> infinite stretch */
1536                 return 1e10f;
1537         
1538         if (f->flag & PFACE_FILLED)
1539                 return 0.0f;
1540
1541         w= 1.0f/(2.0f*area);
1542
1543         /* compute derivatives */
1544         VecCopyf(Ps, v1->co);
1545         VecMulf(Ps, (v2->uv[1] - v3->uv[1]));
1546
1547         VecCopyf(tmp, v2->co);
1548         VecMulf(tmp, (v3->uv[1] - v1->uv[1]));
1549         VecAddf(Ps, Ps, tmp);
1550
1551         VecCopyf(tmp, v3->co);
1552         VecMulf(tmp, (v1->uv[1] - v2->uv[1]));
1553         VecAddf(Ps, Ps, tmp);
1554
1555         VecMulf(Ps, w);
1556
1557         VecCopyf(Pt, v1->co);
1558         VecMulf(Pt, (v3->uv[0] - v2->uv[0]));
1559
1560         VecCopyf(tmp, v2->co);
1561         VecMulf(tmp, (v1->uv[0] - v3->uv[0]));
1562         VecAddf(Pt, Pt, tmp);
1563
1564         VecCopyf(tmp, v3->co);
1565         VecMulf(tmp, (v2->uv[0] - v1->uv[0]));
1566         VecAddf(Pt, Pt, tmp);
1567
1568         VecMulf(Pt, w);
1569
1570         /* Sander Tensor */
1571         a= Inpf(Ps, Ps);
1572         c= Inpf(Pt, Pt);
1573
1574         T = sqrt(0.5f*(a + c)*f->u.area3d);
1575
1576         return T;
1577 }
1578
1579 static float p_stretch_compute_vertex(PVert *v)
1580 {
1581         PEdge *e = v->edge;
1582         float sum = 0.0f;
1583
1584         do {
1585                 sum += p_face_stretch(e->face);
1586                 e = p_wheel_edge_next(e);
1587         } while (e && e != (v->edge));
1588
1589         return sum;
1590 }
1591
1592 static void p_chart_stretch_minimize(PChart *chart, RNG *rng)
1593 {
1594         PVert *v;
1595         PEdge *e;
1596         int j, nedges;
1597         float orig_stretch, low, stretch_low, high, stretch_high, mid, stretch;
1598         float orig_uv[2], dir[2], random_angle, trusted_radius;
1599
1600         for(v=(PVert*)chart->verts->first; v; v=v->link.next) {
1601                 if((v->flag & PVERT_PIN) || !(v->flag & PVERT_SELECT))
1602                         continue;
1603
1604                 orig_stretch = p_stretch_compute_vertex(v);
1605                 orig_uv[0] = v->uv[0];
1606                 orig_uv[1] = v->uv[1];
1607
1608                 /* move vertex in a random direction */
1609                 trusted_radius = 0.0f;
1610                 nedges = 0;
1611                 e = v->edge;
1612
1613                 do {
1614                         trusted_radius += p_edge_uv_length(e);
1615                         nedges++;
1616
1617                         e = p_wheel_edge_next(e);
1618                 } while (e && e != (v->edge));
1619
1620                 trusted_radius /= 2 * nedges;
1621
1622                 random_angle = rng_getFloat(rng) * 2.0 * M_PI;
1623                 dir[0] = trusted_radius * cos(random_angle);
1624                 dir[1] = trusted_radius * sin(random_angle);
1625
1626                 /* calculate old and new stretch */
1627                 low = 0;
1628                 stretch_low = orig_stretch;
1629
1630                 Vec2Addf(v->uv, orig_uv, dir);
1631                 high = 1;
1632                 stretch = stretch_high = p_stretch_compute_vertex(v);
1633
1634                 /* binary search for lowest stretch position */
1635                 for (j = 0; j < P_STRETCH_ITER; j++) {
1636                         mid = 0.5 * (low + high);
1637                         v->uv[0]= orig_uv[0] + mid*dir[0];
1638                         v->uv[1]= orig_uv[1] + mid*dir[1];
1639                         stretch = p_stretch_compute_vertex(v);
1640
1641                         if (stretch_low < stretch_high) {
1642                                 high = mid;
1643                                 stretch_high = stretch;
1644                         }
1645                         else {
1646                                 low = mid;
1647                                 stretch_low = stretch;
1648                         }
1649                 }
1650
1651                 /* no luck, stretch has increased, reset to old values */
1652                 if(stretch >= orig_stretch)
1653                         Vec2Copyf(v->uv, orig_uv);
1654         }
1655 }
1656
1657 void param_stretch_begin(ParamHandle *handle)
1658 {
1659         PHandle *phandle = (PHandle*)handle;
1660         PChart *chart;
1661         PVert *v;
1662         PFace *f;
1663         int i;
1664
1665         param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
1666         phandle->state = PHANDLE_STATE_STRETCH;
1667
1668         phandle->rng = rng_new(31415926);
1669         phandle->blend = 0.0f;
1670
1671         for (i = 0; i < phandle->ncharts; i++) {
1672                 chart = phandle->charts[i];
1673
1674                 for (v=(PVert*)chart->verts->first; v; v=v->link.next)
1675                         p_vert_load_select_uvs(v);
1676
1677                 p_stretch_pin_boundary(chart);
1678
1679                 for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
1680                         p_face_backup_uvs(f);
1681                         f->u.area3d = p_face_area(f);
1682                 }
1683         }
1684 }
1685
1686 void param_stretch_blend(ParamHandle *handle, float blend)
1687 {
1688         PHandle *phandle = (PHandle*)handle;
1689
1690         param_assert(phandle->state == PHANDLE_STATE_STRETCH);
1691         phandle->blend = blend;
1692 }
1693
1694 void param_stretch_iter(ParamHandle *handle)
1695 {
1696         PHandle *phandle = (PHandle*)handle;
1697         PChart *chart;
1698         int i;
1699
1700         param_assert(phandle->state == PHANDLE_STATE_STRETCH);
1701
1702         for (i = 0; i < phandle->ncharts; i++) {
1703                 chart = phandle->charts[i];
1704                 p_chart_stretch_minimize(chart, phandle->rng);
1705         }
1706 }
1707
1708 void param_stretch_end(ParamHandle *handle)
1709 {
1710         PHandle *phandle = (PHandle*)handle;
1711
1712         param_assert(phandle->state == PHANDLE_STATE_STRETCH);
1713         phandle->state = PHANDLE_STATE_CONSTRUCTED;
1714
1715         rng_free(phandle->rng);
1716         phandle->rng = NULL;
1717 }
1718
1719 /* Flushing */
1720
1721 void param_flush(ParamHandle *handle)
1722 {
1723         PHandle *phandle = (PHandle*)handle;
1724         PChart *chart;
1725         int i;
1726
1727         for (i = 0; i < phandle->ncharts; i++) {
1728                 chart = phandle->charts[i];
1729
1730                 if ((phandle->state == PHANDLE_STATE_LSCM) && !chart->u.lscm.context)
1731                         continue;
1732
1733                 if (phandle->blend == 0.0f)
1734                         p_flush_uvs(chart);
1735                 else
1736                         p_flush_uvs_blend(chart, phandle->blend);
1737         }
1738 }
1739
1740 void param_flush_restore(ParamHandle *handle)
1741 {
1742         PHandle *phandle = (PHandle*)handle;
1743         PChart *chart;
1744         PFace *f;
1745         int i;
1746
1747         for (i = 0; i < phandle->ncharts; i++) {
1748                 chart = phandle->charts[i];
1749
1750                 for (f=(PFace*)chart->faces->first; f; f=f->link.next)
1751                         p_face_restore_uvs(f);
1752         }
1753 }
1754
1755 /* Packing */
1756
1757 static int compare_chart_area(const void *a, const void *b)
1758 {
1759         PChart *ca = *((PChart**)a);
1760         PChart *cb = *((PChart**)b);
1761
1762     if (ca->u.pack.area > cb->u.pack.area)
1763                 return -1;
1764         else if (ca->u.pack.area == cb->u.pack.area)
1765                 return 0;
1766         else
1767                 return 1;
1768 }
1769
1770 static PBool p_pack_try(PHandle *handle, float side)
1771 {
1772         PChart *chart;
1773         float packx, packy, rowh, groupw, w, h;
1774         int i;
1775
1776         packx= packy= 0.0;
1777         rowh= 0.0;
1778         groupw= 1.0/sqrt(handle->ncharts);
1779
1780         for (i = 0; i < handle->ncharts; i++) {
1781                 chart = handle->charts[i];
1782
1783                 if (chart->flag & PCHART_NOPACK)
1784                         continue;
1785
1786                 w = chart->u.pack.size[0];
1787                 h = chart->u.pack.size[1];
1788
1789                 if(w <= (side-packx)) {
1790                         chart->u.pack.trans[0] = packx;
1791                         chart->u.pack.trans[1] = packy;
1792
1793                         packx += w;
1794                         rowh= MAX2(rowh, h);
1795                 }
1796                 else {
1797                         packy += rowh;
1798                         packx = w;
1799                         rowh = h;
1800
1801                         chart->u.pack.trans[0] = 0.0;
1802                         chart->u.pack.trans[1] = packy;
1803                 }
1804
1805                 if (packy+rowh > side)
1806                         return P_FALSE;
1807         }
1808
1809         return P_TRUE;
1810 }
1811
1812 #define PACK_SEARCH_DEPTH 15
1813
1814 void param_pack(ParamHandle *handle)
1815 {
1816         PHandle *phandle = (PHandle*)handle;
1817         PChart *chart;
1818         float uv_area, area, trans[2], minside, maxside, totarea, side;
1819         int i;
1820
1821         /* very simple rectangle packing */
1822
1823         if (phandle->ncharts == 0)
1824                 return;
1825
1826         totarea = 0.0f;
1827         maxside = 0.0f;
1828
1829         for (i = 0; i < phandle->ncharts; i++) {
1830                 chart = phandle->charts[i];
1831
1832                 if (chart->flag & PCHART_NOPACK) {
1833                         chart->u.pack.area = 0.0f;
1834                         continue;
1835                 }
1836
1837                 p_chart_area(chart, &uv_area, &area);
1838                 p_chart_uv_bbox(chart, trans, chart->u.pack.size);
1839
1840                 /* translate to origin and make area equal to 3d area */
1841                 chart->u.pack.rescale = (uv_area > 0.0f)? sqrt(area)/sqrt(uv_area): 0.0f;
1842                 chart->u.pack.area = area;
1843                 totarea += area;
1844
1845                 trans[0] = -trans[0];
1846                 trans[1] = -trans[1];
1847                 p_chart_uv_translate(chart, trans);
1848                 p_chart_uv_scale(chart, chart->u.pack.rescale);
1849
1850                 /* compute new dimensions for packing */
1851                 chart->u.pack.size[0] += trans[0];
1852                 chart->u.pack.size[1] += trans[1];
1853                 chart->u.pack.size[0] *= chart->u.pack.rescale;
1854                 chart->u.pack.size[1] *= chart->u.pack.rescale;
1855
1856                 maxside = MAX3(maxside, chart->u.pack.size[0], chart->u.pack.size[1]);
1857         }
1858
1859         /* sort by chart area, largest first */
1860         qsort(phandle->charts, phandle->ncharts, sizeof(PChart*), compare_chart_area);
1861
1862         /* binary search over pack region size */
1863         minside = MAX2(sqrt(totarea), maxside);
1864         maxside = (((int)sqrt(phandle->ncharts-1))+1)*maxside;
1865
1866         if (minside < maxside) { /* should always be true */
1867
1868                 for (i = 0; i < PACK_SEARCH_DEPTH; i++) {
1869                         if (p_pack_try(phandle, (minside+maxside)*0.5f + 1e-5))
1870                                 maxside = (minside+maxside)*0.5f;
1871                         else
1872                                 minside = (minside+maxside)*0.5f;
1873                 }
1874         }
1875
1876         /* do the actual packing */
1877         side = maxside + 1e-5;
1878         if (!p_pack_try(phandle, side))
1879                 param_warning("packing failed.\n");
1880
1881         for (i = 0; i < phandle->ncharts; i++) {
1882                 chart = phandle->charts[i];
1883
1884                 if (chart->flag & PCHART_NOPACK)
1885                         continue;
1886
1887                 p_chart_uv_scale(chart, 1.0f/side);
1888                 trans[0] = chart->u.pack.trans[0]/side;
1889                 trans[1] = chart->u.pack.trans[1]/side;
1890                 p_chart_uv_translate(chart, trans);
1891         }
1892 }
1893