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