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