source/blender/src/parametrizer.c

   1
   2 #include "MEM_guardedalloc.h"
   3
   4 #include "BLI_memarena.h"
   5 #include "BLI_arithb.h"
   6 #include "BLI_rand.h"
   7
   8 #include "BKE_utildefines.h"
   9
  10 #include "BIF_editsima.h"
  11 #include "BIF_toolbox.h"
  12
  13 #include "ONL_opennl.h"
  14
  15 #include "parametrizer.h"
  16 #include "parametrizer_intern.h"
  17
  18 #include <math.h>
  19 #include <stdlib.h>
  20 #include <stdio.h>
  21 #include <string.h>
  22
  23 #if defined(_WIN32)
  24 #define M_PI 3.14159265358979323846
  25 #endif
  26
  27 /* Hash */
  28
  29 static int PHashSizes[] = {
  30         1, 3, 5, 11, 17, 37, 67, 131, 257, 521, 1031, 2053, 4099, 8209,
  31         16411, 32771, 65537, 131101, 262147, 524309, 1048583, 2097169,
  32         4194319, 8388617, 16777259, 33554467, 67108879, 134217757, 268435459
  33 };
  34
  35 #define PHASH_hash(ph, item) (((unsigned long) (item))%((unsigned int) (ph)->cursize))
  36
  37 PHash *phash_new(int sizehint)
  38 {
  39         PHash *ph = (PHash*)MEM_callocN(sizeof(PHash), "PHash");
  40         ph->size = 0;
  41         ph->cursize_id = 0;
  42         ph->first = NULL;
  43
  44         while (PHashSizes[ph->cursize_id] < sizehint)
  45                 ph->cursize_id++;
  46
  47         ph->cursize = PHashSizes[ph->cursize_id];
  48         ph->buckets = (PHashLink**)MEM_callocN(ph->cursize*sizeof(*ph->buckets), "PHashBuckets");
  49
  50         return ph;
  51 }
  52
  53 void phash_delete(PHash *ph)
  54 {
  55         MEM_freeN(ph->buckets);
  56         MEM_freeN(ph);
  57 }
  58
  59 void phash_delete_with_links(PHash *ph)
  60 {
  61         PHashLink *link, *next=NULL;
  62
  63         for (link = ph->first; link; link = next) {
  64                 next = link->next;
  65                 MEM_freeN(link);
  66         }
  67
  68         phash_delete(ph);
  69 }
  70
  71 int phash_size(PHash *ph)
  72 {
  73         return ph->size;
  74 }
  75
  76 void phash_insert(PHash *ph, PHashLink *link)
  77 {
  78         int size = ph->cursize;
  79         int hash = PHASH_hash(ph, link->key);
  80         PHashLink *lookup = ph->buckets[hash];
  81
  82         if (lookup == NULL) {
  83                 /* insert in front of the list */
  84                 ph->buckets[hash] = link;
  85                 link->next = ph->first;
  86                 ph->first = link;
  87         }
  88         else {
  89                 /* insert after existing element */
  90                 link->next = lookup->next;
  91                 lookup->next = link;
  92         }
  93
  94         ph->size++;
  95
  96         if (ph->size > (size*3)) {
  97                 PHashLink *next = NULL, *first = ph->first;
  98
  99                 ph->cursize = PHashSizes[++ph->cursize_id];
 100                 MEM_freeN(ph->buckets);
 101                 ph->buckets = (PHashLink**)MEM_callocN(ph->cursize*sizeof(*ph->buckets), "PHashBuckets");
 102                 ph->size = 0;
 103                 ph->first = NULL;
 104
 105                 for (link = first; link; link = next) {
 106                         next = link->next;
 107                         phash_insert(ph, link);
 108                 }
 109         }
 110 }
 111
 112 PHashLink *phash_lookup(PHash *ph, PHashKey key)
 113 {
 114         PHashLink *link;
 115         int hash = PHASH_hash(ph, key);
 116
 117         for (link = ph->buckets[hash]; link; link = link->next)
 118                 if (link->key == key)
 119                         return link;
 120                 else if (PHASH_hash(ph, link->key) != hash)
 121                         return NULL;
 122
 123         return link;
 124 }
 125
 126 PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
 127 {
 128         int hash = PHASH_hash(ph, key);
 129
 130         for (link = link->next; link; link = link->next)
 131                 if (link->key == key)
 132                         return link;
 133                 else if (PHASH_hash(ph, link->key) != hash)
 134                         return NULL;
 135
 136         return link;
 137 }
 138
 139 /* Heap */
 140
 141 #define PHEAP_PARENT(i) ((i-1)>>1)
 142 #define PHEAP_LEFT(i)   ((i<<1)+1)
 143 #define PHEAP_RIGHT(i)  ((i<<1)+2)
 144 #define PHEAP_COMPARE(a, b) (a->value < b->value)
 145 #define PHEAP_EQUALS(a, b) (a->value == b->value)
 146 #define PHEAP_SWAP(heap, i, j) \
 147         { SWAP(int, heap->tree[i]->index, heap->tree[j]->index); \
 148           SWAP(PHeapLink*, heap->tree[i], heap->tree[j]);  }
 149
 150 static void pheap_down(PHeap *heap, int i)
 151 {
 152         while (P_TRUE) {
 153                 int size = heap->size, smallest;
 154                 int l = PHEAP_LEFT(i);
 155                 int r = PHEAP_RIGHT(i);
 156
 157                 smallest = ((l < size) && PHEAP_COMPARE(heap->tree[l], heap->tree[i]))? l: i;
 158
 159                 if ((r < size) && PHEAP_COMPARE(heap->tree[r], heap->tree[smallest]))
 160                         smallest = r;
 161
 162                 if (smallest == i)
 163                         break;
 164
 165                 PHEAP_SWAP(heap, i, smallest);
 166                 i = smallest;
 167         }
 168 }
 169
 170 static void pheap_up(PHeap *heap, int i)
 171 {
 172         while (i > 0) {
 173                 int p = PHEAP_PARENT(i);
 174
 175                 if (PHEAP_COMPARE(heap->tree[p], heap->tree[i]))
 176                         break;
 177
 178                 PHEAP_SWAP(heap, p, i);
 179                 i = p;
 180         }
 181 }
 182
 183 PHeap *pheap_new()
 184 {
 185         /* TODO: replace mallocN with something faster */
 186
 187         PHeap *heap = (PHeap*)MEM_callocN(sizeof(PHeap), "PHeap");
 188         heap->bufsize = 1;
 189         heap->tree = (PHeapLink**)MEM_mallocN(sizeof(PHeapLink*), "PHeapTree");
 190
 191         return heap;
 192 }
 193
 194 void pheap_delete(PHeap *heap)
 195 {
 196         MEM_freeN(heap->tree);
 197         MEM_freeN(heap);
 198 }
 199
 200 PHeapLink *pheap_insert(PHeap *heap, float value, void *ptr)
 201 {
 202         PHeapLink *link;
 203
 204         if ((heap->size + 1) > heap->bufsize) {
 205                 int newsize = heap->bufsize*2;
 206
 207                 PHeapLink **ntree = (PHeapLink**)MEM_mallocN(newsize*sizeof(PHeapLink*), "PHeapTree");
 208                 memcpy(ntree, heap->tree, sizeof(PHeapLink*)*heap->size);
 209                 MEM_freeN(heap->tree);
 210
 211                 heap->tree = ntree;
 212                 heap->bufsize = newsize;
 213         }
 214
 215         param_assert(heap->size < heap->bufsize);
 216
 217         link = MEM_mallocN(sizeof *link, "PHeapLink");
 218         link->value = value;
 219         link->ptr = ptr;
 220         link->index = heap->size;
 221
 222         heap->tree[link->index] = link;
 223
 224         heap->size++;
 225
 226         pheap_up(heap, heap->size-1);
 227
 228         return link;
 229 }
 230
 231 int pheap_empty(PHeap *heap)
 232 {
 233         return (heap->size == 0);
 234 }
 235
 236 int pheap_size(PHeap *heap)
 237 {
 238         return heap->size;
 239 }
 240
 241 void *pheap_min(PHeap *heap)
 242 {
 243         return heap->tree[0]->ptr;
 244 }
 245
 246 void *pheap_popmin(PHeap *heap)
 247 {
 248         void *ptr = heap->tree[0]->ptr;
 249
 250         MEM_freeN(heap->tree[0]);
 251
 252         if (heap->size == 1)
 253                 heap->size--;
 254         else {
 255                 PHEAP_SWAP(heap, 0, heap->size-1);
 256                 heap->size--;
 257
 258                 pheap_down(heap, 0);
 259         }
 260
 261         return ptr;
 262 }
 263
 264 static void pheap_remove(PHeap *heap, PHeapLink *link)
 265 {
 266         int i = link->index;
 267
 268         while (i > 0) {
 269                 int p = PHEAP_PARENT(i);
 270
 271                 PHEAP_SWAP(heap, p, i);
 272                 i = p;
 273         }
 274
 275         pheap_popmin(heap);
 276 }
 277
 278 /* Construction */
 279
 280 PEdge *p_wheel_edge_next(PEdge *e)
 281 {
 282         return e->next->next->pair;
 283 }
 284
 285 PEdge *p_wheel_edge_prev(PEdge *e)
 286 {
 287         return (e->pair)? e->pair->next: NULL;
 288 }
 289
 290 static PVert *p_vert_add(PChart *chart, PHashKey key, float *co, PEdge *e)
 291 {
 292         PVert *v = (PVert*)BLI_memarena_alloc(chart->handle->arena, sizeof *v);
 293         v->co = co;
 294         v->link.key = key;
 295         v->edge = e;
 296
 297         phash_insert(chart->verts, (PHashLink*)v);
 298
 299         return v;
 300 }
 301
 302 static PVert *p_vert_lookup(PChart *chart, PHashKey key, float *co, PEdge *e)
 303 {
 304         PVert *v = (PVert*)phash_lookup(chart->verts, key);
 305
 306         if (v)
 307                 return v;
 308         else
 309                 return p_vert_add(chart, key, co, e);
 310 }
 311
 312 static PVert *p_vert_copy(PChart *chart, PVert *v)
 313 {
 314         PVert *nv = (PVert*)BLI_memarena_alloc(chart->handle->arena, sizeof *nv);
 315         nv->co = v->co;
 316         nv->uv[0] = v->uv[0];
 317         nv->uv[1] = v->uv[1];
 318         nv->link.key = v->link.key;
 319         nv->edge = v->edge;
 320
 321         phash_insert(chart->verts, (PHashLink*)nv);
 322
 323         return nv;
 324 }
 325
 326 static PEdge *p_edge_lookup(PChart *chart, PHashKey *vkeys)
 327 {
 328         PHashKey key = vkeys[0]^vkeys[1];
 329         PEdge *e = (PEdge*)phash_lookup(chart->edges, key);
 330
 331         while (e) {
 332                 if ((e->vert->link.key == vkeys[0]) && (e->next->vert->link.key == vkeys[1]))
 333                         return e;
 334                 else if ((e->vert->link.key == vkeys[1]) && (e->next->vert->link.key == vkeys[0]))
 335                         return e;
 336
 337                 e = (PEdge*)phash_next(chart->edges, key, (PHashLink*)e);
 338         }
 339
 340         return NULL;
 341 }
 342
 343 static void p_face_flip(PFace *f)
 344 {
 345         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 346         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
 347         int f1 = e1->flag, f2 = e2->flag, f3 = e3->flag;
 348
 349         e1->vert = v2;
 350         e1->next = e3;
 351         e1->flag = (f1 & ~PEDGE_VERTEX_FLAGS) | (f2 & PEDGE_VERTEX_FLAGS);
 352
 353         e2->vert = v3;
 354         e2->next = e1;
 355         e2->flag = (f2 & ~PEDGE_VERTEX_FLAGS) | (f3 & PEDGE_VERTEX_FLAGS);
 356
 357         e3->vert = v1;
 358         e3->next = e2;
 359         e3->flag = (f3 & ~PEDGE_VERTEX_FLAGS) | (f1 & PEDGE_VERTEX_FLAGS);
 360 }
 361
 362 static void p_vert_load_pin_select_uvs(PVert *v)
 363 {
 364         PEdge *e;
 365         int nedges = 0;
 366
 367         v->uv[0] = v->uv[1] = 0.0f;
 368         nedges = 0;
 369         e = v->edge;
 370         do {
 371                 if (e->orig_uv && (e->flag & PEDGE_PIN)) {
 372                         if (e->flag & PEDGE_SELECT)
 373                                 v->flag |= PVERT_SELECT;
 374
 375                         v->flag |= PVERT_PIN;
 376                         v->uv[0] += e->orig_uv[0];
 377                         v->uv[1] += e->orig_uv[1];
 378                         nedges++;
 379                 }
 380
 381                 e = p_wheel_edge_next(e);
 382         } while (e && e != (v->edge));
 383
 384         if (nedges > 0) {
 385                 v->uv[0] /= nedges;
 386                 v->uv[1] /= nedges;
 387         }
 388 }
 389
 390 static void p_vert_load_select_uvs(PVert *v)
 391 {
 392         PEdge *e;
 393         int nedges = 0;
 394
 395         v->uv[0] = v->uv[1] = 0.0f;
 396         nedges = 0;
 397         e = v->edge;
 398         do {
 399                 if (e->orig_uv && (e->flag & PEDGE_SELECT))
 400                         v->flag |= PVERT_SELECT;
 401
 402                 v->uv[0] += e->orig_uv[0];
 403                 v->uv[1] += e->orig_uv[1];
 404                 nedges++;
 405
 406                 e = p_wheel_edge_next(e);
 407         } while (e && e != (v->edge));
 408
 409         if (nedges > 0) {
 410                 v->uv[0] /= nedges;
 411                 v->uv[1] /= nedges;
 412         }
 413 }
 414
 415 static void p_extrema_verts(PChart *chart, PVert **v1, PVert **v2)
 416 {
 417         float minv[3], maxv[3], dirlen;
 418         PVert *v, *minvert[3], *maxvert[3];
 419         int i, dir;
 420
 421         /* find minimum and maximum verts over x/y/z axes */
 422         minv[0] = minv[1] = minv[2] = 1e20;
 423         maxv[0] = maxv[1] = maxv[2] = -1e20;
 424
 425         minvert[0] = minvert[1] = minvert[2] = NULL;
 426         maxvert[0] = maxvert[1] = maxvert[2] = NULL;
 427
 428         for (v = (PVert*)chart->verts->first; v; v=v->link.next) {
 429                 for (i = 0; i < 3; i++) {
 430                         if (v->co[i] < minv[i]) {
 431                                 minv[i] = v->co[i];
 432                                 minvert[i] = v;
 433                         }
 434                         if (v->co[i] > maxv[i]) {
 435                                 maxv[i] = v->co[i];
 436                                 maxvert[i] = v;
 437                         }
 438                 }
 439         }
 440
 441         /* find axes with longest distance */
 442         dir = 0;
 443         dirlen = -1.0;
 444
 445         for (i = 0; i < 3; i++) {
 446                 if (maxv[i] - minv[i] > dirlen) {
 447                         dir = i;
 448                         dirlen = maxv[i] - minv[i];
 449                 }
 450         }
 451
 452         if (minvert[dir] == maxvert[dir]) {
 453                 /* degenerate case */
 454                 PFace *f = (PFace*)chart->faces->first;
 455                 *v1 = f->edge->vert;
 456                 *v2 = f->edge->next->vert;
 457
 458                 (*v1)->uv[0] = 0.0f;
 459                 (*v1)->uv[1] = 0.5f;
 460                 (*v2)->uv[0] = 1.0f;
 461                 (*v2)->uv[1] = 0.5f;
 462         }
 463         else {
 464                 *v1 = minvert[dir];
 465                 *v2 = maxvert[dir];
 466
 467                 (*v1)->uv[0] = (*v1)->co[dir];
 468                 (*v1)->uv[1] = (*v1)->co[(dir+1)%3];
 469                 (*v2)->uv[0] = (*v2)->co[dir];
 470                 (*v2)->uv[1] = (*v2)->co[(dir+1)%3];
 471         }
 472 }
 473
 474 static float p_vec_normalise(float *v)
 475 {
 476    float d;
 477
 478     d = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
 479
 480     if(d != 0.0f) {
 481                 d = 1.0f/d;
 482
 483                 v[0] *= d;
 484                 v[1] *= d;
 485                 v[2] *= d;
 486         }
 487
 488         return d;
 489 }
 490
 491 static float p_vec_angle_cos(float *v1, float *v2, float *v3)
 492 {
 493         float d1[3], d2[3];
 494
 495         d1[0] = v1[0] - v2[0];
 496         d1[1] = v1[1] - v2[1];
 497         d1[2] = v1[2] - v2[2];
 498
 499         d2[0] = v3[0] - v2[0];
 500         d2[1] = v3[1] - v2[1];
 501         d2[2] = v3[2] - v2[2];
 502
 503         p_vec_normalise(d1);
 504         p_vec_normalise(d2);
 505
 506         return d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2];
 507 }
 508
 509 static float p_vec_angle(float *v1, float *v2, float *v3)
 510 {
 511         float dot = p_vec_angle_cos(v1, v2, v3);
 512
 513         if (dot <= -1.0f)
 514                 return (float)M_PI;
 515         else if (dot >= 1.0f)
 516                 return 0.0f;
 517         else
 518                 return (float)acos(dot);
 519 }
 520
 521 static void p_face_angles(PFace *f, float *a1, float *a2, float *a3)
 522 {
 523         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 524         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
 525
 526         *a1 = p_vec_angle(v3->co, v1->co, v2->co);
 527         *a2 = p_vec_angle(v1->co, v2->co, v3->co);
 528         *a3 = M_PI - *a2 - *a1;
 529 }
 530
 531 static float p_face_area(PFace *f)
 532 {
 533         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 534         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
 535
 536         return AreaT3Dfl(v1->co, v2->co, v3->co);
 537 }
 538
 539 static float p_face_uv_area_signed(PFace *f)
 540 {
 541         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 542         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
 543
 544         return 0.5f*(((v2->uv[0]-v1->uv[0]) * (v3->uv[1]-v1->uv[1])) -
 545                     ((v3->uv[0]-v1->uv[0]) * (v2->uv[1]-v1->uv[1])));
 546 }
 547
 548 static float p_face_uv_area(PFace *f)
 549 {
 550         return fabs(p_face_uv_area_signed(f));
 551 }
 552
 553 static void p_chart_area(PChart *chart, float *uv_area, float *area)
 554 {
 555         PFace *f;
 556
 557         *uv_area = *area = 0.0f;
 558
 559         for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
 560                 *uv_area += p_face_uv_area(f);
 561                 *area += p_face_area(f);
 562         }
 563 }
 564
 565 static PChart *p_chart_new(PHandle *handle)
 566 {
 567         PChart *chart = (PChart*)MEM_callocN(sizeof*chart, "PChart");
 568         chart->verts = phash_new(1);
 569         chart->edges = phash_new(1);
 570         chart->faces = phash_new(1);
 571         chart->handle = handle;
 572
 573         return chart;
 574 }
 575
 576 static void p_chart_delete(PChart *chart)
 577 {
 578         /* the actual links are free by memarena */
 579         phash_delete(chart->verts);
 580         phash_delete(chart->edges);
 581         phash_delete(chart->faces);
 582
 583         MEM_freeN(chart);
 584 }
 585
 586 static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
 587 {
 588         float *uv1, *uv2, *uvp1, *uvp2;
 589         float limit[2];
 590
 591         uv1 = e->orig_uv;
 592         uv2 = e->next->orig_uv;
 593
 594         if (e->vert->link.key == ep->vert->link.key) {
 595                 uvp1 = ep->orig_uv;
 596                 uvp2 = ep->next->orig_uv;
 597         }
 598         else {
 599                 uvp1 = ep->next->orig_uv;
 600                 uvp2 = ep->orig_uv;
 601         }
 602
 603         get_connected_limit_tface_uv(limit);
 604
 605         if((fabs(uv1[0]-uvp1[0]) > limit[0]) && (fabs(uv1[1]-uvp1[1]) > limit[1])) {
 606                 e->flag |= PEDGE_SEAM;
 607                 ep->flag |= PEDGE_SEAM;
 608                 return P_TRUE;
 609         }
 610         if((fabs(uv2[0]-uvp2[0]) > limit[0]) && (fabs(uv2[1]-uvp2[1]) > limit[1])) {
 611                 e->flag |= PEDGE_SEAM;
 612                 ep->flag |= PEDGE_SEAM;
 613                 return P_TRUE;
 614         }
 615
 616         return P_FALSE;
 617 }
 618
 619 static PBool p_edge_has_pair(PChart *chart, PEdge *e, PEdge **pair, PBool impl)
 620 {
 621         PHashKey key;
 622         PEdge *pe;
 623         PVert *v1, *v2;
 624         PHashKey key1 = e->vert->link.key;
 625         PHashKey key2 = e->next->vert->link.key;
 626
 627         if (e->flag & PEDGE_SEAM)
 628                 return P_FALSE;
 629
 630         key = key1 ^ key2;
 631         pe = (PEdge*)phash_lookup(chart->edges, key);
 632         *pair = NULL;
 633
 634         while (pe) {
 635                 if (pe != e) {
 636                         v1 = pe->vert;
 637                         v2 = pe->next->vert;
 638
 639                         if (((v1->link.key == key1) && (v2->link.key == key2)) ||
 640                                 ((v1->link.key == key2) && (v2->link.key == key1))) {
 641
 642                                 /* don't connect seams and t-junctions */
 643                                 if ((pe->flag & PEDGE_SEAM) || *pair ||
 644                                     (impl && p_edge_implicit_seam(e, pe))) {
 645                                         *pair = NULL;
 646                                         return P_FALSE;
 647                                 }
 648
 649                                 *pair = pe;
 650                         }
 651                 }
 652
 653                 pe = (PEdge*)phash_next(chart->edges, key, (PHashLink*)pe);
 654         }
 655
 656         if (*pair && (e->vert == (*pair)->vert)) {
 657                 if ((*pair)->next->pair || (*pair)->next->next->pair) {
 658                         /* non unfoldable, maybe mobius ring or klein bottle */
 659                         *pair = NULL;
 660                         return P_FALSE;
 661                 }
 662         }
 663
 664         return (*pair != NULL);
 665 }
 666
 667 static PBool p_edge_connect_pair(PChart *chart, PEdge *e, PEdge ***stack, PBool impl)
 668 {
 669         PEdge *pair;
 670
 671         if(!e->pair && p_edge_has_pair(chart, e, &pair, impl)) {
 672                 if (e->vert == pair->vert)
 673                         p_face_flip(pair->face);
 674
 675                 e->pair = pair;
 676                 pair->pair = e;
 677
 678                 if (!(pair->face->flag & PFACE_CONNECTED)) {
 679                         **stack = pair;
 680                         (*stack)++;
 681                 }
 682         }
 683
 684         return (e->pair != NULL);
 685 }
 686
 687 static int p_connect_pairs(PChart *chart, PBool impl)
 688 {
 689         PEdge **stackbase = MEM_mallocN(sizeof*stackbase * phash_size(chart->faces), "Pstackbase");
 690         PEdge **stack = stackbase;
 691         PFace *f, *first;
 692         PEdge *e, *e1, *e2;
 693         int ncharts = 0;
 694
 695         /* connect pairs, count edges, set vertex-edge pointer to a pairless edge */
 696         for (first=(PFace*)chart->faces->first; first; first=first->link.next) {
 697                 if (first->flag & PFACE_CONNECTED)
 698                         continue;
 699
 700                 *stack = first->edge;
 701                 stack++;
 702
 703                 while (stack != stackbase) {
 704                         stack--;
 705                         e = *stack;
 706                         e1 = e->next;
 707                         e2 = e1->next;
 708
 709                         f = e->face;
 710                         f->flag |= PFACE_CONNECTED;
 711
 712                         /* assign verts to charts so we can sort them later */
 713                         f->u.chart = ncharts;
 714
 715                         if (!p_edge_connect_pair(chart, e, &stack, impl))
 716                                 e->vert->edge = e;
 717                         if (!p_edge_connect_pair(chart, e1, &stack, impl))
 718                                 e1->vert->edge = e1;
 719                         if (!p_edge_connect_pair(chart, e2, &stack, impl))
 720                                 e2->vert->edge = e2;
 721                 }
 722
 723                 ncharts++;
 724         }
 725
 726         MEM_freeN(stackbase);
 727
 728         return ncharts;
 729 }
 730
 731 static void p_split_vert(PChart *chart, PEdge *e)
 732 {
 733         PEdge *we, *lastwe = NULL;
 734         PVert *v = e->vert;
 735         PBool copy = P_TRUE;
 736
 737         if (e->flag & PEDGE_VERTEX_SPLIT)
 738                 return;
 739
 740         /* rewind to start */
 741         lastwe = e;
 742         for (we = p_wheel_edge_prev(e); we && (we != e); we = p_wheel_edge_prev(we))
 743                 lastwe = we;
 744
 745         /* go over all edges in wheel */
 746         for (we = lastwe; we; we = p_wheel_edge_next(we)) {
 747                 if (we->flag & PEDGE_VERTEX_SPLIT)
 748                         break;
 749
 750                 we->flag |= PEDGE_VERTEX_SPLIT;
 751
 752                 if (we == v->edge) {
 753                         /* found it, no need to copy */
 754                         copy = P_FALSE;
 755                         phash_insert(chart->verts, (PHashLink*)v);
 756                 }
 757         }
 758
 759         if (copy) {
 760                 /* not found, copying */
 761                 v = p_vert_copy(chart, v);
 762                 v->edge = lastwe;
 763
 764                 we = lastwe;
 765                 do {
 766                         we->vert = v;
 767                         we = p_wheel_edge_next(we);
 768                 } while (we && (we != lastwe));
 769         }
 770 }
 771
 772 static PChart **p_split_charts(PHandle *handle, PChart *chart, int ncharts)
 773 {
 774         PChart **charts = MEM_mallocN(sizeof*charts * ncharts, "PCharts"), *nchart;
 775         PFace *f, *nextf;
 776         int i;
 777
 778         for (i = 0; i < ncharts; i++)
 779                 charts[i] = p_chart_new(handle);
 780
 781         f = (PFace*)chart->faces->first;
 782         while (f) {
 783                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 784                 nextf = f->link.next;
 785
 786                 nchart = charts[f->u.chart];
 787
 788                 phash_insert(nchart->faces, (PHashLink*)f);
 789                 phash_insert(nchart->edges, (PHashLink*)e1);
 790                 phash_insert(nchart->edges, (PHashLink*)e2);
 791                 phash_insert(nchart->edges, (PHashLink*)e3);
 792
 793                 p_split_vert(nchart, e1);
 794                 p_split_vert(nchart, e2);
 795                 p_split_vert(nchart, e3);
 796
 797                 f = nextf;
 798         }
 799
 800         return charts;
 801 }
 802
 803 static void p_face_backup_uvs(PFace *f)
 804 {
 805         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 806
 807         e1->old_uv[0] = e1->orig_uv[0];
 808         e1->old_uv[1] = e1->orig_uv[1];
 809         e2->old_uv[0] = e2->orig_uv[0];
 810         e2->old_uv[1] = e2->orig_uv[1];
 811         e3->old_uv[0] = e3->orig_uv[0];
 812         e3->old_uv[1] = e3->orig_uv[1];
 813 }
 814
 815 static void p_face_restore_uvs(PFace *f)
 816 {
 817         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 818
 819         e1->orig_uv[0] = e1->old_uv[0];
 820         e1->orig_uv[1] = e1->old_uv[1];
 821         e2->orig_uv[0] = e2->old_uv[0];
 822         e2->orig_uv[1] = e2->old_uv[1];
 823         e3->orig_uv[0] = e3->old_uv[0];
 824         e3->orig_uv[1] = e3->old_uv[1];
 825 }
 826
 827 static PFace *p_face_add(PChart *chart, ParamKey key, ParamKey *vkeys,
 828                          float *co[3], float *uv[3], int i1, int i2, int i3,
 829                          ParamBool *pin, ParamBool *select)
 830 {
 831         PFace *f;
 832         PEdge *e1, *e2, *e3;
 833
 834         /* allocate */
 835         f = (PFace*)BLI_memarena_alloc(chart->handle->arena, sizeof *f);
 836
 837         e1 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e1);
 838         e2 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e2);
 839         e3 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e3);
 840
 841         /* set up edges */
 842         f->edge = e1;
 843         e1->face = e2->face = e3->face = f;
 844
 845         e1->next = e2;
 846         e2->next = e3;
 847         e3->next = e1;
 848
 849         if (co && uv) {
 850                 e1->vert = p_vert_lookup(chart, vkeys[i1], co[i1], e1);
 851                 e2->vert = p_vert_lookup(chart, vkeys[i2], co[i2], e2);
 852                 e3->vert = p_vert_lookup(chart, vkeys[i3], co[i3], e3);
 853
 854                 e1->orig_uv = uv[i1];
 855                 e2->orig_uv = uv[i2];
 856                 e3->orig_uv = uv[i3];
 857         }
 858         else {
 859                 /* internal call to add face */
 860                 e1->vert = e2->vert = e3->vert = NULL;
 861                 e1->orig_uv = e2->orig_uv = e3->orig_uv = NULL;
 862         }
 863
 864         if (pin) {
 865                 if (pin[i1]) e1->flag |= PEDGE_PIN;
 866                 if (pin[i2]) e2->flag |= PEDGE_PIN;
 867                 if (pin[i3]) e3->flag |= PEDGE_PIN;
 868         }
 869
 870         if (select) {
 871                 if (select[i1]) e1->flag |= PEDGE_SELECT;
 872                 if (select[i2]) e2->flag |= PEDGE_SELECT;
 873                 if (select[i3]) e3->flag |= PEDGE_SELECT;
 874         }
 875
 876         /* insert into hash */
 877         f->link.key = key;
 878         phash_insert(chart->faces, (PHashLink*)f);
 879
 880         e1->link.key = vkeys[i1]^vkeys[i2];
 881         e2->link.key = vkeys[i2]^vkeys[i3];
 882         e3->link.key = vkeys[i3]^vkeys[i1];
 883
 884         phash_insert(chart->edges, (PHashLink*)e1);
 885         phash_insert(chart->edges, (PHashLink*)e2);
 886         phash_insert(chart->edges, (PHashLink*)e3);
 887
 888         return f;
 889 }
 890
 891 static PBool p_quad_split_direction(float **co)
 892 {
 893     float a1, a2;
 894
 895         a1 = p_vec_angle_cos(co[0], co[1], co[2]);
 896         a1 += p_vec_angle_cos(co[1], co[0], co[2]);
 897         a1 += p_vec_angle_cos(co[2], co[0], co[1]);
 898
 899         a2 = p_vec_angle_cos(co[0], co[1], co[3]);
 900         a2 += p_vec_angle_cos(co[1], co[0], co[3]);
 901         a2 += p_vec_angle_cos(co[3], co[0], co[1]);
 902
 903         return (a1 > a2);
 904 }
 905
 906 static float p_edge_length(PEdge *e)
 907 {
 908     PVert *v1 = e->vert, *v2 = e->next->vert;
 909     float d[3];
 910
 911     d[0] = v2->co[0] - v1->co[0];
 912     d[1] = v2->co[1] - v1->co[1];
 913     d[2] = v2->co[2] - v1->co[2];
 914
 915     return sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
 916 }
 917
 918 static float p_edge_uv_length(PEdge *e)
 919 {
 920     PVert *v1 = e->vert, *v2 = e->next->vert;
 921     float d[3];
 922
 923     d[0] = v2->uv[0] - v1->uv[0];
 924     d[1] = v2->uv[1] - v1->uv[1];
 925
 926     return sqrt(d[0]*d[0] + d[1]*d[1]);
 927 }
 928
 929 void p_chart_uv_bbox(PChart *chart, float *minv, float *maxv)
 930 {
 931         PVert *v;
 932
 933         INIT_MINMAX2(minv, maxv);
 934
 935         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
 936                 DO_MINMAX2(v->uv, minv, maxv);
 937         }
 938 }
 939
 940 static void p_chart_uv_scale(PChart *chart, float scale)
 941 {
 942         PVert *v;
 943
 944         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
 945                 v->uv[0] *= scale;
 946                 v->uv[1] *= scale;
 947         }
 948 }
 949
 950 static void p_chart_uv_translate(PChart *chart, float trans[2])
 951 {
 952         PVert *v;
 953
 954         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
 955                 v->uv[0] += trans[0];
 956                 v->uv[1] += trans[1];
 957         }
 958 }
 959
 960 static void p_chart_boundaries(PChart *chart, int *nboundaries, PEdge **outer)
 961 {
 962     PEdge *e, *be;
 963     float len, maxlen = -1.0;
 964
 965         *nboundaries = 0;
 966         *outer = NULL;
 967
 968         for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
 969         if (e->pair || (e->flag & PEDGE_DONE))
 970             continue;
 971
 972                 (*nboundaries)++;
 973         len = 0.0f;
 974
 975                 be = e;
 976                 do {
 977             be->flag |= PEDGE_DONE;
 978             len += p_edge_length(be);
 979                         be = be->next->vert->edge;
 980         } while(be != e);
 981
 982         if (len > maxlen) {
 983                         *outer = e;
 984             maxlen = len;
 985         }
 986     }
 987
 988         for (e=(PEdge*)chart->edges->first; e; e=e->link.next)
 989         e->flag &= ~PEDGE_DONE;
 990 }
 991
 992 static float p_edge_boundary_angle(PEdge *e)
 993 {
 994         PEdge *we;
 995         PVert *v, *v1, *v2;
 996         float angle;
 997         int n = 0;
 998
 999         v = e->vert;
1000
1001         /* concave angle check -- could be better */
1002         angle = M_PI;
1003
1004         we = v->edge;
1005         do {
1006                 v1 = we->next->vert;
1007                 v2 = we->next->next->vert;
1008                 angle -= p_vec_angle(v1->co, v->co, v2->co);
1009
1010                 we = we->next->next->pair;
1011                 n++;
1012         } while (we && (we != v->edge));
1013
1014         return angle;
1015 }
1016
1017 static PEdge *p_boundary_edge_next(PEdge *e)
1018 {
1019         return e->next->vert->edge;
1020 }
1021
1022 static PEdge *p_boundary_edge_prev(PEdge *e)
1023 {
1024     PEdge *we = e, *last;
1025
1026         do {
1027                 last = we;
1028                 we = p_wheel_edge_next(we);
1029         } while (we && (we != e));
1030
1031         return last->next->next;
1032 }
1033
1034 static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
1035 {
1036         PEdge *e, *e1, *e2;
1037         PHashKey vkeys[3];
1038         PFace *f;
1039         struct PHeap *heap = pheap_new(nedges);
1040         float angle;
1041
1042         e = be;
1043         do {
1044                 angle = p_edge_boundary_angle(e);
1045                 e->u.heaplink = pheap_insert(heap, angle, e);
1046
1047                 e = e->next->vert->edge;
1048         } while(e != be);
1049
1050         if (nedges == 2) {
1051                 /* no real boundary, but an isolated seam */
1052                 e = be->next->vert->edge;
1053                 e->pair = be;
1054                 be->pair = e;
1055
1056                 pheap_remove(heap, e->u.heaplink);
1057                 pheap_remove(heap, be->u.heaplink);
1058         }
1059         else {
1060                 while (nedges > 2) {
1061                         PEdge *ne, *ne1, *ne2;
1062
1063                         e = pheap_popmin(heap);
1064
1065                         e1 = p_boundary_edge_prev(e);
1066                         e2 = p_boundary_edge_next(e);
1067
1068                         pheap_remove(heap, e1->u.heaplink);
1069                         pheap_remove(heap, e2->u.heaplink);
1070                         e->u.heaplink = e1->u.heaplink = e2->u.heaplink = NULL;
1071
1072                         e->flag |= PEDGE_FILLED;
1073                         e1->flag |= PEDGE_FILLED;
1074
1075                         vkeys[0] = e->vert->link.key;
1076                         vkeys[1] = e1->vert->link.key;
1077                         vkeys[2] = e2->vert->link.key;
1078
1079                         f = p_face_add(chart, -1, vkeys, NULL, NULL, 0, 1, 2, NULL, NULL);
1080                         f->flag |= PFACE_FILLED;
1081
1082                         ne = f->edge->next->next;
1083                         ne1 = f->edge;
1084                         ne2 = f->edge->next;
1085
1086                         ne->flag = ne1->flag = ne2->flag = PEDGE_FILLED;
1087
1088                         e->pair = ne;
1089                         ne->pair = e;
1090                         e1->pair = ne1;
1091                         ne1->pair = e1;
1092
1093                         ne->vert = e2->vert;
1094                         ne1->vert = e->vert;
1095                         ne2->vert = e1->vert;
1096
1097                         if (nedges == 3) {
1098                                 e2->pair = ne2;
1099                                 ne2->pair = e2;
1100                         }
1101                         else {
1102                                 ne2->vert->edge = ne2;
1103
1104                                 ne2->u.heaplink = pheap_insert(heap, p_edge_boundary_angle(ne2), ne2);
1105                                 e2->u.heaplink = pheap_insert(heap, p_edge_boundary_angle(e2), e2);
1106                         }
1107
1108                         nedges--;
1109                 }
1110         }
1111
1112         pheap_delete(heap);
1113 }
1114
1115 static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
1116 {
1117         PEdge *e, *enext, *be;
1118         int nedges;
1119
1120         for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
1121                 enext = e->link.next;
1122
1123         if (e->pair || (e->flag & PEDGE_FILLED))
1124             continue;
1125
1126                 nedges = 0;
1127                 be = e;
1128                 do {
1129                         be->flag |= PEDGE_FILLED;
1130                         be = be->next->vert->edge;
1131                         nedges++;
1132                 } while(be != e);
1133
1134                 if (e != outer)
1135                         p_chart_fill_boundary(chart, e, nedges);
1136     }
1137 }
1138
1139 static void p_flush_uvs(PChart *chart)
1140 {
1141         PEdge *e;
1142
1143         for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
1144                 if (e->orig_uv) {
1145                         e->orig_uv[0] = e->vert->uv[0];
1146                         e->orig_uv[1] = e->vert->uv[1];
1147                 }
1148         }
1149 }
1150
1151 static void p_flush_uvs_blend(PChart *chart, float blend)
1152 {
1153         PEdge *e;
1154         float invblend = 1.0f - blend;
1155
1156         for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
1157                 if (e->orig_uv) {
1158                         e->orig_uv[0] = blend*e->old_uv[0] + invblend*e->vert->uv[0];
1159                         e->orig_uv[1] = blend*e->old_uv[1] + invblend*e->vert->uv[1];
1160                 }
1161         }
1162 }
1163
1164 /* Exported */
1165
1166 ParamHandle *param_construct_begin()
1167 {
1168         PHandle *handle = MEM_callocN(sizeof*handle, "PHandle");
1169         handle->construction_chart = p_chart_new(handle);
1170         handle->state = PHANDLE_STATE_ALLOCATED;
1171         handle->arena = BLI_memarena_new((1<<16));
1172
1173         return (ParamHandle*)handle;
1174 }
1175
1176 void param_delete(ParamHandle *handle)
1177 {
1178         PHandle *phandle = (PHandle*)handle;
1179         int i;
1180
1181         param_assert((phandle->state == PHANDLE_STATE_ALLOCATED) ||
1182                      (phandle->state == PHANDLE_STATE_CONSTRUCTED));
1183
1184         for (i = 0; i < phandle->ncharts; i++)
1185                 p_chart_delete(phandle->charts[i]);
1186
1187         if (phandle->charts)
1188                 MEM_freeN(phandle->charts);
1189
1190         if (phandle->construction_chart)
1191                 p_chart_delete(phandle->construction_chart);
1192
1193         BLI_memarena_free(phandle->arena);
1194         MEM_freeN(phandle);
1195 }
1196
1197 void param_face_add(ParamHandle *handle, ParamKey key, int nverts,
1198                     ParamKey *vkeys, float **co, float **uv,
1199                     ParamBool *pin, ParamBool *select)
1200 {
1201         PHandle *phandle = (PHandle*)handle;
1202         PChart *chart = phandle->construction_chart;
1203
1204         param_assert(phash_lookup(chart->faces, key) == NULL);
1205         param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
1206         param_assert((nverts == 3) || (nverts == 4));
1207
1208         if (nverts == 4) {
1209                 if (!p_quad_split_direction(co)) {
1210                         p_face_add(chart, key, vkeys, co, uv, 0, 1, 2, pin, select);
1211                         p_face_add(chart, key, vkeys, co, uv, 0, 2, 3, pin, select);
1212                 }
1213                 else {
1214                         p_face_add(chart, key, vkeys, co, uv, 0, 1, 3, pin, select);
1215                         p_face_add(chart, key, vkeys, co, uv, 1, 2, 3, pin, select);
1216                 }
1217         }
1218         else
1219                 p_face_add(chart, key, vkeys, co, uv, 0, 1, 2, pin, select);
1220 }
1221
1222 void param_edge_set_seam(ParamHandle *handle, ParamKey *vkeys)
1223 {
1224         PHandle *phandle = (PHandle*)handle;
1225         PChart *chart = phandle->construction_chart;
1226         PEdge *e;
1227
1228         param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
1229
1230         e = p_edge_lookup(chart, vkeys);
1231         if (e)
1232                 e->flag |= PEDGE_SEAM;
1233 }
1234
1235 void param_construct_end(ParamHandle *handle, ParamBool fill, ParamBool impl)
1236 {
1237         PHandle *phandle = (PHandle*)handle;
1238         PChart *chart = phandle->construction_chart;
1239         int i, j, nboundaries = 0;
1240         PEdge *outer;
1241
1242         param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
1243
1244         phandle->ncharts = p_connect_pairs(chart, impl);
1245         phandle->charts = p_split_charts(phandle, chart, phandle->ncharts);
1246
1247         p_chart_delete(chart);
1248         phandle->construction_chart = NULL;
1249
1250         for (i = j = 0; i < phandle->ncharts; i++) {
1251                 p_chart_boundaries(phandle->charts[i], &nboundaries, &outer);
1252
1253                 if (nboundaries == 0) {
1254                         p_chart_delete(phandle->charts[i]);
1255                         continue;
1256                 }
1257
1258                 phandle->charts[j] = phandle->charts[i];
1259                 j++;
1260
1261 #if 0
1262                 if (fill && (nboundaries > 1))
1263                         p_chart_fill_boundaries(phandle->charts[i], outer);
1264 #endif
1265         }
1266
1267         phandle->ncharts = j;
1268
1269         phandle->state = PHANDLE_STATE_CONSTRUCTED;
1270 }
1271
1272 /* Least Squares Conformal Maps */
1273
1274 static void p_chart_lscm_load_solution(PChart *chart)
1275 {
1276         PVert *v;
1277
1278         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
1279                 v->uv[0] = nlGetVariable(2*v->u.index);
1280                 v->uv[1] = nlGetVariable(2*v->u.index + 1);
1281         }
1282 }
1283
1284 static void p_chart_lscm_begin(PChart *chart, PBool live)
1285 {
1286         PVert *v, *pin1, *pin2;
1287         PBool select = P_FALSE;
1288         int npins = 0, id = 0;
1289
1290         /* give vertices matrix indices and count pins */
1291         for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
1292                 p_vert_load_pin_select_uvs(v);
1293
1294                 if (v->flag & PVERT_PIN)
1295                         npins++;
1296
1297                 if (v->flag & PVERT_SELECT)
1298                         select = P_TRUE;
1299
1300                 v->u.index = id++;
1301         }
1302
1303         if ((live && !select) || (npins == 1)) {
1304                 chart->u.lscm.context = NULL;
1305         }
1306         else {
1307                 if (npins <= 1) {
1308                         /* not enough pins, lets find some ourself */
1309                         p_extrema_verts(chart, &pin1, &pin2);
1310
1311                         chart->u.lscm.pin1 = pin1;
1312                         chart->u.lscm.pin2 = pin2;
1313                 }
1314                 else {
1315                         chart->flag |= PCHART_NOPACK;
1316                 }
1317
1318                 nlNewContext();
1319                 nlSolverParameteri(NL_NB_VARIABLES, 2*phash_size(chart->verts));
1320                 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
1321
1322                 chart->u.lscm.context = nlGetCurrent();
1323         }
1324 }
1325
1326 static PBool p_chart_lscm_solve(PChart *chart)
1327 {
1328         PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
1329         PFace *f;
1330
1331         nlMakeCurrent(chart->u.lscm.context);
1332
1333         nlBegin(NL_SYSTEM);
1334
1335         for (v=(PVert*)chart->verts->first; v; v=v->link.next)
1336                 if (v->flag & PVERT_PIN)
1337                         p_vert_load_pin_select_uvs(v);
1338
1339         if (chart->u.lscm.pin1) {
1340                 nlLockVariable(2*pin1->u.index);
1341                 nlLockVariable(2*pin1->u.index + 1);
1342                 nlLockVariable(2*pin2->u.index);
1343                 nlLockVariable(2*pin2->u.index + 1);
1344
1345                 nlSetVariable(2*pin1->u.index, pin1->uv[0]);
1346                 nlSetVariable(2*pin1->u.index + 1, pin1->uv[1]);
1347                 nlSetVariable(2*pin2->u.index, pin2->uv[0]);
1348                 nlSetVariable(2*pin2->u.index + 1, pin2->uv[1]);
1349         }
1350         else {
1351                 /* set and lock the pins */
1352                 for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
1353                         if (v->flag & PVERT_PIN) {
1354                                 nlLockVariable(2*v->u.index);
1355                                 nlLockVariable(2*v->u.index + 1);
1356
1357                                 nlSetVariable(2*v->u.index, v->uv[0]);
1358                                 nlSetVariable(2*v->u.index + 1, v->uv[1]);
1359                         }
1360                 }
1361         }
1362
1363         /* construct matrix */
1364
1365         nlBegin(NL_MATRIX);
1366
1367         for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
1368                 PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1369                 PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
1370                 float a1, a2, a3, ratio, cosine, sine;
1371                 float sina1, sina2, sina3, sinmax;
1372
1373                 if (chart->u.lscm.abf_alpha) {
1374                         /* use abf angles if passed on */
1375                         a1 = *(chart->u.lscm.abf_alpha++);
1376                         a2 = *(chart->u.lscm.abf_alpha++);
1377                         a3 = *(chart->u.lscm.abf_alpha++);
1378                 }
1379                 else
1380                         p_face_angles(f, &a1, &a2, &a3);
1381
1382                 sina1 = sin(a1);
1383                 sina2 = sin(a2);
1384                 sina3 = sin(a3);
1385
1386                 sinmax = MAX3(sina1, sina2, sina3);
1387
1388                 /* shift vertices to find most stable order */
1389                 #define SHIFT3(type, a, b, c) \
1390                         { type tmp; tmp = a; a = c; c = b; b = tmp; }
1391
1392                 if (sina3 != sinmax) {
1393                         SHIFT3(PVert*, v1, v2, v3);
1394                         SHIFT3(float, a1, a2, a3);
1395                         SHIFT3(float, sina1, sina2, sina3);
1396
1397                         if (sina2 == sinmax) {
1398                                 SHIFT3(PVert*, v1, v2, v3);
1399                                 SHIFT3(float, a1, a2, a3);
1400                                 SHIFT3(float, sina1, sina2, sina3);
1401                         }
1402                 }
1403
1404                 /* angle based lscm formulation */
1405                 ratio = (sina3 == 0.0f)? 0.0f: sina2/sina3;
1406                 cosine = cos(a1)*ratio;
1407                 sine = sina1*ratio;
1408
1409                 nlBegin(NL_ROW);
1410                 nlCoefficient(2*v1->u.index,   cosine - 1.0);
1411                 nlCoefficient(2*v1->u.index+1, -sine);
1412                 nlCoefficient(2*v2->u.index,   -cosine);
1413                 nlCoefficient(2*v2->u.index+1, sine);
1414                 nlCoefficient(2*v3->u.index,   1.0);
1415                 nlEnd(NL_ROW);
1416
1417                 nlBegin(NL_ROW);
1418                 nlCoefficient(2*v1->u.index,   sine);
1419                 nlCoefficient(2*v1->u.index+1, cosine - 1.0);
1420                 nlCoefficient(2*v2->u.index,   -sine);
1421                 nlCoefficient(2*v2->u.index+1, -cosine);
1422                 nlCoefficient(2*v3->u.index+1, 1.0);
1423                 nlEnd(NL_ROW);
1424         }
1425
1426         nlEnd(NL_MATRIX);
1427
1428         nlEnd(NL_SYSTEM);
1429
1430         if (nlSolveAdvanced(NULL, NL_TRUE)) {
1431                 p_chart_lscm_load_solution(chart);
1432                 return P_TRUE;
1433         }
1434
1435         return P_FALSE;
1436 }
1437
1438 static void p_chart_lscm_end(PChart *chart)
1439 {
1440         if (chart->u.lscm.context)
1441                 nlDeleteContext(chart->u.lscm.context);
1442
1443         chart->u.lscm.context = NULL;
1444         chart->u.lscm.pin1 = NULL;
1445         chart->u.lscm.pin2 = NULL;
1446 }
1447
1448 void param_lscm_begin(ParamHandle *handle, ParamBool live)
1449 {
1450         PHandle *phandle = (PHandle*)handle;
1451         int i;
1452
1453         param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
1454         phandle->state = PHANDLE_STATE_LSCM;
1455
1456         for (i = 0; i < phandle->ncharts; i++)
1457                 p_chart_lscm_begin(phandle->charts[i], live);
1458 }
1459
1460 void param_lscm_solve(ParamHandle *handle)
1461 {
1462         PHandle *phandle = (PHandle*)handle;
1463         PChart *chart;
1464         int i;
1465         PBool result;
1466
1467         param_assert(phandle->state == PHANDLE_STATE_LSCM);
1468
1469         for (i = 0; i < phandle->ncharts; i++) {
1470                 chart = phandle->charts[i];
1471
1472                 if (chart->u.lscm.context) {
1473                         result = p_chart_lscm_solve(chart);
1474
1475                         if (!result || (chart->u.lscm.pin1))
1476                                 p_chart_lscm_end(chart);
1477                 }
1478         }
1479 }
1480
1481 void param_lscm_end(ParamHandle *handle)
1482 {
1483         PHandle *phandle = (PHandle*)handle;
1484         int i;
1485
1486         param_assert(phandle->state == PHANDLE_STATE_LSCM);
1487
1488         for (i = 0; i < phandle->ncharts; i++)
1489                 p_chart_lscm_end(phandle->charts[i]);
1490
1491         phandle->state = PHANDLE_STATE_CONSTRUCTED;
1492 }
1493
1494 /* Stretch */
1495
1496 #define P_STRETCH_ITER 20
1497
1498 static void p_stretch_pin_boundary(PChart *chart)
1499 {
1500         PVert *v;
1501
1502         for(v=(PVert*)chart->verts->first; v; v=v->link.next)
1503                 if (v->edge->pair == NULL)
1504                         v->flag |= PVERT_PIN;
1505                 else
1506                         v->flag &= ~PVERT_PIN;
1507 }
1508
1509 static float p_face_stretch(PFace *f)
1510 {
1511         float T, w, tmp[3];
1512         float Ps[3], Pt[3];
1513         float a, c, area;
1514         PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
1515         PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
1516
1517         area = p_face_uv_area_signed(f);
1518
1519         if (area <= 0.0f) /* flipped face -> infinite stretch */
1520                 return 1e10f;
1521
1522         if (f->flag & PFACE_FILLED)
1523                 return 0.0f;
1524
1525         w= 1.0f/(2.0f*area);
1526
1527         /* compute derivatives */
1528         VecCopyf(Ps, v1->co);
1529         VecMulf(Ps, (v2->uv[1] - v3->uv[1]));
1530
1531         VecCopyf(tmp, v2->co);
1532         VecMulf(tmp, (v3->uv[1] - v1->uv[1]));
1533         VecAddf(Ps, Ps, tmp);
1534
1535         VecCopyf(tmp, v3->co);
1536         VecMulf(tmp, (v1->uv[1] - v2->uv[1]));
1537         VecAddf(Ps, Ps, tmp);
1538
1539         VecMulf(Ps, w);
1540
1541         VecCopyf(Pt, v1->co);
1542         VecMulf(Pt, (v3->uv[0] - v2->uv[0]));
1543
1544         VecCopyf(tmp, v2->co);
1545         VecMulf(tmp, (v1->uv[0] - v3->uv[0]));
1546         VecAddf(Pt, Pt, tmp);
1547
1548         VecCopyf(tmp, v3->co);
1549         VecMulf(tmp, (v2->uv[0] - v1->uv[0]));
1550         VecAddf(Pt, Pt, tmp);
1551
1552         VecMulf(Pt, w);
1553
1554         /* Sander Tensor */
1555         a= Inpf(Ps, Ps);
1556         c= Inpf(Pt, Pt);
1557
1558         T = sqrt(0.5f*(a + c)*f->u.area3d);
1559
1560         return T;
1561 }
1562
1563 static float p_stretch_compute_vertex(PVert *v)
1564 {
1565         PEdge *e = v->edge;
1566         float sum = 0.0f;
1567
1568         do {
1569                 sum += p_face_stretch(e->face);
1570                 e = p_wheel_edge_next(e);
1571         } while (e && e != (v->edge));
1572
1573         return sum;
1574 }
1575
1576 static void p_chart_stretch_minimize(PChart *chart, RNG *rng)
1577 {
1578         PVert *v;
1579         PEdge *e;
1580         int j, nedges;
1581         float orig_stretch, low, stretch_low, high, stretch_high, mid, stretch;
1582         float orig_uv[2], dir[2], random_angle, trusted_radius;
1583
1584         for(v=(PVert*)chart->verts->first; v; v=v->link.next) {
1585                 if((v->flag & PVERT_PIN) || !(v->flag & PVERT_SELECT))
1586                         continue;
1587
1588                 orig_stretch = p_stretch_compute_vertex(v);
1589                 orig_uv[0] = v->uv[0];
1590                 orig_uv[1] = v->uv[1];
1591
1592                 /* move vertex in a random direction */
1593                 trusted_radius = 0.0f;
1594                 nedges = 0;
1595                 e = v->edge;
1596
1597                 do {
1598                         trusted_radius += p_edge_uv_length(e);
1599                         nedges++;
1600
1601                         e = p_wheel_edge_next(e);
1602                 } while (e && e != (v->edge));
1603
1604                 trusted_radius /= 2 * nedges;
1605
1606                 random_angle = rng_getFloat(rng) * 2.0 * M_PI;
1607                 dir[0] = trusted_radius * cos(random_angle);
1608                 dir[1] = trusted_radius * sin(random_angle);
1609
1610                 /* calculate old and new stretch */
1611                 low = 0;
1612                 stretch_low = orig_stretch;
1613
1614                 Vec2Addf(v->uv, orig_uv, dir);
1615                 high = 1;
1616                 stretch = stretch_high = p_stretch_compute_vertex(v);
1617
1618                 /* binary search for lowest stretch position */
1619                 for (j = 0; j < P_STRETCH_ITER; j++) {
1620                         mid = 0.5 * (low + high);
1621                         v->uv[0]= orig_uv[0] + mid*dir[0];
1622                         v->uv[1]= orig_uv[1] + mid*dir[1];
1623                         stretch = p_stretch_compute_vertex(v);
1624
1625                         if (stretch_low < stretch_high) {
1626                                 high = mid;
1627                                 stretch_high = stretch;
1628                         }
1629                         else {
1630                                 low = mid;
1631                                 stretch_low = stretch;
1632                         }
1633                 }
1634
1635                 /* no luck, stretch has increased, reset to old values */
1636                 if(stretch >= orig_stretch)
1637                         Vec2Copyf(v->uv, orig_uv);
1638         }
1639 }
1640
1641 void param_stretch_begin(ParamHandle *handle)
1642 {
1643         PHandle *phandle = (PHandle*)handle;
1644         PChart *chart;
1645         PVert *v;
1646         PFace *f;
1647         int i;
1648
1649         param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
1650         phandle->state = PHANDLE_STATE_STRETCH;
1651
1652         phandle->rng = rng_new(31415926);
1653         phandle->blend = 0.0f;
1654
1655         for (i = 0; i < phandle->ncharts; i++) {
1656                 chart = phandle->charts[i];
1657
1658                 for (v=(PVert*)chart->verts->first; v; v=v->link.next)
1659                         p_vert_load_select_uvs(v);
1660
1661                 p_stretch_pin_boundary(chart);
1662
1663                 for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
1664                         p_face_backup_uvs(f);
1665                         f->u.area3d = p_face_area(f);
1666                 }
1667         }
1668 }
1669
1670 void param_stretch_blend(ParamHandle *handle, float blend)
1671 {
1672         PHandle *phandle = (PHandle*)handle;
1673
1674         param_assert(phandle->state == PHANDLE_STATE_STRETCH);
1675         phandle->blend = blend;
1676 }
1677
1678 void param_stretch_iter(ParamHandle *handle)
1679 {
1680         PHandle *phandle = (PHandle*)handle;
1681         PChart *chart;
1682         int i;
1683
1684         param_assert(phandle->state == PHANDLE_STATE_STRETCH);
1685
1686         for (i = 0; i < phandle->ncharts; i++) {
1687                 chart = phandle->charts[i];
1688                 p_chart_stretch_minimize(chart, phandle->rng);
1689         }
1690 }
1691
1692 void param_stretch_end(ParamHandle *handle)
1693 {
1694         PHandle *phandle = (PHandle*)handle;
1695
1696         param_assert(phandle->state == PHANDLE_STATE_STRETCH);
1697         phandle->state = PHANDLE_STATE_CONSTRUCTED;
1698
1699         rng_free(phandle->rng);
1700         phandle->rng = NULL;
1701 }
1702
1703 /* Flushing */
1704
1705 void param_flush(ParamHandle *handle)
1706 {
1707         PHandle *phandle = (PHandle*)handle;
1708         PChart *chart;
1709         int i;
1710
1711         for (i = 0; i < phandle->ncharts; i++) {
1712                 chart = phandle->charts[i];
1713
1714                 if ((phandle->state == PHANDLE_STATE_LSCM) && !chart->u.lscm.context)
1715                         continue;
1716
1717                 if (phandle->blend == 0.0f)
1718                         p_flush_uvs(chart);
1719                 else
1720                         p_flush_uvs_blend(chart, phandle->blend);
1721         }
1722 }
1723
1724 void param_flush_restore(ParamHandle *handle)
1725 {
1726         PHandle *phandle = (PHandle*)handle;
1727         PChart *chart;
1728         PFace *f;
1729         int i;
1730
1731         for (i = 0; i < phandle->ncharts; i++) {
1732                 chart = phandle->charts[i];
1733
1734                 for (f=(PFace*)chart->faces->first; f; f=f->link.next)
1735                         p_face_restore_uvs(f);
1736         }
1737 }
1738
1739 /* Packing */
1740
1741 static int compare_chart_area(const void *a, const void *b)
1742 {
1743         PChart *ca = *((PChart**)a);
1744         PChart *cb = *((PChart**)b);
1745
1746     if (ca->u.pack.area > cb->u.pack.area)
1747                 return -1;
1748         else if (ca->u.pack.area == cb->u.pack.area)
1749                 return 0;
1750         else
1751                 return 1;
1752 }
1753
1754 static PBool p_pack_try(PHandle *handle, float side)
1755 {
1756         PChart *chart;
1757         float packx, packy, rowh, groupw, w, h;
1758         int i;
1759
1760         packx= packy= 0.0;
1761         rowh= 0.0;
1762         groupw= 1.0/sqrt(handle->ncharts);
1763
1764         for (i = 0; i < handle->ncharts; i++) {
1765                 chart = handle->charts[i];
1766
1767                 if (chart->flag & PCHART_NOPACK)
1768                         continue;
1769
1770                 w = chart->u.pack.size[0];
1771                 h = chart->u.pack.size[1];
1772
1773                 if(w <= (side-packx)) {
1774                         chart->u.pack.trans[0] = packx;
1775                         chart->u.pack.trans[1] = packy;
1776
1777                         packx += w;
1778                         rowh= MAX2(rowh, h);
1779                 }
1780                 else {
1781                         packy += rowh;
1782                         packx = w;
1783                         rowh = h;
1784
1785                         chart->u.pack.trans[0] = 0.0;
1786                         chart->u.pack.trans[1] = packy;
1787                 }
1788
1789                 if (packy+rowh > side)
1790                         return P_FALSE;
1791         }
1792
1793         return P_TRUE;
1794 }
1795
1796 #define PACK_SEARCH_DEPTH 15
1797
1798 void param_pack(ParamHandle *handle)
1799 {
1800         PHandle *phandle = (PHandle*)handle;
1801         PChart *chart;
1802         float uv_area, area, trans[2], minside, maxside, totarea, side;
1803         int i;
1804
1805         /* very simple rectangle packing */
1806
1807         if (phandle->ncharts == 0)
1808                 return;
1809
1810         totarea = 0.0f;
1811         maxside = 0.0f;
1812
1813         for (i = 0; i < phandle->ncharts; i++) {
1814                 chart = phandle->charts[i];
1815
1816                 if (chart->flag & PCHART_NOPACK) {
1817                         chart->u.pack.area = 0.0f;
1818                         continue;
1819                 }
1820
1821                 p_chart_area(chart, &uv_area, &area);
1822                 p_chart_uv_bbox(chart, trans, chart->u.pack.size);
1823
1824                 /* translate to origin and make area equal to 3d area */
1825                 chart->u.pack.rescale = (uv_area > 0.0f)? sqrt(area)/sqrt(uv_area): 0.0f;
1826                 chart->u.pack.area = area;
1827                 totarea += area;
1828
1829                 trans[0] = -trans[0];
1830                 trans[1] = -trans[1];
1831                 p_chart_uv_translate(chart, trans);
1832                 p_chart_uv_scale(chart, chart->u.pack.rescale);
1833
1834                 /* compute new dimensions for packing */
1835                 chart->u.pack.size[0] += trans[0];
1836                 chart->u.pack.size[1] += trans[1];
1837                 chart->u.pack.size[0] *= chart->u.pack.rescale;
1838                 chart->u.pack.size[1] *= chart->u.pack.rescale;
1839
1840                 maxside = MAX3(maxside, chart->u.pack.size[0], chart->u.pack.size[1]);
1841         }
1842
1843         /* sort by chart area, largest first */
1844         qsort(phandle->charts, phandle->ncharts, sizeof(PChart*), compare_chart_area);
1845
1846         /* binary search over pack region size */
1847         minside = MAX2(sqrt(totarea), maxside);
1848         maxside = (((int)sqrt(phandle->ncharts-1))+1)*maxside;
1849
1850         if (minside < maxside) { /* should always be true */
1851
1852                 for (i = 0; i < PACK_SEARCH_DEPTH; i++) {
1853                         if (p_pack_try(phandle, (minside+maxside)*0.5f + 1e-5))
1854                                 maxside = (minside+maxside)*0.5f;
1855                         else
1856                                 minside = (minside+maxside)*0.5f;
1857                 }
1858         }
1859
1860         /* do the actual packing */
1861         side = maxside + 1e-5;
1862         if (!p_pack_try(phandle, side))
1863                 param_warning("packing failed.\n");
1864
1865         for (i = 0; i < phandle->ncharts; i++) {
1866                 chart = phandle->charts[i];
1867
1868                 if (chart->flag & PCHART_NOPACK)
1869                         continue;
1870
1871                 p_chart_uv_scale(chart, 1.0f/side);
1872                 trans[0] = chart->u.pack.trans[0]/side;
1873                 trans[1] = chart->u.pack.trans[1]/side;
1874                 p_chart_uv_translate(chart, trans);
1875         }
1876 }
1877