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