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