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