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