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