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