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