Mesh Deform Modifier
[blender.git] / source / blender / src / parametrizer.c
index ff5d92e4604c939215bb7e012ba020ca2240121b..d0a51027ad311d88f2aaf435c1374cec562db387 100644 (file)
@@ -4,6 +4,8 @@
 #include "BLI_memarena.h"
 #include "BLI_arithb.h"
 #include "BLI_rand.h"
+#include "BLI_heap.h"
+#include "BLI_boxpack2d.h"
 
 #include "BKE_utildefines.h"
 
 #define M_PI 3.14159265358979323846
 #endif
 
-/* Hash */
+/* PHash
+   - special purpose hash that keeps all its elements in a single linked list.
+   - after construction, this hash is thrown away, and the list remains.
+   - removing elements is not possible efficiently.
+*/
 
 static int PHashSizes[] = {
        1, 3, 5, 11, 17, 37, 67, 131, 257, 521, 1031, 2053, 4099, 8209, 
@@ -33,13 +39,14 @@ static int PHashSizes[] = {
 };
 
 #define PHASH_hash(ph, item) (((unsigned long) (item))%((unsigned int) (ph)->cursize))
+#define PHASH_edge(v1, v2)      ((v1)^(v2))
 
-PHash *phash_new(int sizehint)
+static PHash *phash_new(PHashLink **list, int sizehint)
 {
        PHash *ph = (PHash*)MEM_callocN(sizeof(PHash), "PHash");
        ph->size = 0;
        ph->cursize_id = 0;
-       ph->first = NULL;
+       ph->list = list;
 
        while (PHashSizes[ph->cursize_id] < sizehint)
                ph->cursize_id++;
@@ -50,30 +57,18 @@ PHash *phash_new(int sizehint)
        return ph;
 }
 
-void phash_delete(PHash *ph)
+static void phash_delete(PHash *ph)
 {
        MEM_freeN(ph->buckets);
        MEM_freeN(ph);
 }
 
-void phash_delete_with_links(PHash *ph)
-{
-       PHashLink *link, *next=NULL;
-
-       for (link = ph->first; link; link = next) {
-               next = link->next;
-               MEM_freeN(link);
-       }
-
-       phash_delete(ph);
-}
-
-int phash_size(PHash *ph)
+static int phash_size(PHash *ph)
 {
        return ph->size;
 }
 
-void phash_insert(PHash *ph, PHashLink *link)
+static void phash_insert(PHash *ph, PHashLink *link)
 {
        int size = ph->cursize;
        int hash = PHASH_hash(ph, link->key);
@@ -82,8 +77,8 @@ void phash_insert(PHash *ph, PHashLink *link)
        if (lookup == NULL) {
                /* insert in front of the list */
                ph->buckets[hash] = link;
-               link->next = ph->first;
-               ph->first = link;
+               link->next = *(ph->list);
+               *(ph->list) = link;
        }
        else {
                /* insert after existing element */
@@ -94,13 +89,13 @@ void phash_insert(PHash *ph, PHashLink *link)
        ph->size++;
 
        if (ph->size > (size*3)) {
-               PHashLink *next = NULL, *first = ph->first;
+               PHashLink *next = NULL, *first = *(ph->list);
 
                ph->cursize = PHashSizes[++ph->cursize_id];
                MEM_freeN(ph->buckets);
                ph->buckets = (PHashLink**)MEM_callocN(ph->cursize*sizeof(*ph->buckets), "PHashBuckets");
                ph->size = 0;
-               ph->first = NULL;
+               *(ph->list) = NULL;
 
                for (link = first; link; link = next) {
                        next = link->next;
@@ -109,7 +104,7 @@ void phash_insert(PHash *ph, PHashLink *link)
        }
 }
 
-PHashLink *phash_lookup(PHash *ph, PHashKey key)
+static PHashLink *phash_lookup(PHash *ph, PHashKey key)
 {
        PHashLink *link;
        int hash = PHASH_hash(ph, key);
@@ -123,7 +118,7 @@ PHashLink *phash_lookup(PHash *ph, PHashKey key)
        return link;
 }
 
-PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
+static PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
 {
        int hash = PHASH_hash(ph, key);
 
@@ -136,208 +131,212 @@ PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
        return link;
 }
 
-/* Heap */
+/* Geometry */
 
-#define PHEAP_PARENT(i) ((i-1)>>1)
-#define PHEAP_LEFT(i)   ((i<<1)+1)
-#define PHEAP_RIGHT(i)  ((i<<1)+2)
-#define PHEAP_COMPARE(a, b) (a->value < b->value)
-#define PHEAP_EQUALS(a, b) (a->value == b->value)
-#define PHEAP_SWAP(heap, i, j) \
-       { SWAP(int, heap->tree[i]->index, heap->tree[j]->index); \
-         SWAP(PHeapLink*, heap->tree[i], heap->tree[j]);  }
-
-static void pheap_down(PHeap *heap, int i)
+static float p_vec_angle_cos(float *v1, float *v2, float *v3)
 {
-       while (P_TRUE) {
-               int size = heap->size, smallest;
-               int l = PHEAP_LEFT(i);
-               int r = PHEAP_RIGHT(i);
+       float d1[3], d2[3];
 
-               smallest = ((l < size) && PHEAP_COMPARE(heap->tree[l], heap->tree[i]))? l: i;
+       d1[0] = v1[0] - v2[0];
+       d1[1] = v1[1] - v2[1];
+       d1[2] = v1[2] - v2[2];
 
-               if ((r < size) && PHEAP_COMPARE(heap->tree[r], heap->tree[smallest]))
-                       smallest = r;
-               
-               if (smallest == i)
-                       break;
+       d2[0] = v3[0] - v2[0];
+       d2[1] = v3[1] - v2[1];
+       d2[2] = v3[2] - v2[2];
 
-               PHEAP_SWAP(heap, i, smallest);
-               i = smallest;
-       }
+       Normalize(d1);
+       Normalize(d2);
+
+       return d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2];
 }
 
-static void pheap_up(PHeap *heap, int i)
+static float p_vec_angle(float *v1, float *v2, float *v3)
 {
-       while (i > 0) {
-               int p = PHEAP_PARENT(i);
-
-               if (PHEAP_COMPARE(heap->tree[p], heap->tree[i]))
-                       break;
+       float dot = p_vec_angle_cos(v1, v2, v3);
 
-               PHEAP_SWAP(heap, p, i);
-               i = p;
-       }
+       if (dot <= -1.0f)
+               return (float)M_PI;
+       else if (dot >= 1.0f)
+               return 0.0f;
+       else
+               return (float)acos(dot);
 }
 
-PHeap *pheap_new()
+static float p_vec2_angle(float *v1, float *v2, float *v3)
 {
-       /* TODO: replace mallocN with something faster */
+       float u1[3], u2[3], u3[3];
 
-       PHeap *heap = (PHeap*)MEM_callocN(sizeof(PHeap), "PHeap");
-       heap->bufsize = 1;
-       heap->tree = (PHeapLink**)MEM_mallocN(sizeof(PHeapLink*), "PHeapTree");
+       u1[0] = v1[0]; u1[1] = v1[1]; u1[2] = 0.0f;
+       u2[0] = v2[0]; u2[1] = v2[1]; u2[2] = 0.0f;
+       u3[0] = v3[0]; u3[1] = v3[1]; u3[2] = 0.0f;
 
-       return heap;
+       return p_vec_angle(u1, u2, u3);
 }
 
-void pheap_delete(PHeap *heap)
+static void p_triangle_angles(float *v1, float *v2, float *v3, float *a1, float *a2, float *a3)
 {
-       MEM_freeN(heap->tree);
-       MEM_freeN(heap);
+       *a1 = p_vec_angle(v3, v1, v2);
+       *a2 = p_vec_angle(v1, v2, v3);
+       *a3 = M_PI - *a2 - *a1;
 }
 
-PHeapLink *pheap_insert(PHeap *heap, float value, void *ptr)
+static void p_face_angles(PFace *f, float *a1, float *a2, float *a3)
 {
-       PHeapLink *link;
+       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+       PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+
+       p_triangle_angles(v1->co, v2->co, v3->co, a1, a2, a3);
+}
 
-       if ((heap->size + 1) > heap->bufsize) {
-               int newsize = heap->bufsize*2;
+static float p_face_area(PFace *f)
+{
+       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+       PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
 
-               PHeapLink **ntree = (PHeapLink**)MEM_mallocN(newsize*sizeof(PHeapLink*), "PHeapTree");
-               memcpy(ntree, heap->tree, sizeof(PHeapLink*)*heap->size);
-               MEM_freeN(heap->tree);
+       return AreaT3Dfl(v1->co, v2->co, v3->co);
+}
 
-               heap->tree = ntree;
-               heap->bufsize = newsize;
-       }
+static float p_area_signed(float *v1, float *v2, float *v3)
+{
+       return 0.5f*(((v2[0] - v1[0])*(v3[1] - v1[1])) - 
+                   ((v3[0] - v1[0])*(v2[1] - v1[1])));
+}
 
-       param_assert(heap->size < heap->bufsize);
+static float p_face_uv_area_signed(PFace *f)
+{
+       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+       PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
 
-       link = MEM_mallocN(sizeof *link, "PHeapLink");
-       link->value = value;
-       link->ptr = ptr;
-       link->index = heap->size;
+       return 0.5f*(((v2->uv[0] - v1->uv[0])*(v3->uv[1] - v1->uv[1])) - 
+                   ((v3->uv[0] - v1->uv[0])*(v2->uv[1] - v1->uv[1])));
+}
 
-       heap->tree[link->index] = link;
+static float p_edge_length(PEdge *e)
+{
+    PVert *v1 = e->vert, *v2 = e->next->vert;
+    float d[3];
 
-       heap->size++;
+    d[0] = v2->co[0] - v1->co[0];
+    d[1] = v2->co[1] - v1->co[1];
+    d[2] = v2->co[2] - v1->co[2];
 
-       pheap_up(heap, heap->size-1);
+    return sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
+}
 
-       return link;
+static float p_edge_uv_length(PEdge *e)
+{
+    PVert *v1 = e->vert, *v2 = e->next->vert;
+    float d[3];
+
+    d[0] = v2->uv[0] - v1->uv[0];
+    d[1] = v2->uv[1] - v1->uv[1];
+
+    return sqrt(d[0]*d[0] + d[1]*d[1]);
 }
 
-int pheap_empty(PHeap *heap)
+static void p_chart_uv_bbox(PChart *chart, float *minv, float *maxv)
 {
-       return (heap->size == 0);
+       PVert *v;
+
+       INIT_MINMAX2(minv, maxv);
+
+       for (v=chart->verts; v; v=v->nextlink) {
+               DO_MINMAX2(v->uv, minv, maxv);
+       }
 }
 
-int pheap_size(PHeap *heap)
+static void p_chart_uv_scale(PChart *chart, float scale)
 {
-       return heap->size;
+       PVert *v;
+
+       for (v=chart->verts; v; v=v->nextlink) {
+               v->uv[0] *= scale;
+               v->uv[1] *= scale;
+       }
 }
 
-void *pheap_min(PHeap *heap)
+static void p_chart_uv_translate(PChart *chart, float trans[2])
 {
-       return heap->tree[0]->ptr;
+       PVert *v;
+
+       for (v=chart->verts; v; v=v->nextlink) {
+               v->uv[0] += trans[0];
+               v->uv[1] += trans[1];
+       }
 }
 
-void *pheap_popmin(PHeap *heap)
+static PBool p_intersect_line_2d_dir(float *v1, float *dir1, float *v2, float *dir2, float *isect)
 {
-       void *ptr = heap->tree[0]->ptr;
+       float lmbda, div;
 
-       MEM_freeN(heap->tree[0]);
+       div= dir2[0]*dir1[1] - dir2[1]*dir1[0];
 
-       if (heap->size == 1)
-               heap->size--;
-       else {
-               PHEAP_SWAP(heap, 0, heap->size-1);
-               heap->size--;
+       if (div == 0.0f)
+               return P_FALSE;
 
-               pheap_down(heap, 0);
-       }
+       lmbda= ((v1[1]-v2[1])*dir1[0]-(v1[0]-v2[0])*dir1[1])/div;
+       isect[0] = v1[0] + lmbda*dir2[0];
+       isect[1] = v1[1] + lmbda*dir2[1];
 
-       return ptr;
+       return P_TRUE;
 }
 
-static void pheap_remove(PHeap *heap, PHeapLink *link)
+#if 0
+static PBool p_intersect_line_2d(float *v1, float *v2, float *v3, float *v4, float *isect)
 {
-       int i = link->index;
+       float dir1[2], dir2[2];
 
-       while (i > 0) {
-               int p = PHEAP_PARENT(i);
+       dir1[0] = v4[0] - v3[0];
+       dir1[1] = v4[1] - v3[1];
 
-               PHEAP_SWAP(heap, p, i);
-               i = p;
+       dir2[0] = v2[0] - v1[0];
+       dir2[1] = v2[1] - v1[1];
+
+       if (!p_intersect_line_2d_dir(v1, dir1, v2, dir2, isect)) {
+               /* parallel - should never happen in theory for polygon kernel, but
+                  let's give a point nearby in case things go wrong */
+               isect[0] = (v1[0] + v2[0])*0.5f;
+               isect[1] = (v1[1] + v2[1])*0.5f;
+               return P_FALSE;
        }
 
-       pheap_popmin(heap);
+       return P_TRUE;
 }
+#endif
 
-/* Construction */
+/* Topological Utilities */
 
-PEdge *p_wheel_edge_next(PEdge *e)
+static PEdge *p_wheel_edge_next(PEdge *e)
 {
        return e->next->next->pair;
 }
 
-PEdge *p_wheel_edge_prev(PEdge *e)
+static PEdge *p_wheel_edge_prev(PEdge *e)
 {   
        return (e->pair)? e->pair->next: NULL;
 }
 
-static PVert *p_vert_add(PChart *chart, PHashKey key, float *co, PEdge *e)
-{
-       PVert *v = (PVert*)BLI_memarena_alloc(chart->handle->arena, sizeof *v);
-       v->co = co;
-       v->link.key = key;
-       v->edge = e;
-
-       phash_insert(chart->verts, (PHashLink*)v);
-
-       return v;
-}
-
-static PVert *p_vert_lookup(PChart *chart, PHashKey key, float *co, PEdge *e)
+static PEdge *p_boundary_edge_next(PEdge *e)
 {
-       PVert *v = (PVert*)phash_lookup(chart->verts, key);
-
-       if (v)
-               return v;
-       else
-               return p_vert_add(chart, key, co, e);
+       return e->next->vert->edge;
 }
 
-static PVert *p_vert_copy(PChart *chart, PVert *v)
+static PEdge *p_boundary_edge_prev(PEdge *e)
 {
-       PVert *nv = (PVert*)BLI_memarena_alloc(chart->handle->arena, sizeof *nv);
-       nv->co = v->co;
-       nv->uv[0] = v->uv[0];
-       nv->uv[1] = v->uv[1];
-       nv->link.key = v->link.key;
-       nv->edge = v->edge;
+    PEdge *we = e, *last;
 
-       phash_insert(chart->verts, (PHashLink*)nv);
+       do {
+               last = we;
+               we = p_wheel_edge_next(we);
+       } while (we && (we != e));
 
-       return nv;
+       return last->next->next;
 }
 
-static PEdge *p_edge_lookup(PChart *chart, PHashKey *vkeys)
+static PBool p_vert_interior(PVert *v)
 {
-       PHashKey key = vkeys[0]^vkeys[1];
-       PEdge *e = (PEdge*)phash_lookup(chart->edges, key);
-
-       while (e) {
-               if ((e->vert->link.key == vkeys[0]) && (e->next->vert->link.key == vkeys[1]))
-                       return e;
-               else if ((e->vert->link.key == vkeys[1]) && (e->next->vert->link.key == vkeys[0]))
-                       return e;
-
-               e = (PEdge*)phash_next(chart->edges, key, (PHashLink*)e);
-       }
-
-       return NULL;
+       return (v->edge->pair != NULL);
 }
 
 static void p_face_flip(PFace *f)
@@ -359,202 +358,201 @@ static void p_face_flip(PFace *f)
        e3->flag = (f3 & ~PEDGE_VERTEX_FLAGS) | (f1 & PEDGE_VERTEX_FLAGS);
 }
 
-static void p_vert_load_pin_uvs(PVert *v)
+#if 0
+static void p_chart_topological_sanity_check(PChart *chart)
 {
+       PVert *v;
        PEdge *e;
-       int nedges = 0;
 
-       v->uv[0] = v->uv[1] = 0.0f;
-       nedges = 0;
-       e = v->edge;
-       do {
-               if (e->orig_uv && (e->flag & PEDGE_PIN)) {
-                       v->flag |= PVERT_PIN;
-                       v->uv[0] += e->orig_uv[0];
-                       v->uv[1] += e->orig_uv[1];
-                       nedges++;
+       for (v=chart->verts; v; v=v->nextlink)
+               param_test_equals_ptr("v->edge->vert", v, v->edge->vert);
+       
+       for (e=chart->edges; e; e=e->nextlink) {
+               if (e->pair) {
+                       param_test_equals_ptr("e->pair->pair", e, e->pair->pair);
+                       param_test_equals_ptr("pair->vert", e->vert, e->pair->next->vert);
+                       param_test_equals_ptr("pair->next->vert", e->next->vert, e->pair->vert);
                }
-
-               e = p_wheel_edge_next(e);
-       } while (e && e != (v->edge));
-
-       if (nedges > 0) {
-               v->uv[0] /= nedges;
-               v->uv[1] /= nedges;
        }
 }
+#endif
 
-static void p_vert_load_select_uvs(PVert *v)
+/* Loading / Flushing */
+
+static void p_vert_load_pin_select_uvs(PVert *v)
 {
        PEdge *e;
-       int nedges = 0;
+       int nedges = 0, npins = 0;
+       float pinuv[2];
 
        v->uv[0] = v->uv[1] = 0.0f;
-       nedges = 0;
+       pinuv[0] = pinuv[1] = 0.0f;
        e = v->edge;
        do {
-               if (e->orig_uv && (e->flag & PEDGE_SELECT))
-                       v->flag |= PVERT_SELECT;
+               if (e->orig_uv) {
+                       if (e->flag & PEDGE_SELECT)
+                               v->flag |= PVERT_SELECT;
+
+                       if (e->flag & PEDGE_PIN) {
+                               pinuv[0] += e->orig_uv[0];
+                               pinuv[1] += e->orig_uv[1];
+                               npins++;
+                       }
+                       else {
+                               v->uv[0] += e->orig_uv[0];
+                               v->uv[1] += e->orig_uv[1];
+                       }
 
-               v->uv[0] += e->orig_uv[0];
-               v->uv[1] += e->orig_uv[1];
-               nedges++;
+                       nedges++;
+               }
 
                e = p_wheel_edge_next(e);
        } while (e && e != (v->edge));
 
-       if (nedges > 0) {
+       if (npins > 0) {
+               v->uv[0] = pinuv[0]/npins;
+               v->uv[1] = pinuv[1]/npins;
+               v->flag |= PVERT_PIN;
+       }
+       else if (nedges > 0) {
                v->uv[0] /= nedges;
                v->uv[1] /= nedges;
        }
 }
 
-static void p_extrema_verts(PChart *chart, PVert **v1, PVert **v2)
+static void p_flush_uvs(PChart *chart)
 {
-       float minv[3], maxv[3], dirlen;
-       PVert *v, *minvert[3], *maxvert[3];
-       int i, dir;
-
-       /* find minimum and maximum verts over x/y/z axes */
-       minv[0] = minv[1] = minv[2] = 1e20;
-       maxv[0] = maxv[1] = maxv[2] = -1e20;
-
-       minvert[0] = minvert[1] = minvert[2] = NULL;
-       maxvert[0] = maxvert[1] = maxvert[2] = NULL;
+       PEdge *e;
 
-       for (v = (PVert*)chart->verts->first; v; v=v->link.next) {
-               for (i = 0; i < 3; i++) {
-                       if (v->co[i] < minv[i]) {
-                               minv[i] = v->co[i];
-                               minvert[i] = v;
-                       }
-                       if (v->co[i] > maxv[i]) {
-                               maxv[i] = v->co[i];
-                               maxvert[i] = v;
-                       }
+       for (e=chart->edges; e; e=e->nextlink) {
+               if (e->orig_uv) {
+                       e->orig_uv[0] = e->vert->uv[0];
+                       e->orig_uv[1] = e->vert->uv[1];
                }
        }
+}
 
-       /* find axes with longest distance */
-       dir = 0;
-       dirlen = -1.0;
+static void p_flush_uvs_blend(PChart *chart, float blend)
+{
+       PEdge *e;
+       float invblend = 1.0f - blend;
 
-       for (i = 0; i < 3; i++) {
-               if (maxv[i] - minv[i] > dirlen) {
-                       dir = i;
-                       dirlen = maxv[i] - minv[i];
+       for (e=chart->edges; e; e=e->nextlink) {
+               if (e->orig_uv) {
+                       e->orig_uv[0] = blend*e->old_uv[0] + invblend*e->vert->uv[0];
+                       e->orig_uv[1] = blend*e->old_uv[1] + invblend*e->vert->uv[1];
                }
        }
-
-       if (minvert[dir] == maxvert[dir]) {
-               /* degenerate case */
-               PFace *f = (PFace*)chart->faces->first;
-               *v1 = f->edge->vert;
-               *v2 = f->edge->next->vert;
-       }
-       else {
-               *v1 = minvert[dir];
-               *v2 = maxvert[dir];
-       }
 }
 
-static float p_vec_normalise(float *v)
+static void p_face_backup_uvs(PFace *f)
 {
-   float d;
-   
-    d = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
-
-    if(d != 0.0f) {
-               d = 1.0f/d;
+       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 
-               v[0] *= d;
-               v[1] *= d;
-               v[2] *= d;
+       if (e1->orig_uv && e2->orig_uv && e3->orig_uv) {
+               e1->old_uv[0] = e1->orig_uv[0];
+               e1->old_uv[1] = e1->orig_uv[1];
+               e2->old_uv[0] = e2->orig_uv[0];
+               e2->old_uv[1] = e2->orig_uv[1];
+               e3->old_uv[0] = e3->orig_uv[0];
+               e3->old_uv[1] = e3->orig_uv[1];
        }
-
-       return d;
 }
 
-static float p_vec_angle_cos(float *v1, float *v2, float *v3)
+static void p_face_restore_uvs(PFace *f)
 {
-       float d1[3], d2[3];
+       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 
-       d1[0] = v1[0] - v2[0];
-       d1[1] = v1[1] - v2[1];
-       d1[2] = v1[2] - v2[2];
+       if (e1->orig_uv && e2->orig_uv && e3->orig_uv) {
+               e1->orig_uv[0] = e1->old_uv[0];
+               e1->orig_uv[1] = e1->old_uv[1];
+               e2->orig_uv[0] = e2->old_uv[0];
+               e2->orig_uv[1] = e2->old_uv[1];
+               e3->orig_uv[0] = e3->old_uv[0];
+               e3->orig_uv[1] = e3->old_uv[1];
+       }
+}
 
-       d2[0] = v3[0] - v2[0];
-       d2[1] = v3[1] - v2[1];
-       d2[2] = v3[2] - v2[2];
+/* Construction (use only during construction, relies on u.key being set */
 
-       p_vec_normalise(d1);
-       p_vec_normalise(d2);
+static PVert *p_vert_add(PHandle *handle, PHashKey key, float *co, PEdge *e)
+{
+       PVert *v = (PVert*)BLI_memarena_alloc(handle->arena, sizeof *v);
+       v->co = co;
+       v->u.key = key;
+       v->edge = e;
+       v->flag = 0;
 
-       return d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2];
+       phash_insert(handle->hash_verts, (PHashLink*)v);
+
+       return v;
 }
 
-static float p_vec_angle(float *v1, float *v2, float *v3)
+static PVert *p_vert_lookup(PHandle *handle, PHashKey key, float *co, PEdge *e)
 {
-       float dot = p_vec_angle_cos(v1, v2, v3);
+       PVert *v = (PVert*)phash_lookup(handle->hash_verts, key);
 
-       if (dot <= -1.0f)
-               return (float)M_PI;
-       else if (dot >= 1.0f)
-               return 0.0f;
+       if (v)
+               return v;
        else
-               return (float)acos(dot);
+               return p_vert_add(handle, key, co, e);
 }
 
-static void p_face_angles(PFace *f, float *a1, float *a2, float *a3)
+static PVert *p_vert_copy(PChart *chart, PVert *v)
 {
-       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-       PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+       PVert *nv = (PVert*)BLI_memarena_alloc(chart->handle->arena, sizeof *nv);
 
-       *a1 = p_vec_angle(v3->co, v1->co, v2->co);
-       *a2 = p_vec_angle(v1->co, v2->co, v3->co);
-       *a3 = M_PI - *a2 - *a1;
+       nv->co = v->co;
+       nv->uv[0] = v->uv[0];
+       nv->uv[1] = v->uv[1];
+       nv->u.key = v->u.key;
+       nv->edge = v->edge;
+       nv->flag = v->flag;
+
+       return nv;
 }
 
-static float p_face_area(PFace *f)
+static PEdge *p_edge_lookup(PHandle *handle, PHashKey *vkeys)
 {
-       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-       PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+       PHashKey key = PHASH_edge(vkeys[0], vkeys[1]);
+       PEdge *e = (PEdge*)phash_lookup(handle->hash_edges, key);
 
-       return AreaT3Dfl(v1->co, v2->co, v3->co);
-}
+       while (e) {
+               if ((e->vert->u.key == vkeys[0]) && (e->next->vert->u.key == vkeys[1]))
+                       return e;
+               else if ((e->vert->u.key == vkeys[1]) && (e->next->vert->u.key == vkeys[0]))
+                       return e;
 
-static float p_face_uv_area_signed(PFace *f)
-{
-       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-       PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+               e = (PEdge*)phash_next(handle->hash_edges, key, (PHashLink*)e);
+       }
 
-       return 0.5f*(((v2->uv[0]-v1->uv[0]) * (v3->uv[1]-v1->uv[1])) - 
-                   ((v3->uv[0]-v1->uv[0]) * (v2->uv[1]-v1->uv[1])));
+       return NULL;
 }
 
-static float p_face_uv_area(PFace *f)
+static PBool p_face_exists(PHandle *handle, PHashKey *vkeys, int i1, int i2, int i3)
 {
-       return fabs(p_face_uv_area_signed(f));
-}
+       PHashKey key = PHASH_edge(vkeys[i1], vkeys[i2]);
+       PEdge *e = (PEdge*)phash_lookup(handle->hash_edges, key);
 
-static void p_chart_area(PChart *chart, float *uv_area, float *area)
-{
-       PFace *f;
-
-       *uv_area = *area = 0.0f;
+       while (e) {
+               if ((e->vert->u.key == vkeys[i1]) && (e->next->vert->u.key == vkeys[i2])) {
+                       if (e->next->next->vert->u.key == vkeys[i3])
+                               return P_TRUE;
+               }
+               else if ((e->vert->u.key == vkeys[i2]) && (e->next->vert->u.key == vkeys[i1])) {
+                       if (e->next->next->vert->u.key == vkeys[i3])
+                               return P_TRUE;
+               }
 
-       for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
-               *uv_area += p_face_uv_area(f);
-               *area += p_face_area(f);
+               e = (PEdge*)phash_next(handle->hash_edges, key, (PHashLink*)e);
        }
+
+       return P_FALSE;
 }
 
 static PChart *p_chart_new(PHandle *handle)
 {
        PChart *chart = (PChart*)MEM_callocN(sizeof*chart, "PChart");
-       chart->verts = phash_new(1);
-       chart->edges = phash_new(1);
-       chart->faces = phash_new(1);
        chart->handle = handle;
 
        return chart;
@@ -563,10 +561,6 @@ static PChart *p_chart_new(PHandle *handle)
 static void p_chart_delete(PChart *chart)
 {
        /* the actual links are free by memarena */
-       phash_delete(chart->verts);
-       phash_delete(chart->edges);
-       phash_delete(chart->faces);
-
        MEM_freeN(chart);
 }
 
@@ -575,10 +569,13 @@ static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
        float *uv1, *uv2, *uvp1, *uvp2;
        float limit[2];
 
+       limit[0] = 0.00001;
+       limit[1] = 0.00001;
+
        uv1 = e->orig_uv;
        uv2 = e->next->orig_uv;
 
-       if (e->vert->link.key == ep->vert->link.key) {
+       if (e->vert->u.key == ep->vert->u.key) {
                uvp1 = ep->orig_uv;
                uvp2 = ep->next->orig_uv;
        }
@@ -587,14 +584,12 @@ static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
                uvp2 = ep->orig_uv;
        }
 
-       get_connected_limit_tface_uv(limit);
-
-       if((fabs(uv1[0]-uvp1[0]) > limit[0]) && (fabs(uv1[1]-uvp1[1]) > limit[1])) {
+       if((fabs(uv1[0]-uvp1[0]) > limit[0]) || (fabs(uv1[1]-uvp1[1]) > limit[1])) {
                e->flag |= PEDGE_SEAM;
                ep->flag |= PEDGE_SEAM;
                return P_TRUE;
        }
-       if((fabs(uv2[0]-uvp2[0]) > limit[0]) && (fabs(uv2[1]-uvp2[1]) > limit[1])) {
+       if((fabs(uv2[0]-uvp2[0]) > limit[0]) || (fabs(uv2[1]-uvp2[1]) > limit[1])) {
                e->flag |= PEDGE_SEAM;
                ep->flag |= PEDGE_SEAM;
                return P_TRUE;
@@ -603,19 +598,19 @@ static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
        return P_FALSE;
 }
 
-static PBool p_edge_has_pair(PChart *chart, PEdge *e, PEdge **pair, PBool impl)
+static PBool p_edge_has_pair(PHandle *handle, PEdge *e, PEdge **pair, PBool impl)
 {
        PHashKey key;
        PEdge *pe;
        PVert *v1, *v2;
-       PHashKey key1 = e->vert->link.key;
-       PHashKey key2 = e->next->vert->link.key;
+       PHashKey key1 = e->vert->u.key;
+       PHashKey key2 = e->next->vert->u.key;
 
        if (e->flag & PEDGE_SEAM)
                return P_FALSE;
        
-       key = key1 ^ key2;
-       pe = (PEdge*)phash_lookup(chart->edges, key);
+       key = PHASH_edge(key1, key2);
+       pe = (PEdge*)phash_lookup(handle->hash_edges, key);
        *pair = NULL;
 
        while (pe) {
@@ -623,8 +618,8 @@ static PBool p_edge_has_pair(PChart *chart, PEdge *e, PEdge **pair, PBool impl)
                        v1 = pe->vert;
                        v2 = pe->next->vert;
 
-                       if (((v1->link.key == key1) && (v2->link.key == key2)) ||
-                               ((v1->link.key == key2) && (v2->link.key == key1))) {
+                       if (((v1->u.key == key1) && (v2->u.key == key2)) ||
+                               ((v1->u.key == key2) && (v2->u.key == key1))) {
 
                                /* don't connect seams and t-junctions */
                                if ((pe->flag & PEDGE_SEAM) || *pair ||
@@ -637,7 +632,7 @@ static PBool p_edge_has_pair(PChart *chart, PEdge *e, PEdge **pair, PBool impl)
                        }
                }
 
-               pe = (PEdge*)phash_next(chart->edges, key, (PHashLink*)pe);
+               pe = (PEdge*)phash_next(handle->hash_edges, key, (PHashLink*)pe);
        }
 
        if (*pair && (e->vert == (*pair)->vert)) {
@@ -651,11 +646,11 @@ static PBool p_edge_has_pair(PChart *chart, PEdge *e, PEdge **pair, PBool impl)
        return (*pair != NULL);
 }
 
-static PBool p_edge_connect_pair(PChart *chart, PEdge *e, PEdge ***stack, PBool impl)
+static PBool p_edge_connect_pair(PHandle *handle, PEdge *e, PEdge ***stack, PBool impl)
 {
-       PEdge *pair;
+       PEdge *pair = NULL;
 
-       if(!e->pair && p_edge_has_pair(chart, e, &pair, impl)) {
+       if(!e->pair && p_edge_has_pair(handle, e, &pair, impl)) {
                if (e->vert == pair->vert)
                        p_face_flip(pair->face);
 
@@ -671,16 +666,17 @@ static PBool p_edge_connect_pair(PChart *chart, PEdge *e, PEdge ***stack, PBool
        return (e->pair != NULL);
 }
 
-static int p_connect_pairs(PChart *chart, PBool impl)
+static int p_connect_pairs(PHandle *handle, PBool impl)
 {
-       PEdge **stackbase = MEM_mallocN(sizeof*stackbase * phash_size(chart->faces), "Pstackbase");
+       PEdge **stackbase = MEM_mallocN(sizeof*stackbase*phash_size(handle->hash_faces), "Pstackbase");
        PEdge **stack = stackbase;
        PFace *f, *first;
        PEdge *e, *e1, *e2;
+       PChart *chart = handle->construction_chart;
        int ncharts = 0;
 
        /* connect pairs, count edges, set vertex-edge pointer to a pairless edge */
-       for (first=(PFace*)chart->faces->first; first; first=first->link.next) {
+       for (first=chart->faces; first; first=first->nextlink) {
                if (first->flag & PFACE_CONNECTED)
                        continue;
 
@@ -699,11 +695,11 @@ static int p_connect_pairs(PChart *chart, PBool impl)
                        /* assign verts to charts so we can sort them later */
                        f->u.chart = ncharts;
 
-                       if (!p_edge_connect_pair(chart, e, &stack, impl))
+                       if (!p_edge_connect_pair(handle, e, &stack, impl))
                                e->vert->edge = e;
-                       if (!p_edge_connect_pair(chart, e1, &stack, impl))
+                       if (!p_edge_connect_pair(handle, e1, &stack, impl))
                                e1->vert->edge = e1;
-                       if (!p_edge_connect_pair(chart, e2, &stack, impl))
+                       if (!p_edge_connect_pair(handle, e2, &stack, impl))
                                e2->vert->edge = e2;
                }
 
@@ -739,13 +735,22 @@ static void p_split_vert(PChart *chart, PEdge *e)
                if (we == v->edge) {
                        /* found it, no need to copy */
                        copy = P_FALSE;
-                       phash_insert(chart->verts, (PHashLink*)v);
+                       v->nextlink = chart->verts;
+                       chart->verts = v;
+                       chart->nverts++;
                }
        }
 
        if (copy) {
                /* not found, copying */
+               v->flag |= PVERT_SPLIT;
                v = p_vert_copy(chart, v);
+               v->flag |= PVERT_SPLIT;
+
+               v->nextlink = chart->verts;
+               chart->verts = v;
+               chart->nverts++;
+
                v->edge = lastwe;
 
                we = lastwe;
@@ -765,17 +770,24 @@ static PChart **p_split_charts(PHandle *handle, PChart *chart, int ncharts)
        for (i = 0; i < ncharts; i++)
                charts[i] = p_chart_new(handle);
 
-       f = (PFace*)chart->faces->first;
+       f = chart->faces;
        while (f) {
                PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-               nextf = f->link.next;
+               nextf = f->nextlink;
 
                nchart = charts[f->u.chart];
 
-               phash_insert(nchart->faces, (PHashLink*)f);
-               phash_insert(nchart->edges, (PHashLink*)e1);
-               phash_insert(nchart->edges, (PHashLink*)e2);
-               phash_insert(nchart->edges, (PHashLink*)e3);
+               f->nextlink = nchart->faces;
+               nchart->faces = f;
+               e1->nextlink = nchart->edges;
+               nchart->edges = e1;
+               e2->nextlink = nchart->edges;
+               nchart->edges = e2;
+               e3->nextlink = nchart->edges;
+               nchart->edges = e3;
+
+               nchart->nfaces++;
+               nchart->nedges += 3;
 
                p_split_vert(nchart, e1);
                p_split_vert(nchart, e2);
@@ -787,43 +799,18 @@ static PChart **p_split_charts(PHandle *handle, PChart *chart, int ncharts)
        return charts;
 }
 
-static void p_face_backup_uvs(PFace *f)
-{
-       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-
-       e1->old_uv[0] = e1->orig_uv[0];
-       e1->old_uv[1] = e1->orig_uv[1];
-       e2->old_uv[0] = e2->orig_uv[0];
-       e2->old_uv[1] = e2->orig_uv[1];
-       e3->old_uv[0] = e3->orig_uv[0];
-       e3->old_uv[1] = e3->orig_uv[1];
-}
-
-static void p_face_restore_uvs(PFace *f)
-{
-       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-
-       e1->orig_uv[0] = e1->old_uv[0];
-       e1->orig_uv[1] = e1->old_uv[1];
-       e2->orig_uv[0] = e2->old_uv[0];
-       e2->orig_uv[1] = e2->old_uv[1];
-       e3->orig_uv[0] = e3->old_uv[0];
-       e3->orig_uv[1] = e3->old_uv[1];
-}
-
-static PFace *p_face_add(PChart *chart, ParamKey key, ParamKey *vkeys,
-                         float *co[3], float *uv[3], int i1, int i2, int i3,
-                         ParamBool *pin, ParamBool *select)
+static PFace *p_face_add(PHandle *handle)
 {
        PFace *f;
        PEdge *e1, *e2, *e3;
 
        /* allocate */
-       f = (PFace*)BLI_memarena_alloc(chart->handle->arena, sizeof *f);
+       f = (PFace*)BLI_memarena_alloc(handle->arena, sizeof *f);
+       f->flag=0; // init !
 
-       e1 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e1);
-       e2 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e2);
-       e3 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e3);
+       e1 = (PEdge*)BLI_memarena_alloc(handle->arena, sizeof *e1);
+       e2 = (PEdge*)BLI_memarena_alloc(handle->arena, sizeof *e2);
+       e3 = (PEdge*)BLI_memarena_alloc(handle->arena, sizeof *e3);
 
        /* set up edges */
        f->edge = e1;
@@ -833,20 +820,31 @@ static PFace *p_face_add(PChart *chart, ParamKey key, ParamKey *vkeys,
        e2->next = e3;
        e3->next = e1;
 
-       if (co && uv) {
-               e1->vert = p_vert_lookup(chart, vkeys[i1], co[i1], e1);
-               e2->vert = p_vert_lookup(chart, vkeys[i2], co[i2], e2);
-               e3->vert = p_vert_lookup(chart, vkeys[i3], co[i3], e3);
+       e1->pair = NULL;
+       e2->pair = NULL;
+       e3->pair = NULL;
+   
+       e1->flag =0;
+       e2->flag =0;
+       e3->flag =0;
 
-               e1->orig_uv = uv[i1];
-               e2->orig_uv = uv[i2];
-               e3->orig_uv = uv[i3];
-       }
-       else {
-               /* internal call to add face */
-               e1->vert = e2->vert = e3->vert = NULL;
-               e1->orig_uv = e2->orig_uv = e3->orig_uv = NULL;
-       }
+       return f;
+}
+
+static PFace *p_face_add_construct(PHandle *handle, ParamKey key, ParamKey *vkeys,
+                                   float *co[3], float *uv[3], int i1, int i2, int i3,
+                                   ParamBool *pin, ParamBool *select)
+{
+       PFace *f = p_face_add(handle);
+       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+
+       e1->vert = p_vert_lookup(handle, vkeys[i1], co[i1], e1);
+       e2->vert = p_vert_lookup(handle, vkeys[i2], co[i2], e2);
+       e3->vert = p_vert_lookup(handle, vkeys[i3], co[i3], e3);
+
+       e1->orig_uv = uv[i1];
+       e2->orig_uv = uv[i2];
+       e3->orig_uv = uv[i3];
 
        if (pin) {
                if (pin[i1]) e1->flag |= PEDGE_PIN;
@@ -861,102 +859,86 @@ static PFace *p_face_add(PChart *chart, ParamKey key, ParamKey *vkeys,
        }
 
        /* insert into hash */
-       f->link.key = key;
-       phash_insert(chart->faces, (PHashLink*)f);
+       f->u.key = key;
+       phash_insert(handle->hash_faces, (PHashLink*)f);
 
-       e1->link.key = vkeys[i1]^vkeys[i2];
-       e2->link.key = vkeys[i2]^vkeys[i3];
-       e3->link.key = vkeys[i3]^vkeys[i1];
+       e1->u.key = PHASH_edge(vkeys[i1], vkeys[i2]);
+       e2->u.key = PHASH_edge(vkeys[i2], vkeys[i3]);
+       e3->u.key = PHASH_edge(vkeys[i3], vkeys[i1]);
 
-       phash_insert(chart->edges, (PHashLink*)e1);
-       phash_insert(chart->edges, (PHashLink*)e2);
-       phash_insert(chart->edges, (PHashLink*)e3);
+       phash_insert(handle->hash_edges, (PHashLink*)e1);
+       phash_insert(handle->hash_edges, (PHashLink*)e2);
+       phash_insert(handle->hash_edges, (PHashLink*)e3);
 
        return f;
 }
 
-static PBool p_quad_split_direction(float **co)
-{
-    float a1, a2;
-       
-       a1 = p_vec_angle_cos(co[0], co[1], co[2]);
-       a1 += p_vec_angle_cos(co[1], co[0], co[2]);
-       a1 += p_vec_angle_cos(co[2], co[0], co[1]);
-
-       a2 = p_vec_angle_cos(co[0], co[1], co[3]);
-       a2 += p_vec_angle_cos(co[1], co[0], co[3]);
-       a2 += p_vec_angle_cos(co[3], co[0], co[1]);
-
-       return (a1 > a2);
-}
-
-static float p_edge_length(PEdge *e)
+static PFace *p_face_add_fill(PChart *chart, PVert *v1, PVert *v2, PVert *v3)
 {
-    PVert *v1 = e->vert, *v2 = e->next->vert;
-    float d[3];
+       PFace *f = p_face_add(chart->handle);
+       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 
-    d[0] = v2->co[0] - v1->co[0];
-    d[1] = v2->co[1] - v1->co[1];
-    d[2] = v2->co[2] - v1->co[2];
+       e1->vert = v1;
+       e2->vert = v2;
+       e3->vert = v3;
 
-    return sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
-}
+       e1->orig_uv = e2->orig_uv = e3->orig_uv = NULL;
 
-static float p_edge_uv_length(PEdge *e)
-{
-    PVert *v1 = e->vert, *v2 = e->next->vert;
-    float d[3];
+       f->nextlink = chart->faces;
+       chart->faces = f;
+       e1->nextlink = chart->edges;
+       chart->edges = e1;
+       e2->nextlink = chart->edges;
+       chart->edges = e2;
+       e3->nextlink = chart->edges;
+       chart->edges = e3;
 
-    d[0] = v2->uv[0] - v1->uv[0];
-    d[1] = v2->uv[1] - v1->uv[1];
+       chart->nfaces++;
+       chart->nedges += 3;
 
-    return sqrt(d[0]*d[0] + d[1]*d[1]);
+       return f;
 }
 
-void p_chart_uv_bbox(PChart *chart, float *minv, float *maxv)
+static PBool p_quad_split_direction(PHandle *handle, float **co, PHashKey *vkeys)
 {
-       PVert *v;
-
-       INIT_MINMAX2(minv, maxv);
-
-       for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
-               DO_MINMAX2(v->uv, minv, maxv);
+       float fac= VecLenf(co[0], co[2]) - VecLenf(co[1], co[3]);
+       PBool dir = (fac <= 0.0f);
+
+       /* the face exists check is there because of a special case: when
+          two quads share three vertices, they can each be split into two
+          triangles, resulting in two identical triangles. for example in
+          suzanne's nose. */
+       if (dir) {
+               if (p_face_exists(handle,vkeys,0,1,2) || p_face_exists(handle,vkeys,0,2,3))
+                       return !dir;
        }
-}
-
-static void p_chart_uv_scale(PChart *chart, float scale)
-{
-       PVert *v;
-
-       for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
-               v->uv[0] *= scale;
-               v->uv[1] *= scale;
+       else {
+               if (p_face_exists(handle,vkeys,0,1,3) || p_face_exists(handle,vkeys,1,2,3))
+                       return !dir;
        }
-}
-
-static void p_chart_uv_translate(PChart *chart, float trans[2])
-{
-       PVert *v;
 
-       for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
-               v->uv[0] += trans[0];
-               v->uv[1] += trans[1];
-       }
+       return dir;
 }
 
+/* Construction: boundary filling */
+
 static void p_chart_boundaries(PChart *chart, int *nboundaries, PEdge **outer)
 {   
     PEdge *e, *be;
     float len, maxlen = -1.0;
 
-       *nboundaries = 0;
-       *outer = NULL;
+       if (nboundaries)
+               *nboundaries = 0;
+       if (outer)
+               *outer = NULL;
 
-       for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
+       for (e=chart->edges; e; e=e->nextlink) {
         if (e->pair || (e->flag & PEDGE_DONE))
             continue;
 
-               (*nboundaries)++;
+               if (nboundaries)
+                       (*nboundaries)++;
+
         len = 0.0f;
     
                be = e;
@@ -966,13 +948,13 @@ static void p_chart_boundaries(PChart *chart, int *nboundaries, PEdge **outer)
                        be = be->next->vert->edge;
         } while(be != e);
 
-        if (len > maxlen) {
+        if (outer && (len > maxlen)) {
                        *outer = e;
             maxlen = len;
         }
     }
 
-       for (e=(PEdge*)chart->edges->first; e; e=e->link.next)
+       for (e=chart->edges; e; e=e->nextlink)
         e->flag &= ~PEDGE_DONE;
 }
 
@@ -1001,37 +983,20 @@ static float p_edge_boundary_angle(PEdge *e)
        return angle;
 }
 
-static PEdge *p_boundary_edge_next(PEdge *e)
-{
-       return e->next->vert->edge;
-}
-
-static PEdge *p_boundary_edge_prev(PEdge *e)
-{
-    PEdge *we = e, *last;
-
-       do {
-               last = we;
-               we = p_wheel_edge_next(we);
-       } while (we && (we != e));
-
-       return last->next->next;
-}
-
 static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
 {
        PEdge *e, *e1, *e2;
        PHashKey vkeys[3];
        PFace *f;
-       struct PHeap *heap = pheap_new(nedges);
+       struct Heap *heap = BLI_heap_new();
        float angle;
 
        e = be;
        do {
                angle = p_edge_boundary_angle(e);
-               e->u.heaplink = pheap_insert(heap, angle, e);
+               e->u.heaplink = BLI_heap_insert(heap, angle, e);
 
-               e = e->next->vert->edge;
+               e = p_boundary_edge_next(e);
        } while(e != be);
 
        if (nedges == 2) {
@@ -1040,30 +1005,30 @@ static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
                e->pair = be;
                be->pair = e;
 
-               pheap_remove(heap, e->u.heaplink);
-               pheap_remove(heap, be->u.heaplink);
+               BLI_heap_remove(heap, e->u.heaplink);
+               BLI_heap_remove(heap, be->u.heaplink);
        }
        else {
                while (nedges > 2) {
                        PEdge *ne, *ne1, *ne2;
 
-                       e = pheap_popmin(heap);
+                       e = (PEdge*)BLI_heap_popmin(heap);
 
                        e1 = p_boundary_edge_prev(e);
                        e2 = p_boundary_edge_next(e);
 
-                       pheap_remove(heap, e1->u.heaplink);
-                       pheap_remove(heap, e2->u.heaplink);
+                       BLI_heap_remove(heap, e1->u.heaplink);
+                       BLI_heap_remove(heap, e2->u.heaplink);
                        e->u.heaplink = e1->u.heaplink = e2->u.heaplink = NULL;
 
                        e->flag |= PEDGE_FILLED;
                        e1->flag |= PEDGE_FILLED;
 
-                       vkeys[0] = e->vert->link.key;
-                       vkeys[1] = e1->vert->link.key;
-                       vkeys[2] = e2->vert->link.key;
+                       vkeys[0] = e->vert->u.key;
+                       vkeys[1] = e1->vert->u.key;
+                       vkeys[2] = e2->vert->u.key;
 
-                       f = p_face_add(chart, -1, vkeys, NULL, NULL, 0, 1, 2, NULL, NULL);
+                       f = p_face_add_fill(chart, e->vert, e1->vert, e2->vert);
                        f->flag |= PFACE_FILLED;
 
                        ne = f->edge->next->next;
@@ -1088,15 +1053,15 @@ static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
                        else {
                                ne2->vert->edge = ne2;
                                
-                               ne2->u.heaplink = pheap_insert(heap, p_edge_boundary_angle(ne2), ne2);
-                               e2->u.heaplink = pheap_insert(heap, p_edge_boundary_angle(e2), e2);
+                               ne2->u.heaplink = BLI_heap_insert(heap, p_edge_boundary_angle(ne2), ne2);
+                               e2->u.heaplink = BLI_heap_insert(heap, p_edge_boundary_angle(e2), e2);
                        }
 
                        nedges--;
                }
        }
 
-       pheap_delete(heap);
+       BLI_heap_free(heap, NULL);
 }
 
 static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
@@ -1104,8 +1069,8 @@ static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
        PEdge *e, *enext, *be;
        int nedges;
 
-       for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
-               enext = e->link.next;
+       for (e=chart->edges; e; e=e->nextlink) {
+               enext = e->nextlink;
 
         if (e->pair || (e->flag & PEDGE_FILLED))
             continue;
@@ -1123,308 +1088,2901 @@ static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
     }
 }
 
-static void p_flush_uvs(PChart *chart)
+#if 0
+/* Polygon kernel for inserting uv's non overlapping */
+
+static int p_polygon_point_in(float *cp1, float *cp2, float *p)
 {
-       PEdge *e;
+       if ((cp1[0] == p[0]) && (cp1[1] == p[1]))
+               return 2;
+       else if ((cp2[0] == p[0]) && (cp2[1] == p[1]))
+               return 3;
+       else
+               return (p_area_signed(cp1, cp2, p) >= 0.0f);
+}
 
-       for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
-               if (e->orig_uv) {
-                       e->orig_uv[0] = e->vert->uv[0];
-                       e->orig_uv[1] = e->vert->uv[1];
+static void p_polygon_kernel_clip(float (*oldpoints)[2], int noldpoints, float (*newpoints)[2], int *nnewpoints, float *cp1, float *cp2)
+{
+       float *p2, *p1, isect[2];
+       int i, p2in, p1in;
+
+       p1 = oldpoints[noldpoints-1];
+       p1in = p_polygon_point_in(cp1, cp2, p1);
+       *nnewpoints = 0;
+
+       for (i = 0; i < noldpoints; i++) {
+               p2 = oldpoints[i];
+               p2in = p_polygon_point_in(cp1, cp2, p2);
+
+               if ((p2in >= 2) || (p1in && p2in)) {
+                       newpoints[*nnewpoints][0] = p2[0];
+                       newpoints[*nnewpoints][1] = p2[1];
+                       (*nnewpoints)++;
+               }
+               else if (p1in && !p2in) {
+                       if (p1in != 3) {
+                               p_intersect_line_2d(p1, p2, cp1, cp2, isect);
+                               newpoints[*nnewpoints][0] = isect[0];
+                               newpoints[*nnewpoints][1] = isect[1];
+                               (*nnewpoints)++;
+                       }
+               }
+               else if (!p1in && p2in) {
+                       p_intersect_line_2d(p1, p2, cp1, cp2, isect);
+                       newpoints[*nnewpoints][0] = isect[0];
+                       newpoints[*nnewpoints][1] = isect[1];
+                       (*nnewpoints)++;
+
+                       newpoints[*nnewpoints][0] = p2[0];
+                       newpoints[*nnewpoints][1] = p2[1];
+                       (*nnewpoints)++;
                }
+               
+               p1in = p2in;
+               p1 = p2;
        }
 }
 
-/* Exported */
-
-ParamHandle *param_construct_begin()
+static void p_polygon_kernel_center(float (*points)[2], int npoints, float *center)
 {
-       PHandle *handle = MEM_callocN(sizeof*handle, "PHandle");
-       handle->construction_chart = p_chart_new(handle);
-       handle->state = PHANDLE_STATE_ALLOCATED;
-       handle->arena = BLI_memarena_new((1<<16));
+       int i, size, nnewpoints = npoints;
+       float (*oldpoints)[2], (*newpoints)[2], *p1, *p2;
        
-       return (ParamHandle*)handle;
-}
-
-void param_delete(ParamHandle *handle)
-{
-       PHandle *phandle = (PHandle*)handle;
-       int i;
-
-       param_assert((phandle->state == PHANDLE_STATE_ALLOCATED) ||
-                    (phandle->state == PHANDLE_STATE_CONSTRUCTED));
+       size = npoints*3;
+       oldpoints = MEM_mallocN(sizeof(float)*2*size, "PPolygonOldPoints");
+       newpoints = MEM_mallocN(sizeof(float)*2*size, "PPolygonNewPoints");
 
-       for (i = 0; i < phandle->ncharts; i++)
-               p_chart_delete(phandle->charts[i]);
-       
-       if (phandle->charts)
-               MEM_freeN(phandle->charts);
+       memcpy(oldpoints, points, sizeof(float)*2*npoints);
 
-       if (phandle->construction_chart)
-               p_chart_delete(phandle->construction_chart);
+       for (i = 0; i < npoints; i++) {
+               p1 = points[i];
+               p2 = points[(i+1)%npoints];
+               p_polygon_kernel_clip(oldpoints, nnewpoints, newpoints, &nnewpoints, p1, p2);
 
-       BLI_memarena_free(phandle->arena);
-       MEM_freeN(phandle);
-}
+               if (nnewpoints == 0) {
+                       /* degenerate case, use center of original polygon */
+                       memcpy(oldpoints, points, sizeof(float)*2*npoints);
+                       nnewpoints = npoints;
+                       break;
+               }
+               else if (nnewpoints == 1) {
+                       /* degenerate case, use remaining point */
+                       center[0] = newpoints[0][0];
+                       center[1] = newpoints[0][1];
 
-void param_face_add(ParamHandle *handle, ParamKey key, int nverts,
-                    ParamKey *vkeys, float **co, float **uv,
-                    ParamBool *pin, ParamBool *select)
-{
-       PHandle *phandle = (PHandle*)handle;
-       PChart *chart = phandle->construction_chart;
+                       MEM_freeN(oldpoints);
+                       MEM_freeN(newpoints);
 
-       param_assert(phash_lookup(chart->faces, key) == NULL);
-       param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
-       param_assert((nverts == 3) || (nverts == 4));
+                       return;
+               }
 
-       if (nverts == 4) {
-               if (p_quad_split_direction(co)) {
-                       p_face_add(chart, key, vkeys, co, uv, 0, 1, 2, pin, select);
-                       p_face_add(chart, key, vkeys, co, uv, 0, 2, 3, pin, select);
+               if (nnewpoints*2 > size) {
+                       size *= 2;
+                       free(oldpoints);
+                       oldpoints = malloc(sizeof(float)*2*size);
+                       memcpy(oldpoints, newpoints, sizeof(float)*2*nnewpoints);
+                       free(newpoints);
+                       newpoints = malloc(sizeof(float)*2*size);
                }
                else {
-                       p_face_add(chart, key, vkeys, co, uv, 0, 1, 3, pin, select);
-                       p_face_add(chart, key, vkeys, co, uv, 1, 2, 3, pin, select);
+                       float (*sw_points)[2] = oldpoints;
+                       oldpoints = newpoints;
+                       newpoints = sw_points;
                }
        }
-       else
-               p_face_add(chart, key, vkeys, co, uv, 0, 1, 2, pin, select);
-}
 
-void param_edge_set_seam(ParamHandle *handle, ParamKey *vkeys)
-{
-       PHandle *phandle = (PHandle*)handle;
-       PChart *chart = phandle->construction_chart;
-       PEdge *e;
+       center[0] = center[1] = 0.0f;
 
-       param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
+       for (i = 0; i < nnewpoints; i++) {
+               center[0] += oldpoints[i][0];
+               center[1] += oldpoints[i][1];
+       }
 
-       e = p_edge_lookup(chart, vkeys);
-       if (e)
-               e->flag |= PEDGE_SEAM;
-}
+       center[0] /= nnewpoints;
+       center[1] /= nnewpoints;
 
-void param_construct_end(ParamHandle *handle, ParamBool fill, ParamBool impl)
-{
-       PHandle *phandle = (PHandle*)handle;
-       PChart *chart = phandle->construction_chart;
-       int i, j, nboundaries = 0;
-       PEdge *outer;
+       MEM_freeN(oldpoints);
+       MEM_freeN(newpoints);
+}
+#endif
 
-       param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
+#if 0
+/* Edge Collapser */
 
-       phandle->ncharts = p_connect_pairs(chart, impl);
-       phandle->charts = p_split_charts(phandle, chart, phandle->ncharts);
+int NCOLLAPSE = 1;
+int NCOLLAPSEX = 0;
+       
+static float p_vert_cotan(float *v1, float *v2, float *v3)
+{
+       float a[3], b[3], c[3], clen;
 
-       p_chart_delete(chart);
-       phandle->construction_chart = NULL;
+       VecSubf(a, v2, v1);
+       VecSubf(b, v3, v1);
+       Crossf(c, a, b);
 
-       for (i = j = 0; i < phandle->ncharts; i++) {
-               p_chart_boundaries(phandle->charts[i], &nboundaries, &outer);
+       clen = VecLength(c);
 
-               if (nboundaries > 0) {
-                       phandle->charts[j] = phandle->charts[i];
-                       j++;
+       if (clen == 0.0f)
+               return 0.0f;
+       
+       return Inpf(a, b)/clen;
+}
+       
+static PBool p_vert_flipped_wheel_triangle(PVert *v)
+{
+       PEdge *e = v->edge;
 
-                       if (fill && (nboundaries > 1))
-                               p_chart_fill_boundaries(phandle->charts[i], outer);
-               }
-               else
-                       p_chart_delete(phandle->charts[i]);
-       }
+       do {
+               if (p_face_uv_area_signed(e->face) < 0.0f)
+                       return P_TRUE;
 
-       phandle->ncharts = j;
+               e = p_wheel_edge_next(e);
+       } while (e && (e != v->edge));
 
-       phandle->state = PHANDLE_STATE_CONSTRUCTED;
+       return P_FALSE;
 }
 
-/* Least Squares Conformal Maps */
-
-static void p_chart_lscm_load_solution(PChart *chart)
+static PBool p_vert_map_harmonic_weights(PVert *v)
 {
-       PVert *pin = chart->u.lscm.singlepin, *v;
-       float translation[2];
+       float weightsum, positionsum[2], olduv[2];
 
-       if (pin) {
-               translation[0] = pin->uv[0] - nlGetVariable(2*pin->u.index);
-               translation[1] = pin->uv[1] - nlGetVariable(2*pin->u.index + 1);
-       }
-       else
-               translation[0] = translation[1] = 0.0f;
+       weightsum = 0.0f;
+       positionsum[0] = positionsum[1] = 0.0f;
+
+       if (p_vert_interior(v)) {
+               PEdge *e = v->edge;
 
-       for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
-               v->uv[0] = nlGetVariable(2*v->u.index) + translation[0];
-               v->uv[1] = nlGetVariable(2*v->u.index + 1) + translation[1];
+               do {
+                       float t1, t2, weight;
+                       PVert *v1, *v2;
+                       
+                       v1 = e->next->vert;
+                       v2 = e->next->next->vert;
+                       t1 = p_vert_cotan(v2->co, e->vert->co, v1->co);
+
+                       v1 = e->pair->next->vert;
+                       v2 = e->pair->next->next->vert;
+                       t2 = p_vert_cotan(v2->co, e->pair->vert->co, v1->co);
+
+                       weight = 0.5f*(t1 + t2);
+                       weightsum += weight;
+                       positionsum[0] += weight*e->pair->vert->uv[0];
+                       positionsum[1] += weight*e->pair->vert->uv[1];
+
+                       e = p_wheel_edge_next(e);
+               } while (e && (e != v->edge));
        }
-}
+       else {
+               PEdge *e = v->edge;
 
-static void p_chart_lscm_begin(PChart *chart)
-{
-       PVert *v, *pin1, *pin2;
-       int npins = 0, id = 0;
+               do {
+                       float t1, t2;
+                       PVert *v1, *v2;
 
-       nlNewContext();
-       nlSolverParameteri(NL_NB_VARIABLES, 2*phash_size(chart->verts));
-       nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
+                       v2 = e->next->vert;
+                       v1 = e->next->next->vert;
 
-       /* give vertices matrix indices and count pins */
-       for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
-               p_vert_load_pin_uvs(v);
+                       t1 = p_vert_cotan(v1->co, v->co, v2->co);
+                       t2 = p_vert_cotan(v2->co, v->co, v1->co);
 
-               if (v->flag & PVERT_PIN) {
-                       npins++;
-                       chart->u.lscm.singlepin = v;
-               }
+                       weightsum += t1 + t2;
+                       positionsum[0] += (v2->uv[1] - v1->uv[1]) + (t1*v2->uv[0] + t2*v1->uv[0]);
+                       positionsum[1] += (v1->uv[0] - v2->uv[0]) + (t1*v2->uv[1] + t2*v1->uv[1]);
+               
+                       e = p_wheel_edge_next(e);
+               } while (e && (e != v->edge));
+       }
 
-               v->u.index = id++;
+       if (weightsum != 0.0f) {
+               weightsum = 1.0f/weightsum;
+               positionsum[0] *= weightsum;
+               positionsum[1] *= weightsum;
        }
 
-       if (npins <= 1) {
-               /* not enough pins, lets find some ourself */
-               p_extrema_verts(chart, &pin1, &pin2);
+       olduv[0] = v->uv[0];
+       olduv[1] = v->uv[1];
+       v->uv[0] = positionsum[0];
+       v->uv[1] = positionsum[1];
+
+       if (p_vert_flipped_wheel_triangle(v)) {
+               v->uv[0] = olduv[0];
+               v->uv[1] = olduv[1];
 
-               chart->u.lscm.pin1 = pin1;
-               chart->u.lscm.pin2 = pin2;
+               return P_FALSE;
        }
-       else
-               chart->u.lscm.singlepin = NULL;
 
-       chart->u.lscm.context = nlGetCurrent();
+       return P_TRUE;
 }
 
-static PBool p_chart_lscm_solve(PChart *chart)
+static void p_vert_harmonic_insert(PVert *v)
 {
-       PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
-       PFace *f;
+       PEdge *e;
 
-       nlMakeCurrent(chart->u.lscm.context);
+       if (!p_vert_map_harmonic_weights(v)) {
+               /* do polygon kernel center insertion: this is quite slow, but should
+                  only be needed for 0.01 % of verts or so, when insert with harmonic
+                  weights fails */
 
-       nlBegin(NL_SYSTEM);
+               int npoints = 0, i;
+               float (*points)[2];
 
-       if (chart->u.lscm.pin1) {
-               nlLockVariable(2*pin1->u.index);
-               nlLockVariable(2*pin1->u.index + 1);
-               nlLockVariable(2*pin2->u.index);
-               nlLockVariable(2*pin2->u.index + 1);
-       
-               nlSetVariable(2*pin1->u.index, 0.0f);
-               nlSetVariable(2*pin1->u.index + 1, 0.5f);
-               nlSetVariable(2*pin2->u.index, 1.0f);
-               nlSetVariable(2*pin2->u.index + 1, 0.5f);
-       }
-       else {
-               /* set and lock the pins */
-               for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
-                       if (v->flag & PVERT_PIN) {
-                               nlLockVariable(2*v->u.index);
-                               nlLockVariable(2*v->u.index + 1);
+               e = v->edge;
+               do {
+                       npoints++;      
+                       e = p_wheel_edge_next(e);
+               } while (e && (e != v->edge));
 
-                               nlSetVariable(2*v->u.index, v->uv[0]);
-                               nlSetVariable(2*v->u.index + 1, v->uv[1]);
-                       }
-               }
-       }
+               if (e == NULL)
+                       npoints++;
 
-       /* construct matrix */
+               points = MEM_mallocN(sizeof(float)*2*npoints, "PHarmonicPoints");
 
-       nlBegin(NL_MATRIX);
+               e = v->edge;
+               i = 0;
+               do {
+                       PEdge *nexte = p_wheel_edge_next(e);
 
-       for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
-               PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-               PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-               float a1, a2, a3, ratio, cosine, sine;
-               float sina1, sina2, sina3, sinmax;
+                       points[i][0] = e->next->vert->uv[0]; 
+                       points[i][1] = e->next->vert->uv[1]; 
 
-               if (chart->u.lscm.abf_alpha) {
-                       /* use abf angles if passed on */
-                       a1 = *(chart->u.lscm.abf_alpha++);
-                       a2 = *(chart->u.lscm.abf_alpha++);
-                       a3 = *(chart->u.lscm.abf_alpha++);
-               }
-               else
-                       p_face_angles(f, &a1, &a2, &a3);
+                       if (nexte == NULL) {
+                               i++;
+                               points[i][0] = e->next->next->vert->uv[0]; 
+                               points[i][1] = e->next->next->vert->uv[1]; 
+                               break;
+                       }
 
-               sina1 = sin(a1);
-               sina2 = sin(a2);
-               sina3 = sin(a3);
+                       e = nexte;
+                       i++;
+               } while (e != v->edge);
+               
+               p_polygon_kernel_center(points, npoints, v->uv);
+
+               MEM_freeN(points);
+       }
+
+       e = v->edge;
+       do {
+               if (!(e->next->vert->flag & PVERT_PIN))
+                       p_vert_map_harmonic_weights(e->next->vert);
+               e = p_wheel_edge_next(e);
+       } while (e && (e != v->edge));
+
+       p_vert_map_harmonic_weights(v);
+}
+
+static void p_vert_fix_edge_pointer(PVert *v)
+{
+       PEdge *start = v->edge;
+
+       /* set v->edge pointer to the edge with no pair, if there is one */
+       while (v->edge->pair) {
+               v->edge = p_wheel_edge_prev(v->edge);
+               
+               if (v->edge == start)
+                       break;
+       }
+}
+
+static void p_collapsing_verts(PEdge *edge, PEdge *pair, PVert **newv, PVert **keepv)
+{
+       /* the two vertices that are involved in the collapse */
+       if (edge) {
+               *newv = edge->vert;
+               *keepv = edge->next->vert;
+       }
+       else {
+               *newv = pair->next->vert;
+               *keepv = pair->vert;
+       }
+}
+
+static void p_collapse_edge(PEdge *edge, PEdge *pair)
+{
+       PVert *oldv, *keepv;
+       PEdge *e;
+
+       p_collapsing_verts(edge, pair, &oldv, &keepv);
+
+       /* change e->vert pointers from old vertex to the target vertex */
+       e = oldv->edge;
+       do {
+               if ((e != edge) && !(pair && pair->next == e))
+                       e->vert = keepv;
+
+               e = p_wheel_edge_next(e);
+       } while (e && (e != oldv->edge));
+
+       /* set keepv->edge pointer */
+       if ((edge && (keepv->edge == edge->next)) || (keepv->edge == pair)) {
+               if (edge && edge->next->pair)
+                       keepv->edge = edge->next->pair->next;
+               else if (pair && pair->next->next->pair)
+                       keepv->edge = pair->next->next->pair;
+               else if (edge && edge->next->next->pair)
+                       keepv->edge = edge->next->next->pair;
+               else
+                       keepv->edge = pair->next->pair->next;
+       }
+       
+       /* update pairs and v->edge pointers */
+       if (edge) {
+               PEdge *e1 = edge->next, *e2 = e1->next;
+
+               if (e1->pair)
+                       e1->pair->pair = e2->pair;
+
+               if (e2->pair) {
+                       e2->pair->pair = e1->pair;
+                       e2->vert->edge = p_wheel_edge_prev(e2);
+               }
+               else
+                       e2->vert->edge = p_wheel_edge_next(e2);
+
+               p_vert_fix_edge_pointer(e2->vert);
+       }
+
+       if (pair) {
+               PEdge *e1 = pair->next, *e2 = e1->next;
+
+               if (e1->pair)
+                       e1->pair->pair = e2->pair;
+
+               if (e2->pair) {
+                       e2->pair->pair = e1->pair;
+                       e2->vert->edge = p_wheel_edge_prev(e2);
+               }
+               else
+                       e2->vert->edge = p_wheel_edge_next(e2);
+
+               p_vert_fix_edge_pointer(e2->vert);
+       }
+
+       p_vert_fix_edge_pointer(keepv);
+
+       /* mark for move to collapsed list later */
+       oldv->flag |= PVERT_COLLAPSE;
+
+       if (edge) {
+               PFace *f = edge->face;
+               PEdge *e1 = edge->next, *e2 = e1->next;
+
+               f->flag |= PFACE_COLLAPSE;
+               edge->flag |= PEDGE_COLLAPSE;
+               e1->flag |= PEDGE_COLLAPSE;
+               e2->flag |= PEDGE_COLLAPSE;
+       }
+
+       if (pair) {
+               PFace *f = pair->face;
+               PEdge *e1 = pair->next, *e2 = e1->next;
+
+               f->flag |= PFACE_COLLAPSE;
+               pair->flag |= PEDGE_COLLAPSE;
+               e1->flag |= PEDGE_COLLAPSE;
+               e2->flag |= PEDGE_COLLAPSE;
+       }
+}
+
+static void p_split_vertex(PEdge *edge, PEdge *pair)
+{
+       PVert *newv, *keepv;
+       PEdge *e;
+
+       p_collapsing_verts(edge, pair, &newv, &keepv);
+
+       /* update edge pairs */
+       if (edge) {
+               PEdge *e1 = edge->next, *e2 = e1->next;
+
+               if (e1->pair)
+                       e1->pair->pair = e1;
+               if (e2->pair)
+                       e2->pair->pair = e2;
+
+               e2->vert->edge = e2;
+               p_vert_fix_edge_pointer(e2->vert);
+               keepv->edge = e1;
+       }
+
+       if (pair) {
+               PEdge *e1 = pair->next, *e2 = e1->next;
+
+               if (e1->pair)
+                       e1->pair->pair = e1;
+               if (e2->pair)
+                       e2->pair->pair = e2;
+
+               e2->vert->edge = e2;
+               p_vert_fix_edge_pointer(e2->vert);
+               keepv->edge = pair;
+       }
+
+       p_vert_fix_edge_pointer(keepv);
+
+       /* set e->vert pointers to restored vertex */
+       e = newv->edge;
+       do {
+               e->vert = newv;
+               e = p_wheel_edge_next(e);
+       } while (e && (e != newv->edge));
+}
+
+static PBool p_collapse_allowed_topologic(PEdge *edge, PEdge *pair)
+{
+       PVert *oldv, *keepv;
+
+       p_collapsing_verts(edge, pair, &oldv, &keepv);
+
+       /* boundary edges */
+       if (!edge || !pair) {
+               /* avoid collapsing chart into an edge */
+               if (edge && !edge->next->pair && !edge->next->next->pair)
+                       return P_FALSE;
+               else if (pair && !pair->next->pair && !pair->next->next->pair)
+                       return P_FALSE;
+       }
+       /* avoid merging two boundaries (oldv and keepv are on the 'other side' of
+          the chart) */
+       else if (!p_vert_interior(oldv) && !p_vert_interior(keepv))
+               return P_FALSE;
+       
+       return P_TRUE;
+}
+
+static PBool p_collapse_normal_flipped(float *v1, float *v2, float *vold, float *vnew)
+{
+       float nold[3], nnew[3], sub1[3], sub2[3];
+
+       VecSubf(sub1, vold, v1);
+       VecSubf(sub2, vold, v2);
+       Crossf(nold, sub1, sub2);
+
+       VecSubf(sub1, vnew, v1);
+       VecSubf(sub2, vnew, v2);
+       Crossf(nnew, sub1, sub2);
+
+       return (Inpf(nold, nnew) <= 0.0f);
+}
+
+static PBool p_collapse_allowed_geometric(PEdge *edge, PEdge *pair)
+{
+       PVert *oldv, *keepv;
+       PEdge *e;
+       float angulardefect, angle;
+
+       p_collapsing_verts(edge, pair, &oldv, &keepv);
+
+       angulardefect = 2*M_PI;
+
+       e = oldv->edge;
+       do {
+               float a[3], b[3], minangle, maxangle;
+               PEdge *e1 = e->next, *e2 = e1->next;
+               PVert *v1 = e1->vert, *v2 = e2->vert;
+               int i;
+
+               angle = p_vec_angle(v1->co, oldv->co, v2->co);
+               angulardefect -= angle;
+
+               /* skip collapsing faces */
+               if (v1 == keepv || v2 == keepv) {
+                       e = p_wheel_edge_next(e);
+                       continue;
+               }
+
+               if (p_collapse_normal_flipped(v1->co, v2->co, oldv->co, keepv->co))
+                       return P_FALSE;
+       
+               a[0] = angle;
+               a[1] = p_vec_angle(v2->co, v1->co, oldv->co);
+               a[2] = M_PI - a[0] - a[1];
+
+               b[0] = p_vec_angle(v1->co, keepv->co, v2->co);
+               b[1] = p_vec_angle(v2->co, v1->co, keepv->co);
+               b[2] = M_PI - b[0] - b[1];
+
+               /* abf criterion 1: avoid sharp and obtuse angles */
+               minangle = 15.0f*M_PI/180.0f;
+               maxangle = M_PI - minangle;
+
+               for (i = 0; i < 3; i++) {
+                       if ((b[i] < a[i]) && (b[i] < minangle))
+                               return P_FALSE;
+                       else if ((b[i] > a[i]) && (b[i] > maxangle))
+                               return P_FALSE;
+               }
+
+               e = p_wheel_edge_next(e);
+       } while (e && (e != oldv->edge));
+
+       if (p_vert_interior(oldv)) {
+               /* hlscm criterion: angular defect smaller than threshold */
+               if (fabs(angulardefect) > (M_PI*30.0/180.0))
+                       return P_FALSE;
+       }
+       else {
+               PVert *v1 = p_boundary_edge_next(oldv->edge)->vert;
+               PVert *v2 = p_boundary_edge_prev(oldv->edge)->vert;
+
+               /* abf++ criterion 2: avoid collapsing verts inwards */
+               if (p_vert_interior(keepv))
+                       return P_FALSE;
+               
+               /* don't collapse significant boundary changes */
+               angle = p_vec_angle(v1->co, oldv->co, v2->co);
+               if (angle < (M_PI*160.0/180.0))
+                       return P_FALSE;
+       }
+
+       return P_TRUE;
+}
+
+static PBool p_collapse_allowed(PEdge *edge, PEdge *pair)
+{
+       PVert *oldv, *keepv;
+
+       p_collapsing_verts(edge, pair, &oldv, &keepv);
+
+       if (oldv->flag & PVERT_PIN)
+               return P_FALSE;
+       
+       return (p_collapse_allowed_topologic(edge, pair) &&
+               p_collapse_allowed_geometric(edge, pair));
+}
+
+static float p_collapse_cost(PEdge *edge, PEdge *pair)
+{
+       /* based on volume and boundary optimization from:
+         "Fast and Memory Efficient Polygonal Simplification" P. Lindstrom, G. Turk */
+
+       PVert *oldv, *keepv;
+       PEdge *e;
+       PFace *oldf1, *oldf2;
+       float volumecost = 0.0f, areacost = 0.0f, edgevec[3], cost, weight, elen;
+       float shapecost = 0.0f;
+       float shapeold = 0.0f, shapenew = 0.0f;
+       int nshapeold = 0, nshapenew = 0;
+
+       p_collapsing_verts(edge, pair, &oldv, &keepv);
+       oldf1 = (edge)? edge->face: NULL;
+       oldf2 = (pair)? pair->face: NULL;
+
+       VecSubf(edgevec, keepv->co, oldv->co);
+
+       e = oldv->edge;
+       do {
+               float a1, a2, a3;
+               float *co1 = e->next->vert->co;
+               float *co2 = e->next->next->vert->co;
+
+               if ((e->face != oldf1) && (e->face != oldf2)) {
+                       float tetrav2[3], tetrav3[3], c[3];
+
+                       /* tetrahedron volume = (1/3!)*|a.(b x c)| */
+                       VecSubf(tetrav2, co1, oldv->co);
+                       VecSubf(tetrav3, co2, oldv->co);
+                       Crossf(c, tetrav2, tetrav3);
+
+                       volumecost += fabs(Inpf(edgevec, c)/6.0f);
+#if 0
+                       shapecost += Inpf(co1, keepv->co);
+
+                       if (p_wheel_edge_next(e) == NULL)
+                               shapecost += Inpf(co2, keepv->co);
+#endif
+
+                       p_triangle_angles(oldv->co, co1, co2, &a1, &a2, &a3);
+                       a1 = a1 - M_PI/3.0;
+                       a2 = a2 - M_PI/3.0;
+                       a3 = a3 - M_PI/3.0;
+                       shapeold = (a1*a1 + a2*a2 + a3*a3)/((M_PI/2)*(M_PI/2));
+
+                       nshapeold++;
+               }
+               else {
+                       p_triangle_angles(keepv->co, co1, co2, &a1, &a2, &a3);
+                       a1 = a1 - M_PI/3.0;
+                       a2 = a2 - M_PI/3.0;
+                       a3 = a3 - M_PI/3.0;
+                       shapenew = (a1*a1 + a2*a2 + a3*a3)/((M_PI/2)*(M_PI/2));
+
+                       nshapenew++;
+               }
+
+               e = p_wheel_edge_next(e);
+       } while (e && (e != oldv->edge));
+
+       if (!p_vert_interior(oldv)) {
+               PVert *v1 = p_boundary_edge_prev(oldv->edge)->vert;
+               PVert *v2 = p_boundary_edge_next(oldv->edge)->vert;
+
+               areacost = AreaT3Dfl(oldv->co, v1->co, v2->co);
+       }
+
+       elen = VecLength(edgevec);
+       weight = 1.0f; /* 0.2f */
+       cost = weight*volumecost*volumecost + elen*elen*areacost*areacost;
+#if 0
+       cost += shapecost;
+#else
+       shapeold /= nshapeold;
+       shapenew /= nshapenew;
+       shapecost = (shapeold + 0.00001)/(shapenew + 0.00001);
+
+       cost *= shapecost;
+#endif
+
+       return cost;
+}
+       
+static void p_collapse_cost_vertex(PVert *vert, float *mincost, PEdge **mine)
+{
+       PEdge *e, *enext, *pair;
+
+       *mine = NULL;
+       *mincost = 0.0f;
+       e = vert->edge;
+       do {
+               if (p_collapse_allowed(e, e->pair)) {
+                       float cost = p_collapse_cost(e, e->pair);
+
+                       if ((*mine == NULL) || (cost < *mincost)) {
+                               *mincost = cost;
+                               *mine = e;
+                       }
+               }
+
+               enext = p_wheel_edge_next(e);
+
+               if (enext == NULL) {
+                       /* the other boundary edge, where we only have the pair halfedge */
+                       pair = e->next->next;
+
+                       if (p_collapse_allowed(NULL, pair)) {
+                               float cost = p_collapse_cost(NULL, pair);
+
+                               if ((*mine == NULL) || (cost < *mincost)) {
+                                       *mincost = cost;
+                                       *mine = pair;
+                               }
+                       }
+
+                       break;
+               }
+
+               e = enext;
+       } while (e != vert->edge);
+}
+
+static void p_chart_post_collapse_flush(PChart *chart, PEdge *collapsed)
+{
+       /* move to collapsed_ */
+
+       PVert *v, *nextv = NULL, *verts = chart->verts;
+       PEdge *e, *nexte = NULL, *edges = chart->edges, *laste = NULL;
+       PFace *f, *nextf = NULL, *faces = chart->faces;
+
+       chart->verts = chart->collapsed_verts = NULL;
+       chart->edges = chart->collapsed_edges = NULL;
+       chart->faces = chart->collapsed_faces = NULL;
+
+       chart->nverts = chart->nedges = chart->nfaces = 0;
+
+       for (v=verts; v; v=nextv) {
+               nextv = v->nextlink;
+
+               if (v->flag & PVERT_COLLAPSE) {
+                       v->nextlink = chart->collapsed_verts;
+                       chart->collapsed_verts = v;
+               }
+               else {
+                       v->nextlink = chart->verts;
+                       chart->verts = v;
+                       chart->nverts++;
+               }
+       }
+
+       for (e=edges; e; e=nexte) {
+               nexte = e->nextlink;
+
+               if (!collapsed || !(e->flag & PEDGE_COLLAPSE_EDGE)) {
+                       if (e->flag & PEDGE_COLLAPSE) {
+                               e->nextlink = chart->collapsed_edges;
+                               chart->collapsed_edges = e;
+                       }
+                       else {
+                               e->nextlink = chart->edges;
+                               chart->edges = e;
+                               chart->nedges++;
+                       }
+               }
+       }
+
+       /* these are added last so they can be popped of in the right order
+          for splitting */
+       for (e=collapsed; e; e=e->nextlink) {
+               e->nextlink = e->u.nextcollapse;
+               laste = e;
+       }
+       if (laste) {
+               laste->nextlink = chart->collapsed_edges;
+               chart->collapsed_edges = collapsed;
+       }
+
+       for (f=faces; f; f=nextf) {
+               nextf = f->nextlink;
+
+               if (f->flag & PFACE_COLLAPSE) {
+                       f->nextlink = chart->collapsed_faces;
+                       chart->collapsed_faces = f;
+               }
+               else {
+                       f->nextlink = chart->faces;
+                       chart->faces = f;
+                       chart->nfaces++;
+               }
+       }
+}
+
+static void p_chart_post_split_flush(PChart *chart)
+{
+       /* move from collapsed_ */
+
+       PVert *v, *nextv = NULL;
+       PEdge *e, *nexte = NULL;
+       PFace *f, *nextf = NULL;
+
+       for (v=chart->collapsed_verts; v; v=nextv) {
+               nextv = v->nextlink;
+               v->nextlink = chart->verts;
+               chart->verts = v;
+               chart->nverts++;
+       }
+
+       for (e=chart->collapsed_edges; e; e=nexte) {
+               nexte = e->nextlink;
+               e->nextlink = chart->edges;
+               chart->edges = e;
+               chart->nedges++;
+       }
+
+       for (f=chart->collapsed_faces; f; f=nextf) {
+               nextf = f->nextlink;
+               f->nextlink = chart->faces;
+               chart->faces = f;
+               chart->nfaces++;
+       }
+
+       chart->collapsed_verts = NULL;
+       chart->collapsed_edges = NULL;
+       chart->collapsed_faces = NULL;
+}
+
+static void p_chart_simplify_compute(PChart *chart)
+{
+       /* Computes a list of edge collapses / vertex splits. The collapsed
+          simplices go in the chart->collapsed_* lists, The original and
+          collapsed may then be view as stacks, where the next collapse/split
+          is at the top of the respective lists. */
+
+       Heap *heap = BLI_heap_new();
+       PVert *v, **wheelverts;
+       PEdge *collapsededges = NULL, *e;
+       int nwheelverts, i, ncollapsed = 0;
+
+       wheelverts = MEM_mallocN(sizeof(PVert*)*chart->nverts, "PChartWheelVerts");
+
+       /* insert all potential collapses into heap */
+       for (v=chart->verts; v; v=v->nextlink) {
+               float cost;
+               PEdge *e = NULL;
+               
+               p_collapse_cost_vertex(v, &cost, &e);
+
+               if (e)
+                       v->u.heaplink = BLI_heap_insert(heap, cost, e);
+               else
+                       v->u.heaplink = NULL;
+       }
+
+       for (e=chart->edges; e; e=e->nextlink)
+               e->u.nextcollapse = NULL;
+
+       /* pop edge collapse out of heap one by one */
+       while (!BLI_heap_empty(heap)) {
+               if (ncollapsed == NCOLLAPSE)
+                       break;
+
+               HeapNode *link = BLI_heap_top(heap);
+               PEdge *edge = (PEdge*)BLI_heap_popmin(heap), *pair = edge->pair;
+               PVert *oldv, *keepv;
+               PEdge *wheele, *nexte;
+
+               /* remember the edges we collapsed */
+               edge->u.nextcollapse = collapsededges;
+               collapsededges = edge;
+
+               if (edge->vert->u.heaplink != link) {
+                       edge->flag |= (PEDGE_COLLAPSE_EDGE|PEDGE_COLLAPSE_PAIR);
+                       edge->next->vert->u.heaplink = NULL;
+                       SWAP(PEdge*, edge, pair);
+               }
+               else {
+                       edge->flag |= PEDGE_COLLAPSE_EDGE;
+                       edge->vert->u.heaplink = NULL;
+               }
+
+               p_collapsing_verts(edge, pair, &oldv, &keepv);
+
+               /* gather all wheel verts and remember them before collapse */
+               nwheelverts = 0;
+               wheele = oldv->edge;
+
+               do {
+                       wheelverts[nwheelverts++] = wheele->next->vert;
+                       nexte = p_wheel_edge_next(wheele);
+
+                       if (nexte == NULL)
+                               wheelverts[nwheelverts++] = wheele->next->next->vert;
+
+                       wheele = nexte;
+               } while (wheele && (wheele != oldv->edge));
+
+               /* collapse */
+               p_collapse_edge(edge, pair);
+
+               for (i = 0; i < nwheelverts; i++) {
+                       float cost;
+                       PEdge *collapse = NULL;
+
+                       v = wheelverts[i];
+
+                       if (v->u.heaplink) {
+                               BLI_heap_remove(heap, v->u.heaplink);
+                               v->u.heaplink = NULL;
+                       }
+               
+                       p_collapse_cost_vertex(v, &cost, &collapse);
+
+                       if (collapse)
+                               v->u.heaplink = BLI_heap_insert(heap, cost, collapse);
+               }
+
+               ncollapsed++;
+       }
+
+       MEM_freeN(wheelverts);
+       BLI_heap_free(heap, NULL);
+
+       p_chart_post_collapse_flush(chart, collapsededges);
+}
+
+static void p_chart_complexify(PChart *chart)
+{
+       PEdge *e, *pair, *edge;
+       PVert *newv, *keepv;
+       int x = 0;
+
+       for (e=chart->collapsed_edges; e; e=e->nextlink) {
+               if (!(e->flag & PEDGE_COLLAPSE_EDGE))
+                       break;
+
+               edge = e;
+               pair = e->pair;
+
+               if (edge->flag & PEDGE_COLLAPSE_PAIR) {
+                       SWAP(PEdge*, edge, pair);
+               }
+
+               p_split_vertex(edge, pair);
+               p_collapsing_verts(edge, pair, &newv, &keepv);
+
+               if (x >= NCOLLAPSEX) {
+                       newv->uv[0] = keepv->uv[0];
+                       newv->uv[1] = keepv->uv[1];
+               }
+               else {
+                       p_vert_harmonic_insert(newv);
+                       x++;
+               }
+       }
+
+       p_chart_post_split_flush(chart);
+}
+
+#if 0
+static void p_chart_simplify(PChart *chart)
+{
+       /* Not implemented, needs proper reordering in split_flush. */
+}
+#endif
+#endif
+
+/* ABF */
+
+#define ABF_MAX_ITER 20
+
+typedef struct PAbfSystem {
+       int ninterior, nfaces, nangles;
+       float *alpha, *beta, *sine, *cosine, *weight;
+       float *bAlpha, *bTriangle, *bInterior;
+       float *lambdaTriangle, *lambdaPlanar, *lambdaLength;
+       float (*J2dt)[3], *bstar, *dstar;
+       float minangle, maxangle;
+} PAbfSystem;
+
+static void p_abf_setup_system(PAbfSystem *sys)
+{
+       int i;
+
+       sys->alpha = (float*)MEM_mallocN(sizeof(float)*sys->nangles, "ABFalpha");
+       sys->beta = (float*)MEM_mallocN(sizeof(float)*sys->nangles, "ABFbeta");
+       sys->sine = (float*)MEM_mallocN(sizeof(float)*sys->nangles, "ABFsine");
+       sys->cosine = (float*)MEM_mallocN(sizeof(float)*sys->nangles, "ABFcosine");
+       sys->weight = (float*)MEM_mallocN(sizeof(float)*sys->nangles, "ABFweight");
+
+       sys->bAlpha = (float*)MEM_mallocN(sizeof(float)*sys->nangles, "ABFbalpha");
+       sys->bTriangle = (float*)MEM_mallocN(sizeof(float)*sys->nfaces, "ABFbtriangle");
+       sys->bInterior = (float*)MEM_mallocN(sizeof(float)*2*sys->ninterior, "ABFbinterior");
+
+       sys->lambdaTriangle = (float*)MEM_callocN(sizeof(float)*sys->nfaces, "ABFlambdatri");
+       sys->lambdaPlanar = (float*)MEM_callocN(sizeof(float)*sys->ninterior, "ABFlamdaplane");
+       sys->lambdaLength = (float*)MEM_mallocN(sizeof(float)*sys->ninterior, "ABFlambdalen");
+
+       sys->J2dt = MEM_mallocN(sizeof(float)*sys->nangles*3, "ABFj2dt");
+       sys->bstar = (float*)MEM_mallocN(sizeof(float)*sys->nfaces, "ABFbstar");
+       sys->dstar = (float*)MEM_mallocN(sizeof(float)*sys->nfaces, "ABFdstar");
+
+       for (i = 0; i < sys->ninterior; i++)
+               sys->lambdaLength[i] = 1.0;
+       
+       sys->minangle = 7.5f*M_PI/180.0f;
+       sys->maxangle = M_PI - sys->minangle;
+}
+
+static void p_abf_free_system(PAbfSystem *sys)
+{
+       MEM_freeN(sys->alpha);
+       MEM_freeN(sys->beta);
+       MEM_freeN(sys->sine);
+       MEM_freeN(sys->cosine);
+       MEM_freeN(sys->weight);
+       MEM_freeN(sys->bAlpha);
+       MEM_freeN(sys->bTriangle);
+       MEM_freeN(sys->bInterior);
+       MEM_freeN(sys->lambdaTriangle);
+       MEM_freeN(sys->lambdaPlanar);
+       MEM_freeN(sys->lambdaLength);
+       MEM_freeN(sys->J2dt);
+       MEM_freeN(sys->bstar);
+       MEM_freeN(sys->dstar);
+}
+
+static void p_abf_compute_sines(PAbfSystem *sys)
+{
+       int i;
+       float *sine = sys->sine, *cosine = sys->cosine, *alpha = sys->alpha;
+
+       for (i = 0; i < sys->nangles; i++, sine++, cosine++, alpha++) {
+               *sine = sin(*alpha);
+               *cosine = cos(*alpha);
+       }
+}
+
+static float p_abf_compute_sin_product(PAbfSystem *sys, PVert *v, int aid)
+{
+       PEdge *e, *e1, *e2;
+       float sin1, sin2;
+
+       sin1 = sin2 = 1.0;
+
+       e = v->edge;
+       do {
+               e1 = e->next;
+               e2 = e->next->next;
+
+               if (aid == e1->u.id) {
+                       /* we are computing a derivative for this angle,
+                          so we use cos and drop the other part */
+                       sin1 *= sys->cosine[e1->u.id];
+                       sin2 = 0.0;
+               }
+               else
+                       sin1 *= sys->sine[e1->u.id];
+
+               if (aid == e2->u.id) {
+                       /* see above */
+                       sin1 = 0.0;
+                       sin2 *= sys->cosine[e2->u.id];
+               }
+               else
+                       sin2 *= sys->sine[e2->u.id];
+
+               e = e->next->next->pair;
+       } while (e && (e != v->edge));
+
+       return (sin1 - sin2);
+}
+
+static float p_abf_compute_grad_alpha(PAbfSystem *sys, PFace *f, PEdge *e)
+{
+       PVert *v = e->vert, *v1 = e->next->vert, *v2 = e->next->next->vert;
+       float deriv;
+
+       deriv = (sys->alpha[e->u.id] - sys->beta[e->u.id])*sys->weight[e->u.id];
+       deriv += sys->lambdaTriangle[f->u.id];
+
+       if (v->flag & PVERT_INTERIOR) {
+               deriv += sys->lambdaPlanar[v->u.id];
+       }
+
+       if (v1->flag & PVERT_INTERIOR) {
+               float product = p_abf_compute_sin_product(sys, v1, e->u.id);
+               deriv += sys->lambdaLength[v1->u.id]*product;
+       }
+
+       if (v2->flag & PVERT_INTERIOR) {
+               float product = p_abf_compute_sin_product(sys, v2, e->u.id);
+               deriv += sys->lambdaLength[v2->u.id]*product;
+       }
+
+       return deriv;
+}
+
+static float p_abf_compute_gradient(PAbfSystem *sys, PChart *chart)
+{
+       PFace *f;
+       PEdge *e;
+       PVert *v;
+       float norm = 0.0;
+
+       for (f=chart->faces; f; f=f->nextlink) {
+               PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+               float gtriangle, galpha1, galpha2, galpha3;
+
+               galpha1 = p_abf_compute_grad_alpha(sys, f, e1);
+               galpha2 = p_abf_compute_grad_alpha(sys, f, e2);
+               galpha3 = p_abf_compute_grad_alpha(sys, f, e3);
+
+               sys->bAlpha[e1->u.id] = -galpha1;
+               sys->bAlpha[e2->u.id] = -galpha2;
+               sys->bAlpha[e3->u.id] = -galpha3;
+
+               norm += galpha1*galpha1 + galpha2*galpha2 + galpha3*galpha3;
+
+               gtriangle = sys->alpha[e1->u.id] + sys->alpha[e2->u.id] + sys->alpha[e3->u.id] - M_PI;
+               sys->bTriangle[f->u.id] = -gtriangle;
+               norm += gtriangle*gtriangle;
+       }
+
+       for (v=chart->verts; v; v=v->nextlink) {
+               if (v->flag & PVERT_INTERIOR) {
+                       float gplanar = -2*M_PI, glength;
+
+                       e = v->edge;
+                       do {
+                               gplanar += sys->alpha[e->u.id];
+                               e = e->next->next->pair;
+                       } while (e && (e != v->edge));
+
+                       sys->bInterior[v->u.id] = -gplanar;
+                       norm += gplanar*gplanar;
+
+                       glength = p_abf_compute_sin_product(sys, v, -1);
+                       sys->bInterior[sys->ninterior + v->u.id] = -glength;
+                       norm += glength*glength;
+               }
+       }
+
+       return norm;
+}
+
+static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
+{
+       PFace *f;
+       PEdge *e;
+       int i, j, ninterior = sys->ninterior, nvar = 2*sys->ninterior;
+       PBool success;
+
+       nlNewContext();
+       nlSolverParameteri(NL_NB_VARIABLES, nvar);
+
+       nlBegin(NL_SYSTEM);
+
+       nlBegin(NL_MATRIX);
+
+       for (i = 0; i < nvar; i++)
+               nlRightHandSideAdd(0, i, sys->bInterior[i]);
+
+       for (f=chart->faces; f; f=f->nextlink) {
+               float wi1, wi2, wi3, b, si, beta[3], j2[3][3], W[3][3];
+               float row1[6], row2[6], row3[6];
+               int vid[6];
+               PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+               PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+
+               wi1 = 1.0/sys->weight[e1->u.id];
+               wi2 = 1.0/sys->weight[e2->u.id];
+               wi3 = 1.0/sys->weight[e3->u.id];
+
+               /* bstar1 = (J1*dInv*bAlpha - bTriangle) */
+               b = sys->bAlpha[e1->u.id]*wi1;
+               b += sys->bAlpha[e2->u.id]*wi2;
+               b += sys->bAlpha[e3->u.id]*wi3;
+               b -= sys->bTriangle[f->u.id];
+
+               /* si = J1*d*J1t */
+               si = 1.0/(wi1 + wi2 + wi3);
+
+               /* J1t*si*bstar1 - bAlpha */
+               beta[0] = b*si - sys->bAlpha[e1->u.id];
+               beta[1] = b*si - sys->bAlpha[e2->u.id];
+               beta[2] = b*si - sys->bAlpha[e3->u.id];
+
+               /* use this later for computing other lambda's */
+               sys->bstar[f->u.id] = b;
+               sys->dstar[f->u.id] = si;
+
+               /* set matrix */
+               W[0][0] = si - sys->weight[e1->u.id]; W[0][1] = si; W[0][2] = si;
+               W[1][0] = si; W[1][1] = si - sys->weight[e2->u.id]; W[1][2] = si;
+               W[2][0] = si; W[2][1] = si; W[2][2] = si - sys->weight[e3->u.id];
+
+               vid[0] = vid[1] = vid[2] = vid[3] = vid[4] = vid[5] = -1;
+
+               if (v1->flag & PVERT_INTERIOR) {
+                       vid[0] = v1->u.id;
+                       vid[3] = ninterior + v1->u.id;
+
+                       sys->J2dt[e1->u.id][0] = j2[0][0] = 1.0*wi1;
+                       sys->J2dt[e2->u.id][0] = j2[1][0] = p_abf_compute_sin_product(sys, v1, e2->u.id)*wi2;
+                       sys->J2dt[e3->u.id][0] = j2[2][0] = p_abf_compute_sin_product(sys, v1, e3->u.id)*wi3;
+
+                       nlRightHandSideAdd(0, v1->u.id, j2[0][0]*beta[0]);
+                       nlRightHandSideAdd(0, ninterior + v1->u.id, j2[1][0]*beta[1] + j2[2][0]*beta[2]);
+
+                       row1[0] = j2[0][0]*W[0][0];
+                       row2[0] = j2[0][0]*W[1][0];
+                       row3[0] = j2[0][0]*W[2][0];
+
+                       row1[3] = j2[1][0]*W[0][1] + j2[2][0]*W[0][2];
+                       row2[3] = j2[1][0]*W[1][1] + j2[2][0]*W[1][2];
+                       row3[3] = j2[1][0]*W[2][1] + j2[2][0]*W[2][2];
+               }
+
+               if (v2->flag & PVERT_INTERIOR) {
+                       vid[1] = v2->u.id;
+                       vid[4] = ninterior + v2->u.id;
+
+                       sys->J2dt[e1->u.id][1] = j2[0][1] = p_abf_compute_sin_product(sys, v2, e1->u.id)*wi1;
+                       sys->J2dt[e2->u.id][1] = j2[1][1] = 1.0*wi2;
+                       sys->J2dt[e3->u.id][1] = j2[2][1] = p_abf_compute_sin_product(sys, v2, e3->u.id)*wi3;
+
+                       nlRightHandSideAdd(0, v2->u.id, j2[1][1]*beta[1]);
+                       nlRightHandSideAdd(0, ninterior + v2->u.id, j2[0][1]*beta[0] + j2[2][1]*beta[2]);
+
+                       row1[1] = j2[1][1]*W[0][1];
+                       row2[1] = j2[1][1]*W[1][1];
+                       row3[1] = j2[1][1]*W[2][1];
+
+                       row1[4] = j2[0][1]*W[0][0] + j2[2][1]*W[0][2];
+                       row2[4] = j2[0][1]*W[1][0] + j2[2][1]*W[1][2];
+                       row3[4] = j2[0][1]*W[2][0] + j2[2][1]*W[2][2];
+               }
+
+               if (v3->flag & PVERT_INTERIOR) {
+                       vid[2] = v3->u.id;
+                       vid[5] = ninterior + v3->u.id;
+
+                       sys->J2dt[e1->u.id][2] = j2[0][2] = p_abf_compute_sin_product(sys, v3, e1->u.id)*wi1;
+                       sys->J2dt[e2->u.id][2] = j2[1][2] = p_abf_compute_sin_product(sys, v3, e2->u.id)*wi2;
+                       sys->J2dt[e3->u.id][2] = j2[2][2] = 1.0*wi3;
+
+                       nlRightHandSideAdd(0, v3->u.id, j2[2][2]*beta[2]);
+                       nlRightHandSideAdd(0, ninterior + v3->u.id, j2[0][2]*beta[0] + j2[1][2]*beta[1]);
+
+                       row1[2] = j2[2][2]*W[0][2];
+                       row2[2] = j2[2][2]*W[1][2];
+                       row3[2] = j2[2][2]*W[2][2];
+
+                       row1[5] = j2[0][2]*W[0][0] + j2[1][2]*W[0][1];
+                       row2[5] = j2[0][2]*W[1][0] + j2[1][2]*W[1][1];
+                       row3[5] = j2[0][2]*W[2][0] + j2[1][2]*W[2][1];
+               }
+
+               for (i = 0; i < 3; i++) {
+                       int r = vid[i];
+
+                       if (r == -1)
+                               continue;
+
+                       for (j = 0; j < 6; j++) {
+                               int c = vid[j];
+
+                               if (c == -1)
+                                       continue;
+
+                               if (i == 0)
+                                       nlMatrixAdd(r, c, j2[0][i]*row1[j]);
+                               else
+                                       nlMatrixAdd(r + ninterior, c, j2[0][i]*row1[j]);
+
+                               if (i == 1)
+                                       nlMatrixAdd(r, c, j2[1][i]*row2[j]);
+                               else
+                                       nlMatrixAdd(r + ninterior, c, j2[1][i]*row2[j]);
+
+
+                               if (i == 2)
+                                       nlMatrixAdd(r, c, j2[2][i]*row3[j]);
+                               else
+                                       nlMatrixAdd(r + ninterior, c, j2[2][i]*row3[j]);
+                       }
+               }
+       }
+
+       nlEnd(NL_MATRIX);
+
+       nlEnd(NL_SYSTEM);
+
+       success = nlSolve();
+
+       if (success) {
+               for (f=chart->faces; f; f=f->nextlink) {
+                       float dlambda1, pre[3], dalpha;
+                       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+                       PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+
+                       pre[0] = pre[1] = pre[2] = 0.0;
+
+                       if (v1->flag & PVERT_INTERIOR) {
+                               float x = nlGetVariable(0, v1->u.id);
+                               float x2 = nlGetVariable(0, ninterior + v1->u.id);
+                               pre[0] += sys->J2dt[e1->u.id][0]*x;
+                               pre[1] += sys->J2dt[e2->u.id][0]*x2;
+                               pre[2] += sys->J2dt[e3->u.id][0]*x2;
+                       }
+
+                       if (v2->flag & PVERT_INTERIOR) {
+                               float x = nlGetVariable(0, v2->u.id);
+                               float x2 = nlGetVariable(0, ninterior + v2->u.id);
+                               pre[0] += sys->J2dt[e1->u.id][1]*x2;
+                               pre[1] += sys->J2dt[e2->u.id][1]*x;
+                               pre[2] += sys->J2dt[e3->u.id][1]*x2;
+                       }
+
+                       if (v3->flag & PVERT_INTERIOR) {
+                               float x = nlGetVariable(0, v3->u.id);
+                               float x2 = nlGetVariable(0, ninterior + v3->u.id);
+                               pre[0] += sys->J2dt[e1->u.id][2]*x2;
+                               pre[1] += sys->J2dt[e2->u.id][2]*x2;
+                               pre[2] += sys->J2dt[e3->u.id][2]*x;
+                       }
+
+                       dlambda1 = pre[0] + pre[1] + pre[2];
+                       dlambda1 = sys->dstar[f->u.id]*(sys->bstar[f->u.id] - dlambda1);
+                       
+                       sys->lambdaTriangle[f->u.id] += dlambda1;
+
+                       dalpha = (sys->bAlpha[e1->u.id] - dlambda1);
+                       sys->alpha[e1->u.id] += dalpha/sys->weight[e1->u.id] - pre[0];
+
+                       dalpha = (sys->bAlpha[e2->u.id] - dlambda1);
+                       sys->alpha[e2->u.id] += dalpha/sys->weight[e2->u.id] - pre[1];
+
+                       dalpha = (sys->bAlpha[e3->u.id] - dlambda1);
+                       sys->alpha[e3->u.id] += dalpha/sys->weight[e3->u.id] - pre[2];
+
+                       /* clamp */
+                       e = f->edge;
+                       do {
+                               if (sys->alpha[e->u.id] > M_PI)
+                                       sys->alpha[e->u.id] = M_PI;
+                               else if (sys->alpha[e->u.id] < 0.0f)
+                                       sys->alpha[e->u.id] = 0.0f;
+                       } while (e != f->edge);
+               }
+
+               for (i = 0; i < ninterior; i++) {
+                       sys->lambdaPlanar[i] += nlGetVariable(0, i);
+                       sys->lambdaLength[i] += nlGetVariable(0, ninterior + i);
+               }
+       }
+
+       nlDeleteContext(nlGetCurrent());
+
+       return success;
+}
+
+static PBool p_chart_abf_solve(PChart *chart)
+{
+       PVert *v;
+       PFace *f;
+       PEdge *e, *e1, *e2, *e3;
+       PAbfSystem sys;
+       int i;
+       float lastnorm, limit = (chart->nfaces > 100)? 1.0f: 0.001f;
+
+       /* setup id's */
+       sys.ninterior = sys.nfaces = sys.nangles = 0;
+
+       for (v=chart->verts; v; v=v->nextlink) {
+               if (p_vert_interior(v)) {
+                       v->flag |= PVERT_INTERIOR;
+                       v->u.id = sys.ninterior++;
+               }
+               else
+                       v->flag &= ~PVERT_INTERIOR;
+       }
+
+       for (f=chart->faces; f; f=f->nextlink) {
+               e1 = f->edge; e2 = e1->next; e3 = e2->next;
+               f->u.id = sys.nfaces++;
+
+               /* angle id's are conveniently stored in half edges */
+               e1->u.id = sys.nangles++;
+               e2->u.id = sys.nangles++;
+               e3->u.id = sys.nangles++;
+       }
+
+       p_abf_setup_system(&sys);
+
+       /* compute initial angles */
+       for (f=chart->faces; f; f=f->nextlink) {
+               float a1, a2, a3;
+
+               e1 = f->edge; e2 = e1->next; e3 = e2->next;
+               p_face_angles(f, &a1, &a2, &a3);
+
+               if (a1 < sys.minangle)
+                       a1 = sys.minangle;
+               else if (a1 > sys.maxangle)
+                       a1 = sys.maxangle;
+               if (a2 < sys.minangle)
+                       a2 = sys.minangle;
+               else if (a2 > sys.maxangle)
+                       a2 = sys.maxangle;
+               if (a3 < sys.minangle)
+                       a3 = sys.minangle;
+               else if (a3 > sys.maxangle)
+                       a3 = sys.maxangle;
+
+               sys.alpha[e1->u.id] = sys.beta[e1->u.id] = a1;
+               sys.alpha[e2->u.id] = sys.beta[e2->u.id] = a2;
+               sys.alpha[e3->u.id] = sys.beta[e3->u.id] = a3;
+
+               sys.weight[e1->u.id] = 2.0/(a1*a1);
+               sys.weight[e2->u.id] = 2.0/(a2*a2);
+               sys.weight[e3->u.id] = 2.0/(a3*a3);
+       }
+
+       for (v=chart->verts; v; v=v->nextlink) {
+               if (v->flag & PVERT_INTERIOR) {
+                       float anglesum = 0.0, scale;
+
+                       e = v->edge;
+                       do {
+                               anglesum += sys.beta[e->u.id];
+                               e = e->next->next->pair;
+                       } while (e && (e != v->edge));
+
+                       scale = (anglesum == 0.0f)? 0.0f: 2*M_PI/anglesum;
+
+                       e = v->edge;
+                       do {
+                               sys.beta[e->u.id] = sys.alpha[e->u.id] = sys.beta[e->u.id]*scale;
+                               e = e->next->next->pair;
+                       } while (e && (e != v->edge));
+               }
+       }
+
+       if (sys.ninterior > 0) {
+               p_abf_compute_sines(&sys);
+
+               /* iteration */
+               lastnorm = 1e10;
+
+               for (i = 0; i < ABF_MAX_ITER; i++) {
+                       float norm = p_abf_compute_gradient(&sys, chart);
+
+                       lastnorm = norm;
+
+                       if (norm < limit)
+                               break;
+
+                       if (!p_abf_matrix_invert(&sys, chart)) {
+                               param_warning("ABF failed to invert matrix.");
+                               p_abf_free_system(&sys);
+                               return P_FALSE;
+                       }
+
+                       p_abf_compute_sines(&sys);
+               }
+
+               if (i == ABF_MAX_ITER) {
+                       param_warning("ABF maximum iterations reached.");
+                       p_abf_free_system(&sys);
+                       return P_FALSE;
+               }
+       }
+
+       chart->u.lscm.abf_alpha = MEM_dupallocN(sys.alpha);
+       p_abf_free_system(&sys);
+
+       return P_TRUE;
+}
+
+/* Least Squares Conformal Maps */
+
+static void p_chart_pin_positions(PChart *chart, PVert **pin1, PVert **pin2)
+{
+       if (pin1 == pin2) {
+               /* degenerate case */
+               PFace *f = chart->faces;
+               *pin1 = f->edge->vert;
+               *pin2 = f->edge->next->vert;
+
+               (*pin1)->uv[0] = 0.0f;
+               (*pin1)->uv[1] = 0.5f;
+               (*pin2)->uv[0] = 1.0f;
+               (*pin2)->uv[1] = 0.5f;
+       }
+       else {
+               int diru, dirv, dirx, diry;
+               float sub[3];
+
+               VecSubf(sub, (*pin1)->co, (*pin2)->co);
+               sub[0] = fabs(sub[0]);
+               sub[1] = fabs(sub[1]);
+               sub[2] = fabs(sub[2]);
+
+               if ((sub[0] > sub[1]) && (sub[0] > sub[2])) {
+                       dirx = 0;
+                       diry = (sub[1] > sub[2])? 1: 2;
+               }
+               else if ((sub[1] > sub[0]) && (sub[1] > sub[2])) {
+                       dirx = 1;
+                       diry = (sub[0] > sub[2])? 0: 2;
+               }
+               else {
+                       dirx = 2;
+                       diry = (sub[0] > sub[1])? 0: 1;
+               }
+
+               if (dirx == 2) {
+                       diru = 1;
+                       dirv = 0;
+               }
+               else {
+                       diru = 0;
+                       dirv = 1;
+               }
+
+               (*pin1)->uv[diru] = (*pin1)->co[dirx];
+               (*pin1)->uv[dirv] = (*pin1)->co[diry];
+               (*pin2)->uv[diru] = (*pin2)->co[dirx];
+               (*pin2)->uv[dirv] = (*pin2)->co[diry];
+       }
+}
+
+static PBool p_chart_symmetry_pins(PChart *chart, PEdge *outer, PVert **pin1, PVert **pin2)
+{
+       PEdge *be, *lastbe = NULL, *maxe1 = NULL, *maxe2 = NULL, *be1, *be2;
+       PEdge *cure = NULL, *firste1 = NULL, *firste2 = NULL, *nextbe;
+       float maxlen = 0.0f, curlen = 0.0f, totlen = 0.0f, firstlen = 0.0f;
+       float len1, len2;
+       /* find longest series of verts split in the chart itself, these are
+          marked during construction */
+       be = outer;
+       lastbe = p_boundary_edge_prev(be);
+       do {
+               float len = p_edge_length(be);
+               totlen += len;
+
+               nextbe = p_boundary_edge_next(be);
+
+               if ((be->vert->flag & PVERT_SPLIT) ||
+                   (lastbe->vert->flag & nextbe->vert->flag & PVERT_SPLIT)) {
+                       if (!cure) {
+                               if (be == outer)
+                                       firste1 = be;
+                               cure = be;
+                       }
+                       else
+                               curlen += p_edge_length(lastbe);
+               }
+               else if (cure) {
+                       if (curlen > maxlen) {
+                               maxlen = curlen;
+                               maxe1 = cure;
+                               maxe2 = lastbe;
+                       }
+
+                       if (firste1 == cure) {
+                               firstlen = curlen;
+                               firste2 = lastbe;
+                       }
+
+                       curlen = 0.0f;
+                       cure = NULL;
+               }
+
+               lastbe = be;
+               be = nextbe;
+       } while(be != outer);
+
+       /* make sure we also count a series of splits over the starting point */
+       if (cure && (cure != outer)) {
+               firstlen += curlen + p_edge_length(be);
+
+               if (firstlen > maxlen) {
+                       maxlen = firstlen;
+                       maxe1 = cure;
+                       maxe2 = firste2;
+               }
+       }
+
+       if (!maxe1 || !maxe2 || (maxlen < 0.5f*totlen))
+               return P_FALSE;
+       
+       /* find pin1 in the split vertices */
+       be1 = maxe1;
+       be2 = maxe2;
+       len1 = 0.0f;
+       len2 = 0.0f;
+
+       do {
+               if (len1 < len2) {
+                       len1 += p_edge_length(be1);
+                       be1 = p_boundary_edge_next(be1);
+               }
+               else {
+                       be2 = p_boundary_edge_prev(be2);
+                       len2 += p_edge_length(be2);
+               }
+       } while (be1 != be2);
+
+       *pin1 = be1->vert;
+
+       /* find pin2 outside the split vertices */
+       be1 = maxe1;
+       be2 = maxe2;
+       len1 = 0.0f;
+       len2 = 0.0f;
+
+       do {
+               if (len1 < len2) {
+                       be1 = p_boundary_edge_prev(be1);
+                       len1 += p_edge_length(be1);
+               }
+               else {
+                       len2 += p_edge_length(be2);
+                       be2 = p_boundary_edge_next(be2);
+               }
+       } while (be1 != be2);
+
+       *pin2 = be1->vert;
+
+       p_chart_pin_positions(chart, pin1, pin2);
+
+       return P_TRUE;
+}
+
+static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **pin2)
+{
+       float minv[3], maxv[3], dirlen;
+       PVert *v, *minvert[3], *maxvert[3];
+       int i, dir;
+
+       /* find minimum and maximum verts over x/y/z axes */
+       minv[0] = minv[1] = minv[2] = 1e20;
+       maxv[0] = maxv[1] = maxv[2] = -1e20;
+
+       minvert[0] = minvert[1] = minvert[2] = NULL;
+       maxvert[0] = maxvert[1] = maxvert[2] = NULL;
+
+       for (v = chart->verts; v; v=v->nextlink) {
+               for (i = 0; i < 3; i++) {
+                       if (v->co[i] < minv[i]) {
+                               minv[i] = v->co[i];
+                               minvert[i] = v;
+                       }
+                       if (v->co[i] > maxv[i]) {
+                               maxv[i] = v->co[i];
+                               maxvert[i] = v;
+                       }
+               }
+       }
+
+       /* find axes with longest distance */
+       dir = 0;
+       dirlen = -1.0;
+
+       for (i = 0; i < 3; i++) {
+               if (maxv[i] - minv[i] > dirlen) {
+                       dir = i;
+                       dirlen = maxv[i] - minv[i];
+               }
+       }
+
+       *pin1 = minvert[dir];
+       *pin2 = maxvert[dir];
+
+       p_chart_pin_positions(chart, pin1, pin2);
+}
+
+static void p_chart_lscm_load_solution(PChart *chart)
+{
+       PVert *v;
+
+       for (v=chart->verts; v; v=v->nextlink) {
+               v->uv[0] = nlGetVariable(0, 2*v->u.id);
+               v->uv[1] = nlGetVariable(0, 2*v->u.id + 1);
+       }
+}
+
+static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
+{
+       PVert *v, *pin1, *pin2;
+       PBool select = P_FALSE, deselect = P_FALSE;
+       int npins = 0, id = 0;
+
+       /* give vertices matrix indices and count pins */
+       for (v=chart->verts; v; v=v->nextlink) {
+               if (v->flag & PVERT_PIN) {
+                       npins++;
+                       if (v->flag & PVERT_SELECT)
+                               select = P_TRUE;
+               }
+
+               if (!(v->flag & PVERT_SELECT))
+                       deselect = P_TRUE;
+       }
+
+       if ((live && (!select || !deselect)) || (npins == 1)) {
+               chart->u.lscm.context = NULL;
+       }
+       else {
+#if 0
+               p_chart_simplify_compute(chart);
+               p_chart_topological_sanity_check(chart);
+#endif
+
+               if (abf) {
+                       if (!p_chart_abf_solve(chart))
+                               param_warning("ABF solving failed: falling back to LSCM.\n");
+               }
+
+               if (npins <= 1) {
+                       /* not enough pins, lets find some ourself */
+                       PEdge *outer;
+
+                       p_chart_boundaries(chart, NULL, &outer);
+
+                       if (!p_chart_symmetry_pins(chart, outer, &pin1, &pin2))
+                               p_chart_extrema_verts(chart, &pin1, &pin2);
+
+                       chart->u.lscm.pin1 = pin1;
+                       chart->u.lscm.pin2 = pin2;
+               }
+               else {
+                       chart->flag |= PCHART_NOPACK;
+               }
+
+               for (v=chart->verts; v; v=v->nextlink)
+                       v->u.id = id++;
+
+               nlNewContext();
+               nlSolverParameteri(NL_NB_VARIABLES, 2*chart->nverts);
+               nlSolverParameteri(NL_NB_ROWS, 2*chart->nfaces);
+               nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
+
+               chart->u.lscm.context = nlGetCurrent();
+       }
+}
+
+static PBool p_chart_lscm_solve(PChart *chart)
+{
+       PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
+       PFace *f;
+       float *alpha = chart->u.lscm.abf_alpha;
+       int row;
+
+       nlMakeCurrent(chart->u.lscm.context);
+
+       nlBegin(NL_SYSTEM);
+
+#if 0
+       /* TODO: make loading pins work for simplify/complexify. */
+#endif
+
+       for (v=chart->verts; v; v=v->nextlink)
+               if (v->flag & PVERT_PIN)
+                       p_vert_load_pin_select_uvs(v); /* reload for live */
+
+       if (chart->u.lscm.pin1) {
+               nlLockVariable(2*pin1->u.id);
+               nlLockVariable(2*pin1->u.id + 1);
+               nlLockVariable(2*pin2->u.id);
+               nlLockVariable(2*pin2->u.id + 1);
+       
+               nlSetVariable(0, 2*pin1->u.id, pin1->uv[0]);
+               nlSetVariable(0, 2*pin1->u.id + 1, pin1->uv[1]);
+               nlSetVariable(0, 2*pin2->u.id, pin2->uv[0]);
+               nlSetVariable(0, 2*pin2->u.id + 1, pin2->uv[1]);
+       }
+       else {
+               /* set and lock the pins */
+               for (v=chart->verts; v; v=v->nextlink) {
+                       if (v->flag & PVERT_PIN) {
+                               nlLockVariable(2*v->u.id);
+                               nlLockVariable(2*v->u.id + 1);
+
+                               nlSetVariable(0, 2*v->u.id, v->uv[0]);
+                               nlSetVariable(0, 2*v->u.id + 1, v->uv[1]);
+                       }
+               }
+       }
+
+       /* construct matrix */
+
+       nlBegin(NL_MATRIX);
+
+       row = 0;
+       for (f=chart->faces; f; f=f->nextlink) {
+               PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+               PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+               float a1, a2, a3, ratio, cosine, sine;
+               float sina1, sina2, sina3, sinmax;
+
+               if (alpha) {
+                       /* use abf angles if passed on */
+                       a1 = *(alpha++);
+                       a2 = *(alpha++);
+                       a3 = *(alpha++);
+               }
+               else
+                       p_face_angles(f, &a1, &a2, &a3);
+
+               sina1 = sin(a1);
+               sina2 = sin(a2);
+               sina3 = sin(a3);
 
                sinmax = MAX3(sina1, sina2, sina3);
 
-               /* shift vertices to find most stable order */
-               #define SHIFT3(type, a, b, c) \
-                       { type tmp; tmp = a; a = c; c = b; b = tmp; }
+               /* shift vertices to find most stable order */
+               #define SHIFT3(type, a, b, c) \
+                       { type tmp; tmp = a; a = c; c = b; b = tmp; }
+
+               if (sina3 != sinmax) {
+                       SHIFT3(PVert*, v1, v2, v3);
+                       SHIFT3(float, a1, a2, a3);
+                       SHIFT3(float, sina1, sina2, sina3);
+
+                       if (sina2 == sinmax) {
+                               SHIFT3(PVert*, v1, v2, v3);
+                               SHIFT3(float, a1, a2, a3);
+                               SHIFT3(float, sina1, sina2, sina3);
+                       }
+               }
+
+               /* angle based lscm formulation */
+               ratio = (sina3 == 0.0f)? 1.0f: sina2/sina3;
+               cosine = cos(a1)*ratio;
+               sine = sina1*ratio;
+
+#if 0
+               nlBegin(NL_ROW);
+               nlCoefficient(2*v1->u.id,   cosine - 1.0);
+               nlCoefficient(2*v1->u.id+1, -sine);
+               nlCoefficient(2*v2->u.id,   -cosine);
+               nlCoefficient(2*v2->u.id+1, sine);
+               nlCoefficient(2*v3->u.id,   1.0);
+               nlEnd(NL_ROW);
+
+               nlBegin(NL_ROW);
+               nlCoefficient(2*v1->u.id,   sine);
+               nlCoefficient(2*v1->u.id+1, cosine - 1.0);
+               nlCoefficient(2*v2->u.id,   -sine);
+               nlCoefficient(2*v2->u.id+1, -cosine);
+               nlCoefficient(2*v3->u.id+1, 1.0);
+               nlEnd(NL_ROW);
+#else
+               nlMatrixAdd(row, 2*v1->u.id,   cosine - 1.0);
+               nlMatrixAdd(row, 2*v1->u.id+1, -sine);
+               nlMatrixAdd(row, 2*v2->u.id,   -cosine);
+               nlMatrixAdd(row, 2*v2->u.id+1, sine);
+               nlMatrixAdd(row, 2*v3->u.id,   1.0);
+               row++;
+
+               nlMatrixAdd(row, 2*v1->u.id,   sine);
+               nlMatrixAdd(row, 2*v1->u.id+1, cosine - 1.0);
+               nlMatrixAdd(row, 2*v2->u.id,   -sine);
+               nlMatrixAdd(row, 2*v2->u.id+1, -cosine);
+               nlMatrixAdd(row, 2*v3->u.id+1, 1.0);
+               row++;
+#endif
+       }
+
+       nlEnd(NL_MATRIX);
+
+       nlEnd(NL_SYSTEM);
+
+       if (nlSolveAdvanced(NULL, NL_TRUE)) {
+               p_chart_lscm_load_solution(chart);
+               return P_TRUE;
+       }
+       else {
+               for (v=chart->verts; v; v=v->nextlink) {
+                       v->uv[0] = 0.0f;
+                       v->uv[1] = 0.0f;
+               }
+       }
+
+       return P_FALSE;
+}
+
+static void p_chart_lscm_end(PChart *chart)
+{
+       if (chart->u.lscm.context)
+               nlDeleteContext(chart->u.lscm.context);
+       
+       if (chart->u.lscm.abf_alpha) {
+               MEM_freeN(chart->u.lscm.abf_alpha);
+               chart->u.lscm.abf_alpha = NULL;
+       }
+
+       chart->u.lscm.context = NULL;
+       chart->u.lscm.pin1 = NULL;
+       chart->u.lscm.pin2 = NULL;
+}
+
+/* Stretch */
+
+#define P_STRETCH_ITER 20
+
+static void p_stretch_pin_boundary(PChart *chart)
+{
+       PVert *v;
+
+       for(v=chart->verts; v; v=v->nextlink)
+               if (v->edge->pair == NULL)
+                       v->flag |= PVERT_PIN;
+               else
+                       v->flag &= ~PVERT_PIN;
+}
+
+static float p_face_stretch(PFace *f)
+{
+       float T, w, tmp[3];
+       float Ps[3], Pt[3];
+       float a, c, area;
+       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+       PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+
+       area = p_face_uv_area_signed(f);
+
+       if (area <= 0.0f) /* flipped face -> infinite stretch */
+               return 1e10f;
+       
+       w= 1.0f/(2.0f*area);
+
+       /* compute derivatives */
+       VecCopyf(Ps, v1->co);
+       VecMulf(Ps, (v2->uv[1] - v3->uv[1]));
+
+       VecCopyf(tmp, v2->co);
+       VecMulf(tmp, (v3->uv[1] - v1->uv[1]));
+       VecAddf(Ps, Ps, tmp);
+
+       VecCopyf(tmp, v3->co);
+       VecMulf(tmp, (v1->uv[1] - v2->uv[1]));
+       VecAddf(Ps, Ps, tmp);
+
+       VecMulf(Ps, w);
+
+       VecCopyf(Pt, v1->co);
+       VecMulf(Pt, (v3->uv[0] - v2->uv[0]));
+
+       VecCopyf(tmp, v2->co);
+       VecMulf(tmp, (v1->uv[0] - v3->uv[0]));
+       VecAddf(Pt, Pt, tmp);
+
+       VecCopyf(tmp, v3->co);
+       VecMulf(tmp, (v2->uv[0] - v1->uv[0]));
+       VecAddf(Pt, Pt, tmp);
+
+       VecMulf(Pt, w);
+
+       /* Sander Tensor */
+       a= Inpf(Ps, Ps);
+       c= Inpf(Pt, Pt);
+
+       T =  sqrt(0.5f*(a + c));
+       if (f->flag & PFACE_FILLED)
+               T *= 0.2;
+
+       return T;
+}
+
+static float p_stretch_compute_vertex(PVert *v)
+{
+       PEdge *e = v->edge;
+       float sum = 0.0f;
+
+       do {
+               sum += p_face_stretch(e->face);
+               e = p_wheel_edge_next(e);
+       } while (e && e != (v->edge));
+
+       return sum;
+}
+
+static void p_chart_stretch_minimize(PChart *chart, RNG *rng)
+{
+       PVert *v;
+       PEdge *e;
+       int j, nedges;
+       float orig_stretch, low, stretch_low, high, stretch_high, mid, stretch;
+       float orig_uv[2], dir[2], random_angle, trusted_radius;
+
+       for(v=chart->verts; v; v=v->nextlink) {
+               if((v->flag & PVERT_PIN) || !(v->flag & PVERT_SELECT))
+                       continue;
+
+               orig_stretch = p_stretch_compute_vertex(v);
+               orig_uv[0] = v->uv[0];
+               orig_uv[1] = v->uv[1];
+
+               /* move vertex in a random direction */
+               trusted_radius = 0.0f;
+               nedges = 0;
+               e = v->edge;
+
+               do {
+                       trusted_radius += p_edge_uv_length(e);
+                       nedges++;
+
+                       e = p_wheel_edge_next(e);
+               } while (e && e != (v->edge));
+
+               trusted_radius /= 2 * nedges;
+
+               random_angle = rng_getFloat(rng) * 2.0 * M_PI;
+               dir[0] = trusted_radius * cos(random_angle);
+               dir[1] = trusted_radius * sin(random_angle);
+
+               /* calculate old and new stretch */
+               low = 0;
+               stretch_low = orig_stretch;
+
+               Vec2Addf(v->uv, orig_uv, dir);
+               high = 1;
+               stretch = stretch_high = p_stretch_compute_vertex(v);
+
+               /* binary search for lowest stretch position */
+               for (j = 0; j < P_STRETCH_ITER; j++) {
+                       mid = 0.5 * (low + high);
+                       v->uv[0]= orig_uv[0] + mid*dir[0];
+                       v->uv[1]= orig_uv[1] + mid*dir[1];
+                       stretch = p_stretch_compute_vertex(v);
+
+                       if (stretch_low < stretch_high) {
+                               high = mid;
+                               stretch_high = stretch;
+                       }
+                       else {
+                               low = mid;
+                               stretch_low = stretch;
+                       }
+               }
+
+               /* no luck, stretch has increased, reset to old values */
+               if(stretch >= orig_stretch)
+                       Vec2Copyf(v->uv, orig_uv);
+       }
+}
+
+/* Minimum area enclosing rectangle for packing */
+
+static int p_compare_geometric_uv(const void *a, const void *b)
+{
+       PVert *v1 = *(PVert**)a;
+       PVert *v2 = *(PVert**)b;
+
+       if (v1->uv[0] < v2->uv[0])
+               return -1;
+       else if (v1->uv[0] == v2->uv[0]) {
+               if (v1->uv[1] < v2->uv[1])
+                       return -1;
+               else if (v1->uv[1] == v2->uv[1])
+                       return 0;
+               else
+                       return 1;
+       }
+       else
+               return 1;
+}
+
+static PBool p_chart_convex_hull(PChart *chart, PVert ***verts, int *nverts, int *right)
+{
+       /* Graham algorithm, taken from:
+        * http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/117225 */
+
+       PEdge *be, *e;
+       int npoints = 0, i, ulen, llen;
+       PVert **U, **L, **points, **p;
+
+       p_chart_boundaries(chart, NULL, &be);
+
+       if (!be)
+               return P_FALSE;
+
+       e = be;
+       do {
+               npoints++;
+               e = p_boundary_edge_next(e);
+       } while(e != be);
+
+       p = points = (PVert**)MEM_mallocN(sizeof(PVert*)*npoints*2, "PCHullpoints");
+       U = (PVert**)MEM_mallocN(sizeof(PVert*)*npoints, "PCHullU");
+       L = (PVert**)MEM_mallocN(sizeof(PVert*)*npoints, "PCHullL");
+
+       e = be;
+       do {
+               *p = e->vert;
+               p++;
+               e = p_boundary_edge_next(e);
+       } while(e != be);
+
+       qsort(points, npoints, sizeof(PVert*), p_compare_geometric_uv);
+
+       ulen = llen = 0;
+       for (p=points, i = 0; i < npoints; i++, p++) {
+               while ((ulen > 1) && (p_area_signed(U[ulen-2]->uv, (*p)->uv, U[ulen-1]->uv) <= 0))
+                       ulen--;
+               while ((llen > 1) && (p_area_signed(L[llen-2]->uv, (*p)->uv, L[llen-1]->uv) >= 0))
+                       llen--;
+
+               U[ulen] = *p;
+               ulen++;
+               L[llen] = *p;
+               llen++;
+       }
+
+       npoints = 0;
+       for (p=points, i = 0; i < ulen; i++, p++, npoints++)
+               *p = U[i];
+
+       /* the first and last point in L are left out, since they are also in U */
+       for (i = llen-2; i > 0; i--, p++, npoints++)
+               *p = L[i];
+
+       *verts = points;
+       *nverts = npoints;
+       *right = ulen - 1;
+
+       MEM_freeN(U);
+       MEM_freeN(L);
+
+       return P_TRUE;
+}
+
+float p_rectangle_area(float *p1, float *dir, float *p2, float *p3, float *p4)
+{
+       /* given 4 points on the rectangle edges and the direction of on edge,
+          compute the area of the rectangle */
+
+       float orthodir[2], corner1[2], corner2[2], corner3[2];
+
+       orthodir[0] = dir[1];
+       orthodir[1] = -dir[0];
+
+       if (!p_intersect_line_2d_dir(p1, dir, p2, orthodir, corner1))
+               return 1e10;
+
+       if (!p_intersect_line_2d_dir(p1, dir, p4, orthodir, corner2))
+               return 1e10;
+
+       if (!p_intersect_line_2d_dir(p3, dir, p4, orthodir, corner3))
+               return 1e10;
+
+       return Vec2Lenf(corner1, corner2)*Vec2Lenf(corner2, corner3);
+}
+
+static float p_chart_minimum_area_angle(PChart *chart)
+{
+       /* minimum area enclosing rectangle with rotating calipers, info:
+        * http://cgm.cs.mcgill.ca/~orm/maer.html */
+
+       float rotated, minarea, minangle, area, len;
+       float *angles, miny, maxy, v[2], a[4], mina;
+       int npoints, right, mini, maxi, i, idx[4], nextidx;
+       PVert **points, *p1, *p2, *p3, *p4, *p1n;
+
+       /* compute convex hull */
+       if (!p_chart_convex_hull(chart, &points, &npoints, &right))
+               return 0.0;
+
+       /* find left/top/right/bottom points, and compute angle for each point */
+       angles = MEM_mallocN(sizeof(float)*npoints, "PMinAreaAngles");
+
+       mini = maxi = 0;
+       miny = 1e10;
+       maxy = -1e10;
+
+       for (i = 0; i < npoints; i++) {
+               p1 = (i == 0)? points[npoints-1]: points[i-1];
+               p2 = points[i];
+               p3 = (i == npoints-1)? points[0]: points[i+1];
+
+               angles[i] = M_PI - p_vec2_angle(p1->uv, p2->uv, p3->uv);
+
+               if (points[i]->uv[1] < miny) {
+                       miny = points[i]->uv[1];
+                       mini = i;
+               }
+               if (points[i]->uv[1] > maxy) {
+                       maxy = points[i]->uv[1];
+                       maxi = i;
+               }
+       }
+
+       /* left, top, right, bottom */
+       idx[0] = 0;
+       idx[1] = maxi;
+       idx[2] = right;
+       idx[3] = mini;
+
+       v[0] = points[idx[0]]->uv[0];
+       v[1] = points[idx[0]]->uv[1] + 1.0f;
+       a[0] = p_vec2_angle(points[(idx[0]+1)%npoints]->uv, points[idx[0]]->uv, v);
+
+       v[0] = points[idx[1]]->uv[0] + 1.0f;
+       v[1] = points[idx[1]]->uv[1];
+       a[1] = p_vec2_angle(points[(idx[1]+1)%npoints]->uv, points[idx[1]]->uv, v);
+
+       v[0] = points[idx[2]]->uv[0];
+       v[1] = points[idx[2]]->uv[1] - 1.0f;
+       a[2] = p_vec2_angle(points[(idx[2]+1)%npoints]->uv, points[idx[2]]->uv, v);
+
+       v[0] = points[idx[3]]->uv[0] - 1.0f;
+       v[1] = points[idx[3]]->uv[1];
+       a[3] = p_vec2_angle(points[(idx[3]+1)%npoints]->uv, points[idx[3]]->uv, v);
+
+       /* 4 rotating calipers */
+
+       rotated = 0.0;
+       minarea = 1e10;
+       minangle = 0.0;
+
+       while (rotated <= M_PI/2) { /* INVESTIGATE: how far to rotate? */
+               /* rotate with the smallest angle */
+               mini = 0;
+               mina = 1e10;
+
+               for (i = 0; i < 4; i++)
+                       if (a[i] < mina) {
+                               mina = a[i];
+                               mini = i;
+                       }
+
+               rotated += mina;
+               nextidx = (idx[mini]+1)%npoints;
+
+               a[mini] = angles[nextidx];
+               a[(mini+1)%4] = a[(mini+1)%4] - mina;
+               a[(mini+2)%4] = a[(mini+2)%4] - mina;
+               a[(mini+3)%4] = a[(mini+3)%4] - mina;
+
+               /* compute area */
+               p1 = points[idx[mini]];
+               p1n = points[nextidx];
+               p2 = points[idx[(mini+1)%4]];
+               p3 = points[idx[(mini+2)%4]];
+               p4 = points[idx[(mini+3)%4]];
+
+               len = Vec2Lenf(p1->uv, p1n->uv);
+
+               if (len > 0.0f) {
+                       len = 1.0/len;
+                       v[0] = (p1n->uv[0] - p1->uv[0])*len;
+                       v[1] = (p1n->uv[1] - p1->uv[1])*len;
+
+                       area = p_rectangle_area(p1->uv, v, p2->uv, p3->uv, p4->uv);
+
+                       /* remember smallest area */
+                       if (area < minarea) {
+                               minarea = area;
+                               minangle = rotated;
+                       }
+               }
+
+               idx[mini] = nextidx;
+       }
+
+       /* try keeping rotation as small as possible */
+       if (minangle > M_PI/4)
+               minangle -= M_PI/2;
+
+       MEM_freeN(angles);
+       MEM_freeN(points);
+
+       return minangle;
+}
+
+void p_chart_rotate_minimum_area(PChart *chart)
+{
+       float angle = p_chart_minimum_area_angle(chart);
+       float sine = sin(angle);
+       float cosine = cos(angle);
+       PVert *v;
+
+       for (v = chart->verts; v; v=v->nextlink) {
+               float oldu = v->uv[0], oldv = v->uv[1];
+               v->uv[0] = cosine*oldu - sine*oldv;
+               v->uv[1] = sine*oldu + cosine*oldv;
+       }
+}
+
+/* Area Smoothing */
+
+/* 2d bsp tree for inverse mapping - that's a bit silly */
+
+typedef struct SmoothTriangle {
+       float co1[2], co2[2], co3[2];
+       float oco1[2], oco2[2], oco3[2];
+} SmoothTriangle;
+
+typedef struct SmoothNode {
+       struct SmoothNode *c1, *c2;
+       SmoothTriangle **tri;
+       float split;
+       int axis, ntri;
+} SmoothNode;
+
+static void p_barycentric_2d(float *v1, float *v2, float *v3, float *p, float *b)
+{
+       float a[2], c[2], h[2], div;
+
+       a[0] = v2[0] - v1[0];
+       a[1] = v2[1] - v1[1];
+       c[0] = v3[0] - v1[0];
+       c[1] = v3[1] - v1[1];
+
+       div = a[0]*c[1] - a[1]*c[0];
+
+       if (div == 0.0f) {
+               b[0] = 1.0f/3.0f;
+               b[1] = 1.0f/3.0f;
+               b[2] = 1.0f/3.0f;
+       }
+       else {
+               h[0] = p[0] - v1[0];
+               h[1] = p[1] - v1[1];
+
+               div = 1.0f/div;
 
-               if (sina3 != sinmax) {
-                       SHIFT3(PVert*, v1, v2, v3);
-                       SHIFT3(float, a1, a2, a3);
-                       SHIFT3(float, sina1, sina2, sina3);
+               b[1] = (h[0]*c[1] - h[1]*c[0])*div;
+               b[2] = (a[0]*h[1] - a[1]*h[0])*div;
+               b[0] = 1.0 - b[1] - b[2];
+       }
+}
 
-                       if (sina2 == sinmax) {
-                               SHIFT3(PVert*, v1, v2, v3);
-                               SHIFT3(float, a1, a2, a3);
-                               SHIFT3(float, sina1, sina2, sina3);
+static PBool p_triangle_inside(SmoothTriangle *t, float *co)
+{
+       float b[3];
+
+       p_barycentric_2d(t->co1, t->co2, t->co3, co, b);
+
+       if ((b[0] >= 0.0) && (b[1] >= 0.0) && (b[2] >= 0.0f)) {
+               co[0] = t->oco1[0]*b[0] + t->oco2[0]*b[1] + t->oco3[0]*b[2];
+               co[1] = t->oco1[1]*b[0] + t->oco2[1]*b[1] + t->oco3[1]*b[2];
+               return P_TRUE;
+       }
+
+       return P_FALSE;
+}
+
+static SmoothNode *p_node_new(MemArena *arena, SmoothTriangle **tri, int ntri, float *bmin, float *bmax, int depth)
+{
+       SmoothNode *node = BLI_memarena_alloc(arena, sizeof *node);
+       int axis, i, t1size = 0, t2size = 0;
+       float split, mi, mx;
+       SmoothTriangle **t1, **t2, *t;
+
+       node->tri = tri;
+       node->ntri = ntri;
+
+       if (ntri <= 10 || depth >= 15)
+               return node;
+       
+       t1 = MEM_mallocN(sizeof(SmoothTriangle)*ntri, "PNodeTri1");
+       t2 = MEM_mallocN(sizeof(SmoothTriangle)*ntri, "PNodeTri1");
+
+       axis = (bmax[0] - bmin[0] > bmax[1] - bmin[1])? 0: 1;
+       split = 0.5f*(bmin[axis] + bmax[axis]);
+
+       for (i = 0; i < ntri; i++) {
+               t = tri[i];
+
+               if ((t->co1[axis] <= split) || (t->co2[axis] <= split) || (t->co3[axis] <= split)) {
+                       t1[t1size] = t;
+                       t1size++;
+               }
+               if ((t->co1[axis] >= split) || (t->co2[axis] >= split) || (t->co3[axis] >= split)) {
+                       t2[t2size] = t;
+                       t2size++;
+               }
+       }
+
+       if ((t1size == t2size) && (t1size == ntri)) {
+               MEM_freeN(t1);
+               MEM_freeN(t2);
+               return node;
+       }
+       
+       node->tri = NULL;
+       node->ntri = 0;
+       MEM_freeN(tri);
+
+       node->axis = axis;
+       node->split = split;
+
+       mi = bmin[axis];
+       mx = bmax[axis];
+       bmax[axis] = split;
+       node->c1 = p_node_new(arena, t1, t1size, bmin, bmax, depth+1);
+
+       bmin[axis] = bmax[axis];
+       bmax[axis] = mx;
+       node->c2 = p_node_new(arena, t2, t2size, bmin, bmax, depth+1);
+
+       return node;
+}
+
+static void p_node_delete(SmoothNode *node)
+{
+       if (node->c1)
+               p_node_delete(node->c1);
+       if (node->c2)
+               p_node_delete(node->c2);
+       if (node->tri)
+               MEM_freeN(node->tri);
+}
+
+static PBool p_node_intersect(SmoothNode *node, float *co)
+{
+       int i;
+
+       if (node->tri) {
+               for (i = 0; i < node->ntri; i++)
+                       if (p_triangle_inside(node->tri[i], co))
+                               return P_TRUE;
+
+               return P_FALSE;
+       }
+       else {
+               if (co[node->axis] < node->split)
+                       return p_node_intersect(node->c1, co);
+               else
+                       return p_node_intersect(node->c2, co);
+       }
+
+}
+
+/* smooothing */
+
+static int p_compare_float(const void *a, const void *b)
+{
+       if (*((float*)a) < *((float*)b))
+               return -1;
+       else if (*((float*)a) == *((float*)b))
+               return 0;
+       else
+               return 1;
+}
+
+static float p_smooth_median_edge_length(PChart *chart)
+{
+       PEdge *e;
+       float *lengths = MEM_mallocN(sizeof(chart->edges)*chart->nedges, "PMedianLength");
+       float median;
+       int i;
+
+       /* ok, so i'm lazy */
+       for (i=0, e=chart->edges; e; e=e->nextlink, i++)
+               lengths[i] = p_edge_length(e);
+       
+       qsort(lengths, i, sizeof(float), p_compare_float);
+
+       median = lengths[i/2];
+       MEM_freeN(lengths);
+
+       return median;
+}
+
+static float p_smooth_distortion(PEdge *e, float avg2d, float avg3d)
+{
+       float len2d = p_edge_uv_length(e)*avg3d;
+       float len3d = p_edge_length(e)*avg2d;
+
+       return (len3d == 0.0f)? 0.0f: len2d/len3d;
+}
+
+void p_smooth(PChart *chart)
+{
+       PEdge *e;
+       PVert *v;
+       PFace *f;
+       int j, it2, maxiter2, it;
+       int nedges = chart->nedges, nwheel, gridx, gridy;
+       int edgesx, edgesy, nsize, esize, i, x, y, maxiter, totiter;
+       float minv[2], maxv[2], median, invmedian, distortion, avglen2d, avglen3d;
+       float center[2], dx, dy, *nodes, dlimit, d, *oldnodesx, *oldnodesy;
+       float *nodesx, *nodesy, *hedges, *vedges, climit, moved, padding;
+       SmoothTriangle *triangles, *t, *t2, **tri, **trip;
+       SmoothNode *root;
+       MemArena *arena;
+
+       if (nedges == 0)
+               return;
+
+       p_chart_uv_bbox(chart, minv, maxv);
+       median = p_smooth_median_edge_length(chart)*0.10f;
+
+       if (median == 0.0)
+               return;
+
+       invmedian = 1.0/median;
+
+       /* compute edge distortion */
+       distortion = 0.0;
+       avglen2d = avglen3d = 0.0;
+
+       for (e=chart->edges; e; e=e->nextlink) {
+               avglen2d += p_edge_uv_length(e);
+               avglen3d += p_edge_length(e);
+       }
+
+       avglen2d /= nedges;
+       avglen3d /= nedges;
+
+       for (v=chart->verts; v; v=v->nextlink) {
+               v->u.distortion = 0.0;
+               nwheel = 0;
+
+               e = v->edge;
+               do {
+                       v->u.distortion += p_smooth_distortion(e, avglen2d, avglen3d);
+                       nwheel++;
+
+                       e = e->next->next->pair;
+               } while(e && (e != v->edge));
+
+               v->u.distortion /= nwheel;
+       }
+
+       /* need to do excessive grid size checking still */
+       center[0] = 0.5f*(minv[0] + maxv[0]);
+       center[1] = 0.5f*(minv[1] + maxv[1]);
+
+       dx = 0.5f*(maxv[0] - minv[0]);
+       dy = 0.5f*(maxv[1] - minv[1]);
+
+       padding = 0.15f;
+       dx += padding*dx + 2.0f*median;
+       dy += padding*dy + 2.0f*median;
+
+       gridx = (int)(dx*invmedian);
+       gridy = (int)(dy*invmedian);
+
+       minv[0] = center[0] - median*gridx;
+       minv[1] = center[1] - median*gridy;
+       maxv[0] = center[0] + median*gridx;
+       maxv[1] = center[1] + median*gridy;
+
+       /* create grid */
+       gridx = gridx*2 + 1;
+       gridy = gridy*2 + 1;
+
+       if ((gridx <= 2) || (gridy <= 2))
+               return;
+       
+       edgesx = gridx-1;
+       edgesy = gridy-1;
+       nsize = gridx*gridy;
+       esize = edgesx*edgesy;
+       
+       nodes = MEM_mallocN(sizeof(float)*nsize, "PSmoothNodes");
+       nodesx = MEM_mallocN(sizeof(float)*nsize, "PSmoothNodesX");
+       nodesy = MEM_mallocN(sizeof(float)*nsize, "PSmoothNodesY");
+       oldnodesx = MEM_mallocN(sizeof(float)*nsize, "PSmoothOldNodesX");
+       oldnodesy = MEM_mallocN(sizeof(float)*nsize, "PSmoothOldNodesY");
+       hedges = MEM_mallocN(sizeof(float)*esize, "PSmoothHEdges");
+       vedges = MEM_mallocN(sizeof(float)*esize, "PSmoothVEdges");
+
+       if (!nodes || !nodesx || !nodesy || !oldnodesx || !oldnodesy || !hedges || !vedges) {
+               if (nodes) MEM_freeN(nodes);
+               if (nodesx) MEM_freeN(nodesx);
+               if (nodesy) MEM_freeN(nodesy);
+               if (oldnodesx) MEM_freeN(oldnodesx);
+               if (oldnodesy) MEM_freeN(oldnodesy);
+               if (hedges) MEM_freeN(hedges);
+               if (vedges) MEM_freeN(vedges);
+
+               error("Not enough memory for area smoothing grid.");
+               return;
+       }
+
+       for (x = 0; x < gridx; x++) {
+               for (y = 0; y < gridy; y++) {
+                       i = x + y*gridx;
+
+                       nodesx[i] = minv[0] + median*x;
+                       nodesy[i] = minv[1] + median*y;
+
+                       nodes[i] = 1.0f;
+               }
+       }
+
+       /* embed in grid */
+       for (f=chart->faces; f; f=f->nextlink) {
+               PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+               float fmin[2], fmax[2];
+               int bx1, by1, bx2, by2;
+
+               INIT_MINMAX2(fmin, fmax);
+
+               DO_MINMAX2(e1->vert->uv, fmin, fmax);
+               DO_MINMAX2(e2->vert->uv, fmin, fmax);
+               DO_MINMAX2(e3->vert->uv, fmin, fmax);
+
+               bx1 = (int)((fmin[0] - minv[0])*invmedian);
+               by1 = (int)((fmin[1] - minv[1])*invmedian);
+               bx2 = (int)((fmax[0] - minv[0])*invmedian + 2);
+               by2 = (int)((fmax[1] - minv[1])*invmedian + 2);
+
+               for (x = bx1; x < bx2; x++) {
+                       for (y = by1; y < by2; y++) {
+                               float p[2], b[3];
+
+                               i = x + y*gridx;
+               
+                               p[0] = nodesx[i];
+                               p[1] = nodesy[i];
+
+                               p_barycentric_2d(e1->vert->uv, e2->vert->uv, e3->vert->uv, p, b);
+
+                               if ((b[0] > 0.0) && (b[1] > 0.0) && (b[2] > 0.0)) {
+                                       nodes[i] = e1->vert->u.distortion*b[0];
+                                       nodes[i] += e2->vert->u.distortion*b[1];
+                                       nodes[i] += e3->vert->u.distortion*b[2];
+                               }
                        }
                }
+       }
 
-               /* angle based lscm formulation */
-               ratio = sina2/sina3;
-               cosine = cos(a1)*ratio;
-               sine = sina1*ratio;
+       /* smooth the grid */
+       maxiter = 10;
+       totiter = 0;
+       climit = 0.00001f*nsize;
+
+       for (it = 0; it < maxiter; it++) {
+               moved = 0.0f;
+
+               for (x = 0; x < edgesx; x++) {
+                       for (y = 0; y < edgesy; y++) {
+                               i = x + y*gridx;
+                               j = x + y*edgesx;
+
+                               hedges[j] = (nodes[i] + nodes[i+1])*0.5f;
+                               vedges[j] = (nodes[i] + nodes[i+gridx])*0.5f;
+
+                               /* we do *inverse* mapping */
+                               hedges[j] = 1.0f/hedges[j];
+                               vedges[j] = 1.0f/vedges[j];
+                       }
+               }
+
+               maxiter2 = 50;
+               dlimit = 0.0001f;
+
+               for (it2 = 0; it2 < maxiter2; it2++) {
+                       d = 0.0f;
+                       totiter += 1;
+                       
+                       memcpy(oldnodesx, nodesx, sizeof(float)*nsize);
+                       memcpy(oldnodesy, nodesy, sizeof(float)*nsize);
+
+                       for (x=1; x < gridx-1; x++) {
+                               for (y=1; y < gridy-1; y++) {
+                                       float p[2], oldp[2], sum1, sum2, diff[2], length;
+
+                                       i = x + gridx*y;
+                                       j = x + edgesx*y;
+
+                                       oldp[0] = oldnodesx[i];
+                                       oldp[1] = oldnodesy[i];
+
+                                       sum1 = hedges[j-1]*oldnodesx[i-1];
+                                       sum1 += hedges[j]*oldnodesx[i+1];
+                                       sum1 += vedges[j-edgesx]*oldnodesx[i-gridx];
+                                       sum1 += vedges[j]*oldnodesx[i+gridx];
+
+                                       sum2 = hedges[j-1];
+                                       sum2 += hedges[j];
+                                       sum2 += vedges[j-edgesx];
+                                       sum2 += vedges[j];
+
+                                       nodesx[i] = sum1/sum2;
+
+                                       sum1 = hedges[j-1]*oldnodesy[i-1];
+                                       sum1 += hedges[j]*oldnodesy[i+1];
+                                       sum1 += vedges[j-edgesx]*oldnodesy[i-gridx];
+                                       sum1 += vedges[j]*oldnodesy[i+gridx];
+
+                                       nodesy[i] = sum1/sum2;
+                                       
+                                       p[0] = nodesx[i];
+                                       p[1] = nodesy[i];
+
+                                       diff[0] = p[0] - oldp[0];
+                                       diff[1] = p[1] - oldp[1];
+
+                                       length = sqrt(diff[0]*diff[0] + diff[1]*diff[1]);
+                                       d = MAX2(d, length);
+                                       moved += length;
+                               }
+                       }
+
+                       if (d < dlimit)
+                               break;
+               }
+
+               if (moved < climit)
+                       break;
+       }
+
+       MEM_freeN(oldnodesx);
+       MEM_freeN(oldnodesy);
+       MEM_freeN(hedges);
+       MEM_freeN(vedges);
+
+       /* create bsp */
+       t = triangles = MEM_mallocN(sizeof(SmoothTriangle)*esize*2, "PSmoothTris");
+       trip = tri = MEM_mallocN(sizeof(SmoothTriangle*)*esize*2, "PSmoothTriP");
+
+       if (!triangles || !tri) {
+               MEM_freeN(nodes);
+               MEM_freeN(nodesx);
+               MEM_freeN(nodesy);
+
+               if (triangles) MEM_freeN(triangles);
+               if (tri) MEM_freeN(tri);
+
+               error("Not enough memory for area smoothing grid.");
+               return;
+       }
+
+       for (x = 0; x < edgesx; x++) {
+               for (y = 0; y < edgesy; y++) {
+                       i = x + y*gridx;
+
+                       t->co1[0] = nodesx[i];
+                       t->co1[1] = nodesy[i];
+
+                       t->co2[0] = nodesx[i+1];
+                       t->co2[1] = nodesy[i+1];
+
+                       t->co3[0] = nodesx[i+gridx];
+                       t->co3[1] = nodesy[i+gridx];
+
+                       t->oco1[0] = minv[0] + x*median;
+                       t->oco1[1] = minv[1] + y*median;
+
+                       t->oco2[0] = minv[0] + (x+1)*median;
+                       t->oco2[1] = minv[1] + y*median;
+
+                       t->oco3[0] = minv[0] + x*median;
+                       t->oco3[1] = minv[1] + (y+1)*median;
+
+                       t2 = t+1;
+
+                       t2->co1[0] = nodesx[i+gridx+1];
+                       t2->co1[1] = nodesy[i+gridx+1];
+
+                       t2->oco1[0] = minv[0] + (x+1)*median;
+                       t2->oco1[1] = minv[1] + (y+1)*median;
+
+                       t2->co2[0] = t->co2[0]; t2->co2[1] = t->co2[1];
+                       t2->oco2[0] = t->oco2[0]; t2->oco2[1] = t->oco2[1];
+
+                       t2->co3[0] = t->co3[0]; t2->co3[1] = t->co3[1];
+                       t2->oco3[0] = t->oco3[0]; t2->oco3[1] = t->oco3[1];
+
+                       *trip = t; trip++; t++; 
+                       *trip = t; trip++; t++; 
+               }
+       }
+
+       MEM_freeN(nodes);
+       MEM_freeN(nodesx);
+       MEM_freeN(nodesy);
+
+       arena = BLI_memarena_new(1<<16);
+       root = p_node_new(arena, tri, esize*2, minv, maxv, 0);
+
+       for (v=chart->verts; v; v=v->nextlink)
+               if (!p_node_intersect(root, v->uv))
+                       param_warning("area smoothing error: couldn't find mapping triangle\n");
+
+       p_node_delete(root);
+       BLI_memarena_free(arena);
+       
+       MEM_freeN(triangles);
+}
+
+/* Exported */
+
+ParamHandle *param_construct_begin()
+{
+       PHandle *handle = MEM_callocN(sizeof*handle, "PHandle");
+       handle->construction_chart = p_chart_new(handle);
+       handle->state = PHANDLE_STATE_ALLOCATED;
+       handle->arena = BLI_memarena_new((1<<16));
+
+       handle->hash_verts = phash_new((PHashLink**)&handle->construction_chart->verts, 1);
+       handle->hash_edges = phash_new((PHashLink**)&handle->construction_chart->edges, 1);
+       handle->hash_faces = phash_new((PHashLink**)&handle->construction_chart->faces, 1);
+
+       return (ParamHandle*)handle;
+}
+
+void param_delete(ParamHandle *handle)
+{
+       PHandle *phandle = (PHandle*)handle;
+       int i;
+
+       param_assert((phandle->state == PHANDLE_STATE_ALLOCATED) ||
+                    (phandle->state == PHANDLE_STATE_CONSTRUCTED));
+
+       for (i = 0; i < phandle->ncharts; i++)
+               p_chart_delete(phandle->charts[i]);
+       
+       if (phandle->charts)
+               MEM_freeN(phandle->charts);
+
+       if (phandle->construction_chart) {
+               p_chart_delete(phandle->construction_chart);
+
+               phash_delete(phandle->hash_verts);
+               phash_delete(phandle->hash_edges);
+               phash_delete(phandle->hash_faces);
+       }
+
+       BLI_memarena_free(phandle->arena);
+       MEM_freeN(phandle);
+}
+
+void param_face_add(ParamHandle *handle, ParamKey key, int nverts,
+                    ParamKey *vkeys, float **co, float **uv,
+                    ParamBool *pin, ParamBool *select)
+{
+       PHandle *phandle = (PHandle*)handle;
+
+       param_assert(phash_lookup(phandle->hash_faces, key) == NULL);
+       param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
+       param_assert((nverts == 3) || (nverts == 4));
+
+       if (nverts == 4) {
+               if (p_quad_split_direction(phandle, co, vkeys)) {
+                       p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 2, pin, select);
+                       p_face_add_construct(phandle, key, vkeys, co, uv, 0, 2, 3, pin, select);
+               }
+               else {
+                       p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 3, pin, select);
+                       p_face_add_construct(phandle, key, vkeys, co, uv, 1, 2, 3, pin, select);
+               }
+       }
+       else
+               p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 2, pin, select);
+}
+
+void param_edge_set_seam(ParamHandle *handle, ParamKey *vkeys)
+{
+       PHandle *phandle = (PHandle*)handle;
+       PEdge *e;
+
+       param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
+
+       e = p_edge_lookup(phandle, vkeys);
+       if (e)
+               e->flag |= PEDGE_SEAM;
+}
+
+void param_construct_end(ParamHandle *handle, ParamBool fill, ParamBool impl)
+{
+       PHandle *phandle = (PHandle*)handle;
+       PChart *chart = phandle->construction_chart;
+       int i, j, nboundaries = 0;
+       PEdge *outer;
+
+       param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
+
+       phandle->ncharts = p_connect_pairs(phandle, impl);
+       phandle->charts = p_split_charts(phandle, chart, phandle->ncharts);
+
+       p_chart_delete(phandle->construction_chart);
+       phandle->construction_chart = NULL;
+
+       phash_delete(phandle->hash_verts);
+       phash_delete(phandle->hash_edges);
+       phash_delete(phandle->hash_faces);
+       phandle->hash_verts = phandle->hash_edges = phandle->hash_faces = NULL;
 
-               nlBegin(NL_ROW);
-               nlCoefficient(2*v1->u.index,   cosine - 1.0);
-               nlCoefficient(2*v1->u.index+1, -sine);
-               nlCoefficient(2*v2->u.index,   -cosine);
-               nlCoefficient(2*v2->u.index+1, sine);
-               nlCoefficient(2*v3->u.index,   1.0);
-               nlEnd(NL_ROW);
+       for (i = j = 0; i < phandle->ncharts; i++) {
+               PVert *v;
+               PChart *chart = phandle->charts[i];
 
-               nlBegin(NL_ROW);
-               nlCoefficient(2*v1->u.index,   sine);
-               nlCoefficient(2*v1->u.index+1, cosine - 1.0);
-               nlCoefficient(2*v2->u.index,   -sine);
-               nlCoefficient(2*v2->u.index+1, -cosine);
-               nlCoefficient(2*v3->u.index+1, 1.0);
-               nlEnd(NL_ROW);
-       }
+               p_chart_boundaries(chart, &nboundaries, &outer);
 
-       nlEnd(NL_MATRIX);
+               if (nboundaries == 0) {
+                       p_chart_delete(chart);
+                       continue;
+               }
 
-       nlEnd(NL_SYSTEM);
+               phandle->charts[j] = chart;
+               j++;
 
-       if (nlSolveAdvanced(NULL, NL_TRUE)) {
-               p_chart_lscm_load_solution(chart);
-               p_flush_uvs(chart);
+               if (fill && (nboundaries > 1))
+                       p_chart_fill_boundaries(chart, outer);
 
-               return P_TRUE;
+               for (v=chart->verts; v; v=v->nextlink)
+                       p_vert_load_pin_select_uvs(v);
        }
 
-       return P_FALSE;
-}
-
-static void p_chart_lscm_end(PChart *chart)
-{
-       if (chart->u.lscm.context)
-               nlDeleteContext(chart->u.lscm.context);
+       phandle->ncharts = j;
 
-       chart->u.lscm.context = NULL;
-       chart->u.lscm.singlepin = NULL;
-       chart->u.lscm.pin1 = NULL;
-       chart->u.lscm.pin2 = NULL;
+       phandle->state = PHANDLE_STATE_CONSTRUCTED;
 }
 
-void param_lscm_begin(ParamHandle *handle)
+void param_lscm_begin(ParamHandle *handle, ParamBool live, ParamBool abf)
 {
        PHandle *phandle = (PHandle*)handle;
+       PFace *f;
        int i;
 
        param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
        phandle->state = PHANDLE_STATE_LSCM;
 
-       for (i = 0; i < phandle->ncharts; i++)
-               p_chart_lscm_begin(phandle->charts[i]);
+       for (i = 0; i < phandle->ncharts; i++) {
+               for (f=phandle->charts[i]->faces; f; f=f->nextlink)
+                       p_face_backup_uvs(f);
+               p_chart_lscm_begin(phandle->charts[i], live, abf);
+       }
 }
 
 void param_lscm_solve(ParamHandle *handle)
@@ -1442,6 +4000,9 @@ void param_lscm_solve(ParamHandle *handle)
                if (chart->u.lscm.context) {
                        result = p_chart_lscm_solve(chart);
 
+                       if (result && !(chart->flag & PCHART_NOPACK))
+                               p_chart_rotate_minimum_area(chart);
+
                        if (!result || (chart->u.lscm.pin1))
                                p_chart_lscm_end(chart);
                }
@@ -1455,159 +4016,16 @@ void param_lscm_end(ParamHandle *handle)
 
        param_assert(phandle->state == PHANDLE_STATE_LSCM);
 
-       for (i = 0; i < phandle->ncharts; i++)
+       for (i = 0; i < phandle->ncharts; i++) {
                p_chart_lscm_end(phandle->charts[i]);
+#if 0
+               p_chart_complexify(phandle->charts[i]);
+#endif
+       }
 
        phandle->state = PHANDLE_STATE_CONSTRUCTED;
 }
 
-/* Stretch */
-
-#define P_STRETCH_ITER 20
-
-static void p_stretch_pin_boundary(PChart *chart)
-{
-       PVert *v;
-
-       for(v=(PVert*)chart->verts->first; v; v=v->link.next)
-               if (v->edge->pair == NULL)
-                       v->flag |= PVERT_PIN;
-               else
-                       v->flag &= ~PVERT_PIN;
-}
-
-static float p_face_stretch(PFace *f)
-{
-       float T, w, tmp[3];
-       float Ps[3], Pt[3];
-       float a, c, area;
-       PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-       PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-
-       area = p_face_uv_area_signed(f);
-
-       if (area <= 0.0f) /* flipped face -> infinite stretch */
-               return 1e10f;
-       
-       if (f->flag & PFACE_FILLED)
-               return 0.0f;
-
-       w= 1.0f/(2.0f*area);
-
-       /* compute derivatives */
-       VecCopyf(Ps, v1->co);
-       VecMulf(Ps, (v2->uv[1] - v3->uv[1]));
-
-       VecCopyf(tmp, v2->co);
-       VecMulf(tmp, (v3->uv[1] - v1->uv[1]));
-       VecAddf(Ps, Ps, tmp);
-
-       VecCopyf(tmp, v3->co);
-       VecMulf(tmp, (v1->uv[1] - v2->uv[1]));
-       VecAddf(Ps, Ps, tmp);
-
-       VecMulf(Ps, w);
-
-       VecCopyf(Pt, v1->co);
-       VecMulf(Pt, (v3->uv[0] - v2->uv[0]));
-
-       VecCopyf(tmp, v2->co);
-       VecMulf(tmp, (v1->uv[0] - v3->uv[0]));
-       VecAddf(Pt, Pt, tmp);
-
-       VecCopyf(tmp, v3->co);
-       VecMulf(tmp, (v2->uv[0] - v1->uv[0]));
-       VecAddf(Pt, Pt, tmp);
-
-       VecMulf(Pt, w);
-
-       /* Sander Tensor */
-       a= Inpf(Ps, Ps);
-       c= Inpf(Pt, Pt);
-
-       T = sqrt(0.5f*(a + c)*f->u.area3d);
-
-       return T;
-}
-
-static float p_stretch_compute_vertex(PVert *v)
-{
-       PEdge *e = v->edge;
-       float sum = 0.0f;
-
-       do {
-               sum += p_face_stretch(e->face);
-               e = p_wheel_edge_next(e);
-       } while (e && e != (v->edge));
-
-       return sum;
-}
-
-static void p_chart_stretch_minimize(PChart *chart, RNG *rng)
-{
-       PVert *v;
-       PEdge *e;
-       int j, nedges;
-       float orig_stretch, low, stretch_low, high, stretch_high, mid, stretch;
-       float orig_uv[2], dir[2], random_angle, trusted_radius;
-
-       for(v=(PVert*)chart->verts->first; v; v=v->link.next) {
-               if((v->flag & PVERT_PIN) || !(v->flag & PVERT_SELECT))
-                       continue;
-
-               orig_stretch = p_stretch_compute_vertex(v);
-               orig_uv[0] = v->uv[0];
-               orig_uv[1] = v->uv[1];
-
-               /* move vertex in a random direction */
-               trusted_radius = 0.0f;
-               nedges = 0;
-               e = v->edge;
-
-               do {
-                       trusted_radius += p_edge_uv_length(e);
-                       nedges++;
-
-                       e = p_wheel_edge_next(e);
-               } while (e && e != (v->edge));
-
-               trusted_radius /= 2 * nedges;
-
-               random_angle = rng_getFloat(rng) * 2.0 * M_PI;
-               dir[0] = trusted_radius * cos(random_angle);
-               dir[1] = trusted_radius * sin(random_angle);
-
-               /* calculate old and new stretch */
-               low = 0;
-               stretch_low = orig_stretch;
-
-               Vec2Addf(v->uv, orig_uv, dir);
-               high = 1;
-               stretch = stretch_high = p_stretch_compute_vertex(v);
-
-               /* binary search for lowest stretch position */
-               for (j = 0; j < P_STRETCH_ITER; j++) {
-                       mid = 0.5 * (low + high);
-                       v->uv[0]= orig_uv[0] + mid*dir[0];
-                       v->uv[1]= orig_uv[1] + mid*dir[1];
-                       stretch = p_stretch_compute_vertex(v);
-
-                       if (stretch_low < stretch_high) {
-                               high = mid;
-                               stretch_high = stretch;
-                       }
-                       else {
-                               low = mid;
-                               stretch_low = stretch;
-                       }
-               }
-
-               /* no luck, stretch has increased, reset to old values */
-               if(stretch >= orig_stretch)
-                       Vec2Copyf(v->uv, orig_uv);
-       }
-}
-
 void param_stretch_begin(ParamHandle *handle)
 {
        PHandle *phandle = (PHandle*)handle;
@@ -1620,22 +4038,31 @@ void param_stretch_begin(ParamHandle *handle)
        phandle->state = PHANDLE_STATE_STRETCH;
 
        phandle->rng = rng_new(31415926);
+       phandle->blend = 0.0f;
 
        for (i = 0; i < phandle->ncharts; i++) {
                chart = phandle->charts[i];
 
-               for (v=(PVert*)chart->verts->first; v; v=v->link.next)
-                       p_vert_load_select_uvs(v);
+               for (v=chart->verts; v; v=v->nextlink)
+                       v->flag &= ~PVERT_PIN; /* don't use user-defined pins */
 
                p_stretch_pin_boundary(chart);
 
-               for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
+               for (f=chart->faces; f; f=f->nextlink) {
                        p_face_backup_uvs(f);
                        f->u.area3d = p_face_area(f);
                }
        }
 }
 
+void param_stretch_blend(ParamHandle *handle, float blend)
+{
+       PHandle *phandle = (PHandle*)handle;
+
+       param_assert(phandle->state == PHANDLE_STATE_STRETCH);
+       phandle->blend = blend;
+}
+
 void param_stretch_iter(ParamHandle *handle)
 {
        PHandle *phandle = (PHandle*)handle;
@@ -1646,140 +4073,131 @@ void param_stretch_iter(ParamHandle *handle)
 
        for (i = 0; i < phandle->ncharts; i++) {
                chart = phandle->charts[i];
-
                p_chart_stretch_minimize(chart, phandle->rng);
-               p_flush_uvs(chart);
        }
 }
 
-void param_stretch_end(ParamHandle *handle, ParamBool restore)
+void param_stretch_end(ParamHandle *handle)
 {
        PHandle *phandle = (PHandle*)handle;
-       PChart *chart;
-       PFace *f;
-       int i;
 
        param_assert(phandle->state == PHANDLE_STATE_STRETCH);
        phandle->state = PHANDLE_STATE_CONSTRUCTED;
 
        rng_free(phandle->rng);
        phandle->rng = NULL;
-
-       if (restore) {
-               for (i = 0; i < phandle->ncharts; i++) {
-                       chart = phandle->charts[i];
-
-                       for (f=(PFace*)chart->faces->first; f; f=f->link.next)
-                               p_face_restore_uvs(f);
-               }
-       }
 }
 
-/* Packing */
-
-static PBool p_pack_try(PHandle *handle, float side)
+void param_smooth_area(ParamHandle *handle)
 {
-       PChart *chart;
-       float packx, packy, rowh, groupw, w, h;
+       PHandle *phandle = (PHandle*)handle;
        int i;
 
-       packx= packy= 0.0;
-       rowh= 0.0;
-       groupw= 1.0/sqrt(handle->ncharts);
-
-       for (i = 0; i < handle->ncharts; i++) {
-               chart = handle->charts[i];
-               w = chart->u.pack.size[0];
-               h = chart->u.pack.size[1];
-
-               if(w <= (1.0-packx)) {
-                       chart->u.pack.trans[0] = packx;
-                       chart->u.pack.trans[1] = packy;
+       param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
 
-                       packx += w;
-                       rowh= MAX2(rowh, h);
-               }
-               else {
-                       packy += rowh;
-                       packx = w;
-                       rowh = h;
+       for (i = 0; i < phandle->ncharts; i++) {
+               PChart *chart = phandle->charts[i];
+               PVert *v;
 
-                       chart->u.pack.trans[0] = 0.0;
-                       chart->u.pack.trans[1] = packy;
-               }
+               for (v=chart->verts; v; v=v->nextlink)
+                       v->flag &= ~PVERT_PIN;
 
-               if (rowh > side)
-                       return P_FALSE;
+               p_smooth(chart);
        }
-
-       return P_TRUE;
 }
-
-#define PACK_SEARCH_DEPTH 7
-
 void param_pack(ParamHandle *handle)
 {
-       PHandle *phandle = (PHandle*)handle;
+       /* box packing variables */
+       boxPack *boxarray, *box;
+       float tot_width, tot_height, scale;
+        
        PChart *chart;
-       float uv_area, area, trans[2], minside, maxside, totarea, side;
-       int i;
-
-       /* very simple rectangle packing */
-
+       int i, unpacked=0;
+       float trans[2];
+       
+       PHandle *phandle = (PHandle*)handle;
+       
+       
        if (phandle->ncharts == 0)
                return;
-
-       totarea = 0.0f;
-       maxside = 0.0f;
-
+       
+       /* we may not use all these boxes */
+       boxarray = MEM_mallocN( phandle->ncharts*sizeof(boxPack), "boxPack box");
+       
        for (i = 0; i < phandle->ncharts; i++) {
                chart = phandle->charts[i];
-
-               p_chart_area(chart, &uv_area, &area);
-               chart->u.pack.rescale = (uv_area > 0.0f)? area/uv_area: 0.0f;
-               totarea += uv_area*chart->u.pack.rescale;
-
+               
+               
+               if (chart->flag & PCHART_NOPACK) {
+                       unpacked++;
+                       continue;
+               }
+               
+               box = boxarray+(i-unpacked);
+               
                p_chart_uv_bbox(chart, trans, chart->u.pack.size);
+               
                trans[0] = -trans[0];
                trans[1] = -trans[1];
+               
                p_chart_uv_translate(chart, trans);
-               p_chart_uv_scale(chart, chart->u.pack.rescale);
-
-               chart->u.pack.size[0] -= trans[0];
-               chart->u.pack.size[1] -= trans[1];
-               chart->u.pack.size[0] *= chart->u.pack.rescale;
-               chart->u.pack.size[1] *= chart->u.pack.rescale;
-
-               maxside = MAX3(maxside, chart->u.pack.size[0], chart->u.pack.size[1]);
+               
+               box->w =  chart->u.pack.size[0] + trans[0];
+               box->h =  chart->u.pack.size[1] + trans[1];
+               box->index = i; /* warning this index skips PCHART_NOPACK boxes */
+       }
+       
+       boxPack2D(boxarray, phandle->ncharts-unpacked, &tot_width, &tot_height);
+       
+       if (tot_height>tot_width)
+               scale = 1.0/tot_height;
+       else
+               scale = 1.0/tot_width;
+       
+       for (i = 0; i < phandle->ncharts-unpacked; i++) {
+               box = boxarray+i;
+               trans[0] = box->x;
+               trans[1] = box->y;
+               
+               chart = phandle->charts[box->index];
+               p_chart_uv_translate(chart, trans);
+               p_chart_uv_scale(chart, scale);
        }
+       MEM_freeN(boxarray);
+}
 
-       return;
+void param_flush(ParamHandle *handle)
+{
+       PHandle *phandle = (PHandle*)handle;
+       PChart *chart;
+       int i;
 
-       printf("%f\n", maxside);
+       for (i = 0; i < phandle->ncharts; i++) {
+               chart = phandle->charts[i];
 
-       /* binary search over pack region size */
-       minside = sqrt(totarea);
-       maxside = (((int)sqrt(phandle->ncharts-1))+1)*maxside;
+               if ((phandle->state == PHANDLE_STATE_LSCM) && !chart->u.lscm.context)
+                       continue;
 
-       for (i = 0; i < PACK_SEARCH_DEPTH; i++) {
-               printf("%f %f\n", minside, maxside);
-               if (p_pack_try(phandle, (minside+maxside)*0.5f))
-                       minside = (minside+maxside)*0.5f;
+               if (phandle->blend == 0.0f)
+                       p_flush_uvs(chart);
                else
-                       maxside = (minside+maxside)*0.5f;
+                       p_flush_uvs_blend(chart, phandle->blend);
        }
+}
 
-       side = maxside + 1e-5; /* prevent floating point errors */
-       if (!p_pack_try(phandle, side))
-               printf("impossible\n");
+void param_flush_restore(ParamHandle *handle)
+{
+       PHandle *phandle = (PHandle*)handle;
+       PChart *chart;
+       PFace *f;
+       int i;
 
        for (i = 0; i < phandle->ncharts; i++) {
                chart = phandle->charts[i];
 
-               p_chart_uv_scale(chart, chart->u.pack.rescale/side);
-               trans[0] = chart->u.pack.trans[0]/side;
-               trans[1] = chart->u.pack.trans[1]/side;
-               p_chart_uv_translate(chart, trans);
+               for (f=chart->faces; f; f=f->nextlink)
+                       p_face_restore_uvs(f);
        }
 }