Cleanup: shadowing (unwrap)
[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         NLContext *context;
2475
2476         context = nlNewContext();
2477         nlSolverParameteri(context, NL_NB_VARIABLES, nvar);
2478
2479         nlBegin(context, NL_SYSTEM);
2480
2481         nlBegin(context, NL_MATRIX);
2482
2483         for (i = 0; i < nvar; i++)
2484                 nlRightHandSideAdd(context, 0, i, sys->bInterior[i]);
2485
2486         for (f = chart->faces; f; f = f->nextlink) {
2487                 float wi1, wi2, wi3, b, si, beta[3], j2[3][3], W[3][3];
2488                 float row1[6], row2[6], row3[6];
2489                 int vid[6];
2490                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
2491                 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
2492
2493                 wi1 = 1.0f / sys->weight[e1->u.id];
2494                 wi2 = 1.0f / sys->weight[e2->u.id];
2495                 wi3 = 1.0f / sys->weight[e3->u.id];
2496
2497                 /* bstar1 = (J1*dInv*bAlpha - bTriangle) */
2498                 b = sys->bAlpha[e1->u.id] * wi1;
2499                 b += sys->bAlpha[e2->u.id] * wi2;
2500                 b += sys->bAlpha[e3->u.id] * wi3;
2501                 b -= sys->bTriangle[f->u.id];
2502
2503                 /* si = J1*d*J1t */
2504                 si = 1.0f / (wi1 + wi2 + wi3);
2505
2506                 /* J1t*si*bstar1 - bAlpha */
2507                 beta[0] = b * si - sys->bAlpha[e1->u.id];
2508                 beta[1] = b * si - sys->bAlpha[e2->u.id];
2509                 beta[2] = b * si - sys->bAlpha[e3->u.id];
2510
2511                 /* use this later for computing other lambda's */
2512                 sys->bstar[f->u.id] = b;
2513                 sys->dstar[f->u.id] = si;
2514
2515                 /* set matrix */
2516                 W[0][0] = si - sys->weight[e1->u.id]; W[0][1] = si; W[0][2] = si;
2517                 W[1][0] = si; W[1][1] = si - sys->weight[e2->u.id]; W[1][2] = si;
2518                 W[2][0] = si; W[2][1] = si; W[2][2] = si - sys->weight[e3->u.id];
2519
2520                 vid[0] = vid[1] = vid[2] = vid[3] = vid[4] = vid[5] = -1;
2521
2522                 if (v1->flag & PVERT_INTERIOR) {
2523                         vid[0] = v1->u.id;
2524                         vid[3] = ninterior + v1->u.id;
2525
2526                         sys->J2dt[e1->u.id][0] = j2[0][0] = 1.0f * wi1;
2527                         sys->J2dt[e2->u.id][0] = j2[1][0] = p_abf_compute_sin_product(sys, v1, e2->u.id) * wi2;
2528                         sys->J2dt[e3->u.id][0] = j2[2][0] = p_abf_compute_sin_product(sys, v1, e3->u.id) * wi3;
2529
2530                         nlRightHandSideAdd(context, 0, v1->u.id, j2[0][0] * beta[0]);
2531                         nlRightHandSideAdd(context, 0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
2532
2533                         row1[0] = j2[0][0] * W[0][0];
2534                         row2[0] = j2[0][0] * W[1][0];
2535                         row3[0] = j2[0][0] * W[2][0];
2536
2537                         row1[3] = j2[1][0] * W[0][1] + j2[2][0] * W[0][2];
2538                         row2[3] = j2[1][0] * W[1][1] + j2[2][0] * W[1][2];
2539                         row3[3] = j2[1][0] * W[2][1] + j2[2][0] * W[2][2];
2540                 }
2541
2542                 if (v2->flag & PVERT_INTERIOR) {
2543                         vid[1] = v2->u.id;
2544                         vid[4] = ninterior + v2->u.id;
2545
2546                         sys->J2dt[e1->u.id][1] = j2[0][1] = p_abf_compute_sin_product(sys, v2, e1->u.id) * wi1;
2547                         sys->J2dt[e2->u.id][1] = j2[1][1] = 1.0f * wi2;
2548                         sys->J2dt[e3->u.id][1] = j2[2][1] = p_abf_compute_sin_product(sys, v2, e3->u.id) * wi3;
2549
2550                         nlRightHandSideAdd(context, 0, v2->u.id, j2[1][1] * beta[1]);
2551                         nlRightHandSideAdd(context, 0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
2552
2553                         row1[1] = j2[1][1] * W[0][1];
2554                         row2[1] = j2[1][1] * W[1][1];
2555                         row3[1] = j2[1][1] * W[2][1];
2556
2557                         row1[4] = j2[0][1] * W[0][0] + j2[2][1] * W[0][2];
2558                         row2[4] = j2[0][1] * W[1][0] + j2[2][1] * W[1][2];
2559                         row3[4] = j2[0][1] * W[2][0] + j2[2][1] * W[2][2];
2560                 }
2561
2562                 if (v3->flag & PVERT_INTERIOR) {
2563                         vid[2] = v3->u.id;
2564                         vid[5] = ninterior + v3->u.id;
2565
2566                         sys->J2dt[e1->u.id][2] = j2[0][2] = p_abf_compute_sin_product(sys, v3, e1->u.id) * wi1;
2567                         sys->J2dt[e2->u.id][2] = j2[1][2] = p_abf_compute_sin_product(sys, v3, e2->u.id) * wi2;
2568                         sys->J2dt[e3->u.id][2] = j2[2][2] = 1.0f * wi3;
2569
2570                         nlRightHandSideAdd(context, 0, v3->u.id, j2[2][2] * beta[2]);
2571                         nlRightHandSideAdd(context, 0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
2572
2573                         row1[2] = j2[2][2] * W[0][2];
2574                         row2[2] = j2[2][2] * W[1][2];
2575                         row3[2] = j2[2][2] * W[2][2];
2576
2577                         row1[5] = j2[0][2] * W[0][0] + j2[1][2] * W[0][1];
2578                         row2[5] = j2[0][2] * W[1][0] + j2[1][2] * W[1][1];
2579                         row3[5] = j2[0][2] * W[2][0] + j2[1][2] * W[2][1];
2580                 }
2581
2582                 for (i = 0; i < 3; i++) {
2583                         int r = vid[i];
2584
2585                         if (r == -1)
2586                                 continue;
2587
2588                         for (j = 0; j < 6; j++) {
2589                                 int c = vid[j];
2590
2591                                 if (c == -1)
2592                                         continue;
2593
2594                                 if (i == 0)
2595                                         nlMatrixAdd(context, r, c, j2[0][i] * row1[j]);
2596                                 else
2597                                         nlMatrixAdd(context, r + ninterior, c, j2[0][i] * row1[j]);
2598
2599                                 if (i == 1)
2600                                         nlMatrixAdd(context, r, c, j2[1][i] * row2[j]);
2601                                 else
2602                                         nlMatrixAdd(context, r + ninterior, c, j2[1][i] * row2[j]);
2603
2604
2605                                 if (i == 2)
2606                                         nlMatrixAdd(context, r, c, j2[2][i] * row3[j]);
2607                                 else
2608                                         nlMatrixAdd(context, r + ninterior, c, j2[2][i] * row3[j]);
2609                         }
2610                 }
2611         }
2612
2613         nlEnd(context, NL_MATRIX);
2614
2615         nlEnd(context, NL_SYSTEM);
2616
2617         success = nlSolve(context, NL_FALSE);
2618
2619         if (success) {
2620                 for (f = chart->faces; f; f = f->nextlink) {
2621                         float dlambda1, pre[3], dalpha;
2622                         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
2623                         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
2624
2625                         pre[0] = pre[1] = pre[2] = 0.0;
2626
2627                         if (v1->flag & PVERT_INTERIOR) {
2628                                 float x = nlGetVariable(context, 0, v1->u.id);
2629                                 float x2 = nlGetVariable(context, 0, ninterior + v1->u.id);
2630                                 pre[0] += sys->J2dt[e1->u.id][0] * x;
2631                                 pre[1] += sys->J2dt[e2->u.id][0] * x2;
2632                                 pre[2] += sys->J2dt[e3->u.id][0] * x2;
2633                         }
2634
2635                         if (v2->flag & PVERT_INTERIOR) {
2636                                 float x = nlGetVariable(context, 0, v2->u.id);
2637                                 float x2 = nlGetVariable(context, 0, ninterior + v2->u.id);
2638                                 pre[0] += sys->J2dt[e1->u.id][1] * x2;
2639                                 pre[1] += sys->J2dt[e2->u.id][1] * x;
2640                                 pre[2] += sys->J2dt[e3->u.id][1] * x2;
2641                         }
2642
2643                         if (v3->flag & PVERT_INTERIOR) {
2644                                 float x = nlGetVariable(context, 0, v3->u.id);
2645                                 float x2 = nlGetVariable(context, 0, ninterior + v3->u.id);
2646                                 pre[0] += sys->J2dt[e1->u.id][2] * x2;
2647                                 pre[1] += sys->J2dt[e2->u.id][2] * x2;
2648                                 pre[2] += sys->J2dt[e3->u.id][2] * x;
2649                         }
2650
2651                         dlambda1 = pre[0] + pre[1] + pre[2];
2652                         dlambda1 = sys->dstar[f->u.id] * (sys->bstar[f->u.id] - dlambda1);
2653                         
2654                         sys->lambdaTriangle[f->u.id] += dlambda1;
2655
2656                         dalpha = (sys->bAlpha[e1->u.id] - dlambda1);
2657                         sys->alpha[e1->u.id] += dalpha / sys->weight[e1->u.id] - pre[0];
2658
2659                         dalpha = (sys->bAlpha[e2->u.id] - dlambda1);
2660                         sys->alpha[e2->u.id] += dalpha / sys->weight[e2->u.id] - pre[1];
2661
2662                         dalpha = (sys->bAlpha[e3->u.id] - dlambda1);
2663                         sys->alpha[e3->u.id] += dalpha / sys->weight[e3->u.id] - pre[2];
2664
2665                         /* clamp */
2666                         e = f->edge;
2667                         do {
2668                                 if (sys->alpha[e->u.id] > (float)M_PI)
2669                                         sys->alpha[e->u.id] = (float)M_PI;
2670                                 else if (sys->alpha[e->u.id] < 0.0f)
2671                                         sys->alpha[e->u.id] = 0.0f;
2672                         } while (e != f->edge);
2673                 }
2674
2675                 for (i = 0; i < ninterior; i++) {
2676                         sys->lambdaPlanar[i] += (float)nlGetVariable(context, 0, i);
2677                         sys->lambdaLength[i] += (float)nlGetVariable(context, 0, ninterior + i);
2678                 }
2679         }
2680
2681         nlDeleteContext(context);
2682
2683         return success;
2684 }
2685
2686 static PBool p_chart_abf_solve(PChart *chart)
2687 {
2688         PVert *v;
2689         PFace *f;
2690         PEdge *e, *e1, *e2, *e3;
2691         PAbfSystem sys;
2692         int i;
2693         float /* lastnorm, */ /* UNUSED */ limit = (chart->nfaces > 100) ? 1.0f : 0.001f;
2694
2695         /* setup id's */
2696         sys.ninterior = sys.nfaces = sys.nangles = 0;
2697
2698         for (v = chart->verts; v; v = v->nextlink) {
2699                 if (p_vert_interior(v)) {
2700                         v->flag |= PVERT_INTERIOR;
2701                         v->u.id = sys.ninterior++;
2702                 }
2703                 else
2704                         v->flag &= ~PVERT_INTERIOR;
2705         }
2706
2707         for (f = chart->faces; f; f = f->nextlink) {
2708                 e1 = f->edge; e2 = e1->next; e3 = e2->next;
2709                 f->u.id = sys.nfaces++;
2710
2711                 /* angle id's are conveniently stored in half edges */
2712                 e1->u.id = sys.nangles++;
2713                 e2->u.id = sys.nangles++;
2714                 e3->u.id = sys.nangles++;
2715         }
2716
2717         p_abf_setup_system(&sys);
2718
2719         /* compute initial angles */
2720         for (f = chart->faces; f; f = f->nextlink) {
2721                 float a1, a2, a3;
2722
2723                 e1 = f->edge; e2 = e1->next; e3 = e2->next;
2724                 p_face_angles(f, &a1, &a2, &a3);
2725
2726                 if (a1 < sys.minangle)
2727                         a1 = sys.minangle;
2728                 else if (a1 > sys.maxangle)
2729                         a1 = sys.maxangle;
2730                 if (a2 < sys.minangle)
2731                         a2 = sys.minangle;
2732                 else if (a2 > sys.maxangle)
2733                         a2 = sys.maxangle;
2734                 if (a3 < sys.minangle)
2735                         a3 = sys.minangle;
2736                 else if (a3 > sys.maxangle)
2737                         a3 = sys.maxangle;
2738
2739                 sys.alpha[e1->u.id] = sys.beta[e1->u.id] = a1;
2740                 sys.alpha[e2->u.id] = sys.beta[e2->u.id] = a2;
2741                 sys.alpha[e3->u.id] = sys.beta[e3->u.id] = a3;
2742
2743                 sys.weight[e1->u.id] = 2.0f / (a1 * a1);
2744                 sys.weight[e2->u.id] = 2.0f / (a2 * a2);
2745                 sys.weight[e3->u.id] = 2.0f / (a3 * a3);
2746         }
2747
2748         for (v = chart->verts; v; v = v->nextlink) {
2749                 if (v->flag & PVERT_INTERIOR) {
2750                         float anglesum = 0.0, scale;
2751
2752                         e = v->edge;
2753                         do {
2754                                 anglesum += sys.beta[e->u.id];
2755                                 e = e->next->next->pair;
2756                         } while (e && (e != v->edge));
2757
2758                         scale = (anglesum == 0.0f) ? 0.0f : 2.0f * (float)M_PI / anglesum;
2759
2760                         e = v->edge;
2761                         do {
2762                                 sys.beta[e->u.id] = sys.alpha[e->u.id] = sys.beta[e->u.id] * scale;
2763                                 e = e->next->next->pair;
2764                         } while (e && (e != v->edge));
2765                 }
2766         }
2767
2768         if (sys.ninterior > 0) {
2769                 p_abf_compute_sines(&sys);
2770
2771                 /* iteration */
2772                 /* lastnorm = 1e10; */ /* UNUSED */
2773
2774                 for (i = 0; i < ABF_MAX_ITER; i++) {
2775                         float norm = p_abf_compute_gradient(&sys, chart);
2776
2777                         /* lastnorm = norm; */ /* UNUSED */
2778
2779                         if (norm < limit)
2780                                 break;
2781
2782                         if (!p_abf_matrix_invert(&sys, chart)) {
2783                                 param_warning("ABF failed to invert matrix");
2784                                 p_abf_free_system(&sys);
2785                                 return P_FALSE;
2786                         }
2787
2788                         p_abf_compute_sines(&sys);
2789                 }
2790
2791                 if (i == ABF_MAX_ITER) {
2792                         param_warning("ABF maximum iterations reached");
2793                         p_abf_free_system(&sys);
2794                         return P_FALSE;
2795                 }
2796         }
2797
2798         chart->u.lscm.abf_alpha = MEM_dupallocN(sys.alpha);
2799         p_abf_free_system(&sys);
2800
2801         return P_TRUE;
2802 }
2803
2804 /* Least Squares Conformal Maps */
2805
2806 static void p_chart_pin_positions(PChart *chart, PVert **pin1, PVert **pin2)
2807 {
2808         if (pin1 == pin2) {
2809                 /* degenerate case */
2810                 PFace *f = chart->faces;
2811                 *pin1 = f->edge->vert;
2812                 *pin2 = f->edge->next->vert;
2813
2814                 (*pin1)->uv[0] = 0.0f;
2815                 (*pin1)->uv[1] = 0.5f;
2816                 (*pin2)->uv[0] = 1.0f;
2817                 (*pin2)->uv[1] = 0.5f;
2818         }
2819         else {
2820                 int diru, dirv, dirx, diry;
2821                 float sub[3];
2822
2823                 sub_v3_v3v3(sub, (*pin1)->co, (*pin2)->co);
2824                 sub[0] = fabsf(sub[0]);
2825                 sub[1] = fabsf(sub[1]);
2826                 sub[2] = fabsf(sub[2]);
2827
2828                 if ((sub[0] > sub[1]) && (sub[0] > sub[2])) {
2829                         dirx = 0;
2830                         diry = (sub[1] > sub[2]) ? 1 : 2;
2831                 }
2832                 else if ((sub[1] > sub[0]) && (sub[1] > sub[2])) {
2833                         dirx = 1;
2834                         diry = (sub[0] > sub[2]) ? 0 : 2;
2835                 }
2836                 else {
2837                         dirx = 2;
2838                         diry = (sub[0] > sub[1]) ? 0 : 1;
2839                 }
2840
2841                 if (dirx == 2) {
2842                         diru = 1;
2843                         dirv = 0;
2844                 }
2845                 else {
2846                         diru = 0;
2847                         dirv = 1;
2848                 }
2849
2850                 (*pin1)->uv[diru] = (*pin1)->co[dirx];
2851                 (*pin1)->uv[dirv] = (*pin1)->co[diry];
2852                 (*pin2)->uv[diru] = (*pin2)->co[dirx];
2853                 (*pin2)->uv[dirv] = (*pin2)->co[diry];
2854         }
2855 }
2856
2857 static PBool p_chart_symmetry_pins(PChart *chart, PEdge *outer, PVert **pin1, PVert **pin2)
2858 {
2859         PEdge *be, *lastbe = NULL, *maxe1 = NULL, *maxe2 = NULL, *be1, *be2;
2860         PEdge *cure = NULL, *firste1 = NULL, *firste2 = NULL, *nextbe;
2861         float maxlen = 0.0f, curlen = 0.0f, totlen = 0.0f, firstlen = 0.0f;
2862         float len1, len2;
2863  
2864         /* find longest series of verts split in the chart itself, these are
2865          * marked during construction */
2866         be = outer;
2867         lastbe = p_boundary_edge_prev(be);
2868         do {
2869                 float len = p_edge_length(be);
2870                 totlen += len;
2871
2872                 nextbe = p_boundary_edge_next(be);
2873
2874                 if ((be->vert->flag & PVERT_SPLIT) ||
2875                     (lastbe->vert->flag & nextbe->vert->flag & PVERT_SPLIT))
2876                 {
2877                         if (!cure) {
2878                                 if (be == outer)
2879                                         firste1 = be;
2880                                 cure = be;
2881                         }
2882                         else
2883                                 curlen += p_edge_length(lastbe);
2884                 }
2885                 else if (cure) {
2886                         if (curlen > maxlen) {
2887                                 maxlen = curlen;
2888                                 maxe1 = cure;
2889                                 maxe2 = lastbe;
2890                         }
2891
2892                         if (firste1 == cure) {
2893                                 firstlen = curlen;
2894                                 firste2 = lastbe;
2895                         }
2896
2897                         curlen = 0.0f;
2898                         cure = NULL;
2899                 }
2900
2901                 lastbe = be;
2902                 be = nextbe;
2903         } while (be != outer);
2904
2905         /* make sure we also count a series of splits over the starting point */
2906         if (cure && (cure != outer)) {
2907                 firstlen += curlen + p_edge_length(be);
2908
2909                 if (firstlen > maxlen) {
2910                         maxlen = firstlen;
2911                         maxe1 = cure;
2912                         maxe2 = firste2;
2913                 }
2914         }
2915
2916         if (!maxe1 || !maxe2 || (maxlen < 0.5f * totlen))
2917                 return P_FALSE;
2918         
2919         /* find pin1 in the split vertices */
2920         be1 = maxe1;
2921         be2 = maxe2;
2922         len1 = 0.0f;
2923         len2 = 0.0f;
2924
2925         do {
2926                 if (len1 < len2) {
2927                         len1 += p_edge_length(be1);
2928                         be1 = p_boundary_edge_next(be1);
2929                 }
2930                 else {
2931                         be2 = p_boundary_edge_prev(be2);
2932                         len2 += p_edge_length(be2);
2933                 }
2934         } while (be1 != be2);
2935
2936         *pin1 = be1->vert;
2937
2938         /* find pin2 outside the split vertices */
2939         be1 = maxe1;
2940         be2 = maxe2;
2941         len1 = 0.0f;
2942         len2 = 0.0f;
2943
2944         do {
2945                 if (len1 < len2) {
2946                         be1 = p_boundary_edge_prev(be1);
2947                         len1 += p_edge_length(be1);
2948                 }
2949                 else {
2950                         len2 += p_edge_length(be2);
2951                         be2 = p_boundary_edge_next(be2);
2952                 }
2953         } while (be1 != be2);
2954
2955         *pin2 = be1->vert;
2956
2957         p_chart_pin_positions(chart, pin1, pin2);
2958
2959         return !equals_v3v3((*pin1)->co, (*pin2)->co);
2960 }
2961
2962 static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **pin2)
2963 {
2964         float minv[3], maxv[3], dirlen;
2965         PVert *v, *minvert[3], *maxvert[3];
2966         int i, dir;
2967
2968         /* find minimum and maximum verts over x/y/z axes */
2969         minv[0] = minv[1] = minv[2] = 1e20;
2970         maxv[0] = maxv[1] = maxv[2] = -1e20;
2971
2972         minvert[0] = minvert[1] = minvert[2] = NULL;
2973         maxvert[0] = maxvert[1] = maxvert[2] = NULL;
2974
2975         for (v = chart->verts; v; v = v->nextlink) {
2976                 for (i = 0; i < 3; i++) {
2977                         if (v->co[i] < minv[i]) {
2978                                 minv[i] = v->co[i];
2979                                 minvert[i] = v;
2980                         }
2981                         if (v->co[i] > maxv[i]) {
2982                                 maxv[i] = v->co[i];
2983                                 maxvert[i] = v;
2984                         }
2985                 }
2986         }
2987
2988         /* find axes with longest distance */
2989         dir = 0;
2990         dirlen = -1.0;
2991
2992         for (i = 0; i < 3; i++) {
2993                 if (maxv[i] - minv[i] > dirlen) {
2994                         dir = i;
2995                         dirlen = maxv[i] - minv[i];
2996                 }
2997         }
2998
2999         *pin1 = minvert[dir];
3000         *pin2 = maxvert[dir];
3001
3002         p_chart_pin_positions(chart, pin1, pin2);
3003 }
3004
3005 static void p_chart_lscm_load_solution(PChart *chart)
3006 {
3007         NLContext *context = chart->u.lscm.context;
3008         PVert *v;
3009
3010         for (v = chart->verts; v; v = v->nextlink) {
3011                 v->uv[0] = nlGetVariable(context, 0, 2 * v->u.id);
3012                 v->uv[1] = nlGetVariable(context, 0, 2 * v->u.id + 1);
3013         }
3014 }
3015
3016 static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
3017 {
3018         PVert *v, *pin1, *pin2;
3019         PBool select = P_FALSE, deselect = P_FALSE;
3020         int npins = 0, id = 0;
3021
3022         /* give vertices matrix indices and count pins */
3023         for (v = chart->verts; v; v = v->nextlink) {
3024                 if (v->flag & PVERT_PIN) {
3025                         npins++;
3026                         if (v->flag & PVERT_SELECT)
3027                                 select = P_TRUE;
3028                 }
3029
3030                 if (!(v->flag & PVERT_SELECT))
3031                         deselect = P_TRUE;
3032         }
3033
3034         if ((live && (!select || !deselect)) || (npins == 1)) {
3035                 chart->u.lscm.context = NULL;
3036         }
3037         else {
3038 #if 0
3039                 p_chart_simplify_compute(chart);
3040                 p_chart_topological_sanity_check(chart);
3041 #endif
3042
3043                 if (abf) {
3044                         if (!p_chart_abf_solve(chart))
3045                                 param_warning("ABF solving failed: falling back to LSCM.\n");
3046                 }
3047
3048                 if (npins <= 1) {
3049                         /* not enough pins, lets find some ourself */
3050                         PEdge *outer;
3051
3052                         p_chart_boundaries(chart, NULL, &outer);
3053
3054                         if (!p_chart_symmetry_pins(chart, outer, &pin1, &pin2))
3055                                 p_chart_extrema_verts(chart, &pin1, &pin2);
3056
3057                         chart->u.lscm.pin1 = pin1;
3058                         chart->u.lscm.pin2 = pin2;
3059                 }
3060                 else {
3061                         chart->flag |= PCHART_NOPACK;
3062                 }
3063
3064                 for (v = chart->verts; v; v = v->nextlink)
3065                         v->u.id = id++;
3066
3067                 chart->u.lscm.context = nlNewContext();
3068                 nlSolverParameteri(chart->u.lscm.context, NL_NB_VARIABLES, 2 * chart->nverts);
3069                 nlSolverParameteri(chart->u.lscm.context, NL_NB_ROWS, 2 * chart->nfaces);
3070                 nlSolverParameteri(chart->u.lscm.context, NL_LEAST_SQUARES, NL_TRUE);
3071         }
3072 }
3073
3074 static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
3075 {
3076         NLContext *context = chart->u.lscm.context;
3077         PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
3078         PFace *f;
3079         float *alpha = chart->u.lscm.abf_alpha;
3080         float area_pinned_up, area_pinned_down;
3081         bool flip_faces;
3082         int row;
3083
3084         nlBegin(context, NL_SYSTEM);
3085
3086 #if 0
3087         /* TODO: make loading pins work for simplify/complexify. */
3088 #endif
3089
3090         for (v = chart->verts; v; v = v->nextlink)
3091                 if (v->flag & PVERT_PIN)
3092                         p_vert_load_pin_select_uvs(handle, v);  /* reload for live */
3093
3094         if (chart->u.lscm.pin1) {
3095                 nlLockVariable(context, 2 * pin1->u.id);
3096                 nlLockVariable(context, 2 * pin1->u.id + 1);
3097                 nlLockVariable(context, 2 * pin2->u.id);
3098                 nlLockVariable(context, 2 * pin2->u.id + 1);
3099
3100                 nlSetVariable(context, 0, 2 * pin1->u.id, pin1->uv[0]);
3101                 nlSetVariable(context, 0, 2 * pin1->u.id + 1, pin1->uv[1]);
3102                 nlSetVariable(context, 0, 2 * pin2->u.id, pin2->uv[0]);
3103                 nlSetVariable(context, 0, 2 * pin2->u.id + 1, pin2->uv[1]);
3104         }
3105         else {
3106                 /* set and lock the pins */
3107                 for (v = chart->verts; v; v = v->nextlink) {
3108                         if (v->flag & PVERT_PIN) {
3109                                 nlLockVariable(context, 2 * v->u.id);
3110                                 nlLockVariable(context, 2 * v->u.id + 1);
3111
3112                                 nlSetVariable(context, 0, 2 * v->u.id, v->uv[0]);
3113                                 nlSetVariable(context, 0, 2 * v->u.id + 1, v->uv[1]);
3114                         }
3115                 }
3116         }
3117
3118         /* detect up direction based on pinned vertices */
3119         area_pinned_up = 0.0f;
3120         area_pinned_down = 0.0f;
3121
3122         for (f = chart->faces; f; f = f->nextlink) {
3123                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
3124                 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
3125
3126                 if ((v1->flag & PVERT_PIN) && (v2->flag & PVERT_PIN) && (v3->flag & PVERT_PIN)) {
3127                         float area = p_face_uv_area_signed(f);
3128
3129                         if (area > 0.0f)
3130                                 area_pinned_up += area;
3131                         else
3132                                 area_pinned_down -= area;
3133                 }
3134         }
3135
3136         flip_faces = (area_pinned_down > area_pinned_up);
3137
3138         /* construct matrix */
3139
3140         nlBegin(context, NL_MATRIX);
3141
3142         row = 0;
3143         for (f = chart->faces; f; f = f->nextlink) {
3144                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
3145                 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
3146                 float a1, a2, a3, ratio, cosine, sine;
3147                 float sina1, sina2, sina3, sinmax;
3148
3149                 if (alpha) {
3150                         /* use abf angles if passed on */
3151                         a1 = *(alpha++);
3152                         a2 = *(alpha++);
3153                         a3 = *(alpha++);
3154                 }
3155                 else
3156                         p_face_angles(f, &a1, &a2, &a3);
3157
3158                 if (flip_faces) {
3159                         SWAP(float, a2, a3);
3160                         SWAP(PEdge *, e2, e3);
3161                         SWAP(PVert *, v2, v3);
3162                 }
3163
3164                 sina1 = sinf(a1);
3165                 sina2 = sinf(a2);
3166                 sina3 = sinf(a3);
3167
3168                 sinmax = max_fff(sina1, sina2, sina3);
3169
3170                 /* shift vertices to find most stable order */
3171                 if (sina3 != sinmax) {
3172                         SHIFT3(PVert *, v1, v2, v3);
3173                         SHIFT3(float, a1, a2, a3);
3174                         SHIFT3(float, sina1, sina2, sina3);
3175
3176                         if (sina2 == sinmax) {
3177                                 SHIFT3(PVert *, v1, v2, v3);
3178                                 SHIFT3(float, a1, a2, a3);
3179                                 SHIFT3(float, sina1, sina2, sina3);
3180                         }
3181                 }
3182
3183                 /* angle based lscm formulation */
3184                 ratio = (sina3 == 0.0f) ? 1.0f : sina2 / sina3;
3185                 cosine = cosf(a1) * ratio;
3186                 sine = sina1 * ratio;
3187
3188                 nlMatrixAdd(context, row, 2 * v1->u.id,   cosine - 1.0f);
3189                 nlMatrixAdd(context, row, 2 * v1->u.id + 1, -sine);
3190                 nlMatrixAdd(context, row, 2 * v2->u.id,   -cosine);
3191                 nlMatrixAdd(context, row, 2 * v2->u.id + 1, sine);
3192                 nlMatrixAdd(context, row, 2 * v3->u.id,   1.0);
3193                 row++;
3194
3195                 nlMatrixAdd(context, row, 2 * v1->u.id,   sine);
3196                 nlMatrixAdd(context, row, 2 * v1->u.id + 1, cosine - 1.0f);
3197                 nlMatrixAdd(context, row, 2 * v2->u.id,   -sine);
3198                 nlMatrixAdd(context, row, 2 * v2->u.id + 1, -cosine);
3199                 nlMatrixAdd(context, row, 2 * v3->u.id + 1, 1.0);
3200                 row++;
3201         }
3202
3203         nlEnd(context, NL_MATRIX);
3204
3205         nlEnd(context, NL_SYSTEM);
3206
3207         if (nlSolve(context, NL_TRUE)) {
3208                 p_chart_lscm_load_solution(chart);
3209                 return P_TRUE;
3210         }
3211         else {
3212                 for (v = chart->verts; v; v = v->nextlink) {
3213                         v->uv[0] = 0.0f;
3214                         v->uv[1] = 0.0f;
3215                 }
3216         }
3217
3218         return P_FALSE;
3219 }
3220
3221 static void p_chart_lscm_end(PChart *chart)
3222 {
3223         if (chart->u.lscm.context)
3224                 nlDeleteContext(chart->u.lscm.context);
3225         
3226         if (chart->u.lscm.abf_alpha) {
3227                 MEM_freeN(chart->u.lscm.abf_alpha);
3228                 chart->u.lscm.abf_alpha = NULL;
3229         }