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