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