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