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