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