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