Mesh Deform Modifier
[blender.git] / source / blender / src / parametrizer.c
index c06c6e9dc4597b83ae1c3b41c2f64ae186382d92..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"
 
@@ -129,148 +131,6 @@ static PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
        return link;
 }
 
-/* Heap */
-
-#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)
-{
-       while (P_TRUE) {
-               int size = heap->size, smallest;
-               int l = PHEAP_LEFT(i);
-               int r = PHEAP_RIGHT(i);
-
-               smallest = ((l < size) && PHEAP_COMPARE(heap->tree[l], heap->tree[i]))? l: i;
-
-               if ((r < size) && PHEAP_COMPARE(heap->tree[r], heap->tree[smallest]))
-                       smallest = r;
-               
-               if (smallest == i)
-                       break;
-
-               PHEAP_SWAP(heap, i, smallest);
-               i = smallest;
-       }
-}
-
-static void pheap_up(PHeap *heap, int i)
-{
-       while (i > 0) {
-               int p = PHEAP_PARENT(i);
-
-               if (PHEAP_COMPARE(heap->tree[p], heap->tree[i]))
-                       break;
-
-               PHEAP_SWAP(heap, p, i);
-               i = p;
-       }
-}
-
-static PHeap *pheap_new()
-{
-       PHeap *heap = (PHeap*)MEM_callocN(sizeof(PHeap), "PHeap");
-       heap->bufsize = 1;
-       heap->tree = (PHeapLink**)MEM_mallocN(sizeof(PHeapLink*), "PHeapTree");
-       heap->arena = BLI_memarena_new(1<<16);
-
-       return heap;
-}
-
-static void pheap_delete(PHeap *heap)
-{
-       MEM_freeN(heap->tree);
-       BLI_memarena_free(heap->arena);
-       MEM_freeN(heap);
-}
-
-static PHeapLink *pheap_insert(PHeap *heap, float value, void *ptr)
-{
-       PHeapLink *link;
-
-       if ((heap->size + 1) > heap->bufsize) {
-               int newsize = heap->bufsize*2;
-
-               PHeapLink **ntree = (PHeapLink**)MEM_mallocN(newsize*sizeof(PHeapLink*), "PHeapTree");
-               memcpy(ntree, heap->tree, sizeof(PHeapLink*)*heap->size);
-               MEM_freeN(heap->tree);
-
-               heap->tree = ntree;
-               heap->bufsize = newsize;
-       }
-
-       param_assert(heap->size < heap->bufsize);
-
-       if (heap->freelinks) {
-               link = heap->freelinks;
-               heap->freelinks = (PHeapLink*)(((PHeapLink*)heap->freelinks)->ptr);
-       }
-       else
-               link = (PHeapLink*)BLI_memarena_alloc(heap->arena, sizeof *link);
-       link->value = value;
-       link->ptr = ptr;
-       link->index = heap->size;
-
-       heap->tree[link->index] = link;
-
-       heap->size++;
-
-       pheap_up(heap, heap->size-1);
-
-       return link;
-}
-
-#if 0
-static int pheap_empty(PHeap *heap)
-{
-       return (heap->size == 0);
-}
-
-static PHeapLink *pheap_toplink(PHeap *heap)
-{
-       return heap->tree[0];
-}
-#endif
-
-static void *pheap_popmin(PHeap *heap)
-{
-       void *ptr = heap->tree[0]->ptr;
-
-       heap->tree[0]->ptr = heap->freelinks;
-       heap->freelinks = heap->tree[0];
-
-       if (heap->size == 1)
-               heap->size--;
-       else {
-               PHEAP_SWAP(heap, 0, heap->size-1);
-               heap->size--;
-
-               pheap_down(heap, 0);
-       }
-
-       return ptr;
-}
-
-static void pheap_remove(PHeap *heap, PHeapLink *link)
-{
-       int i = link->index;
-
-       while (i > 0) {
-               int p = PHEAP_PARENT(i);
-
-               PHEAP_SWAP(heap, p, i);
-               i = p;
-       }
-
-       pheap_popmin(heap);
-}
-
 /* Geometry */
 
 static float p_vec_angle_cos(float *v1, float *v2, float *v3)
@@ -285,8 +145,8 @@ static float p_vec_angle_cos(float *v1, float *v2, float *v3)
        d2[1] = v3[1] - v2[1];
        d2[2] = v3[2] - v2[2];
 
-       Normalise(d1);
-       Normalise(d2);
+       Normalize(d1);
+       Normalize(d2);
 
        return d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2];
 }
@@ -352,23 +212,6 @@ static float p_face_uv_area_signed(PFace *f)
                    ((v3->uv[0] - v1->uv[0])*(v2->uv[1] - v1->uv[1])));
 }
 
-static float p_face_uv_area(PFace *f)
-{
-       return fabs(p_face_uv_area_signed(f));
-}
-
-static void p_chart_area(PChart *chart, float *uv_area, float *area)
-{
-       PFace *f;
-
-       *uv_area = *area = 0.0f;
-
-       for (f=chart->faces; f; f=f->nextlink) {
-               *uv_area += p_face_uv_area(f);
-               *area += p_face_area(f);
-       }
-}
-
 static float p_edge_length(PEdge *e)
 {
     PVert *v1 = e->vert, *v2 = e->next->vert;
@@ -686,6 +529,27 @@ static PEdge *p_edge_lookup(PHandle *handle, PHashKey *vkeys)
        return NULL;
 }
 
+static PBool p_face_exists(PHandle *handle, PHashKey *vkeys, int i1, int i2, int i3)
+{
+       PHashKey key = PHASH_edge(vkeys[i1], vkeys[i2]);
+       PEdge *e = (PEdge*)phash_lookup(handle->hash_edges, key);
+
+       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;
+               }
+
+               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");
@@ -720,12 +584,12 @@ static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
                uvp2 = ep->orig_uv;
        }
 
-       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;
@@ -879,7 +743,9 @@ static void p_split_vert(PChart *chart, PEdge *e)
 
        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;
@@ -1033,19 +899,25 @@ static PFace *p_face_add_fill(PChart *chart, PVert *v1, PVert *v2, PVert *v3)
        return f;
 }
 
-static PBool p_quad_split_direction(float **co)
+static PBool p_quad_split_direction(PHandle *handle, float **co, PHashKey *vkeys)
 {
-    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]);
+       float fac= VecLenf(co[0], co[2]) - VecLenf(co[1], co[3]);
+       PBool dir = (fac <= 0.0f);
 
-       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]);
+       /* 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;
+       }
+       else {
+               if (p_face_exists(handle,vkeys,0,1,3) || p_face_exists(handle,vkeys,1,2,3))
+                       return !dir;
+       }
 
-       return (a1 > a2);
+       return dir;
 }
 
 /* Construction: boundary filling */
@@ -1116,13 +988,13 @@ 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 = p_boundary_edge_next(e);
        } while(e != be);
@@ -1133,20 +1005,20 @@ 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 = (PEdge*)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;
@@ -1181,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)
@@ -1216,9 +1088,9 @@ static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
     }
 }
 
+#if 0
 /* Polygon kernel for inserting uv's non overlapping */
 
-#if 0
 static int p_polygon_point_in(float *cp1, float *cp2, float *p)
 {
        if ((cp1[0] == p[0]) && (cp1[1] == p[1]))
@@ -1332,7 +1204,9 @@ static void p_polygon_kernel_center(float (*points)[2], int npoints, float *cent
        MEM_freeN(oldpoints);
        MEM_freeN(newpoints);
 }
+#endif
 
+#if 0
 /* Edge Collapser */
 
 int NCOLLAPSE = 1;
@@ -2012,7 +1886,7 @@ static void p_chart_simplify_compute(PChart *chart)
           collapsed may then be view as stacks, where the next collapse/split
           is at the top of the respective lists. */
 
-       PHeap *heap = pheap_new();
+       Heap *heap = BLI_heap_new();
        PVert *v, **wheelverts;
        PEdge *collapsededges = NULL, *e;
        int nwheelverts, i, ncollapsed = 0;
@@ -2027,7 +1901,7 @@ static void p_chart_simplify_compute(PChart *chart)
                p_collapse_cost_vertex(v, &cost, &e);
 
                if (e)
-                       v->u.heaplink = pheap_insert(heap, cost, e);
+                       v->u.heaplink = BLI_heap_insert(heap, cost, e);
                else
                        v->u.heaplink = NULL;
        }
@@ -2036,12 +1910,12 @@ static void p_chart_simplify_compute(PChart *chart)
                e->u.nextcollapse = NULL;
 
        /* pop edge collapse out of heap one by one */
-       while (!pheap_empty(heap)) {
+       while (!BLI_heap_empty(heap)) {
                if (ncollapsed == NCOLLAPSE)
                        break;
 
-               PHeapLink *link = pheap_toplink(heap);
-               PEdge *edge = (PEdge*)pheap_popmin(heap), *pair = edge->pair;
+               HeapNode *link = BLI_heap_top(heap);
+               PEdge *edge = (PEdge*)BLI_heap_popmin(heap), *pair = edge->pair;
                PVert *oldv, *keepv;
                PEdge *wheele, *nexte;
 
@@ -2085,21 +1959,21 @@ static void p_chart_simplify_compute(PChart *chart)
                        v = wheelverts[i];
 
                        if (v->u.heaplink) {
-                               pheap_remove(heap, 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 = pheap_insert(heap, cost, collapse);
+                               v->u.heaplink = BLI_heap_insert(heap, cost, collapse);
                }
 
                ncollapsed++;
        }
 
        MEM_freeN(wheelverts);
-       pheap_delete(heap);
+       BLI_heap_free(heap, NULL);
 
        p_chart_post_collapse_flush(chart, collapsededges);
 }
@@ -2136,7 +2010,6 @@ static void p_chart_complexify(PChart *chart)
 
        p_chart_post_split_flush(chart);
 }
-#endif
 
 #if 0
 static void p_chart_simplify(PChart *chart)
@@ -2144,6 +2017,7 @@ static void p_chart_simplify(PChart *chart)
        /* Not implemented, needs proper reordering in split_flush. */
 }
 #endif
+#endif
 
 /* ABF */
 
@@ -2327,19 +2201,19 @@ static float p_abf_compute_gradient(PAbfSystem *sys, PChart *chart)
 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);
-       /*nlSolverParameteri(NL_ABF, NL_TRUE); should set symmetric */
 
        nlBegin(NL_SYSTEM);
 
        nlBegin(NL_MATRIX);
 
        for (i = 0; i < nvar; i++)
-               nlRightHandSideAdd(i, sys->bInterior[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];
@@ -2385,8 +2259,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        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(v1->u.id, j2[0][0]*beta[0]);
-                       nlRightHandSideAdd(ninterior + v1->u.id, j2[1][0]*beta[1] + j2[2][0]*beta[2]);
+                       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];
@@ -2405,8 +2279,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        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(v2->u.id, j2[1][1]*beta[1]);
-                       nlRightHandSideAdd(ninterior + v2->u.id, j2[0][1]*beta[0] + j2[2][1]*beta[2]);
+                       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];
@@ -2425,8 +2299,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        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(v3->u.id, j2[2][2]*beta[2]);
-                       nlRightHandSideAdd(ninterior + v3->u.id, j2[0][2]*beta[0] + j2[1][2]*beta[1]);
+                       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];
@@ -2483,24 +2357,24 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        pre[0] = pre[1] = pre[2] = 0.0;
 
                        if (v1->flag & PVERT_INTERIOR) {
-                               float x = nlGetVariable(v1->u.id);
-                               float x2 = nlGetVariable(ninterior + v1->u.id);
+                               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(v2->u.id);
-                               float x2 = nlGetVariable(ninterior + v2->u.id);
+                               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(v3->u.id);
-                               float x2 = nlGetVariable(ninterior + v3->u.id);
+                               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;
@@ -2519,11 +2393,20 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 
                        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(i);
-                       sys->lambdaLength[i] += nlGetVariable(ninterior + i);
+                       sys->lambdaPlanar[i] += nlGetVariable(0, i);
+                       sys->lambdaLength[i] += nlGetVariable(0, ninterior + i);
                }
        }
 
@@ -2623,13 +2506,6 @@ static PBool p_chart_abf_solve(PChart *chart)
                for (i = 0; i < ABF_MAX_ITER; i++) {
                        float norm = p_abf_compute_gradient(&sys, chart);
 
-#if 0
-                       if (norm > lastnorm)
-                               printf("convergence error: %f => %f\n", lastnorm, norm);
-                       else
-                               printf("norm: %f\n", norm);
-#endif
-
                        lastnorm = norm;
 
                        if (norm < limit)
@@ -2659,7 +2535,162 @@ static PBool p_chart_abf_solve(PChart *chart)
 
 /* Least Squares Conformal Maps */
 
-static void p_chart_extrema_verts(PChart *chart, PVert **v1, PVert **v2)
+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];
@@ -2696,37 +2727,10 @@ static void p_chart_extrema_verts(PChart *chart, PVert **v1, PVert **v2)
                }
        }
 
-       if (minvert[dir] == maxvert[dir]) {
-               /* degenerate case */
-               PFace *f = chart->faces;
-               *v1 = f->edge->vert;
-               *v2 = f->edge->next->vert;
-
-               (*v1)->uv[0] = 0.0f;
-               (*v1)->uv[1] = 0.5f;
-               (*v2)->uv[0] = 1.0f;
-               (*v2)->uv[1] = 0.5f;
-       }
-       else {
-               int diru, dirv;
-
-               *v1 = minvert[dir];
-               *v2 = maxvert[dir];
-
-               if (dir == 2) {
-                       diru = 1;
-                       dirv = 0;
-               }
-               else {
-                       diru = 0;
-                       dirv = 1;
-               }
+       *pin1 = minvert[dir];
+       *pin2 = maxvert[dir];
 
-               (*v1)->uv[diru] = (*v1)->co[dir];
-               (*v1)->uv[dirv] = (*v1)->co[(dir+1)%3];
-               (*v2)->uv[diru] = (*v2)->co[dir];
-               (*v2)->uv[dirv] = (*v2)->co[(dir+1)%3];
-       }
+       p_chart_pin_positions(chart, pin1, pin2);
 }
 
 static void p_chart_lscm_load_solution(PChart *chart)
@@ -2734,8 +2738,8 @@ static void p_chart_lscm_load_solution(PChart *chart)
        PVert *v;
 
        for (v=chart->verts; v; v=v->nextlink) {
-               v->uv[0] = nlGetVariable(2*v->u.id);
-               v->uv[1] = nlGetVariable(2*v->u.id + 1);
+               v->uv[0] = nlGetVariable(0, 2*v->u.id);
+               v->uv[1] = nlGetVariable(0, 2*v->u.id + 1);
        }
 }
 
@@ -2767,15 +2771,18 @@ static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
 #endif
 
                if (abf) {
-                       PBool p_chart_abf_solve(PChart *chart);
-
                        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 */
-                       p_chart_extrema_verts(chart, &pin1, &pin2);
+                       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;
@@ -2789,6 +2796,7 @@ static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
 
                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();
@@ -2800,6 +2808,7 @@ 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);
 
@@ -2819,10 +2828,10 @@ static PBool p_chart_lscm_solve(PChart *chart)
                nlLockVariable(2*pin2->u.id);
                nlLockVariable(2*pin2->u.id + 1);
        
-               nlSetVariable(2*pin1->u.id, pin1->uv[0]);
-               nlSetVariable(2*pin1->u.id + 1, pin1->uv[1]);
-               nlSetVariable(2*pin2->u.id, pin2->uv[0]);
-               nlSetVariable(2*pin2->u.id + 1, pin2->uv[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 */
@@ -2831,8 +2840,8 @@ static PBool p_chart_lscm_solve(PChart *chart)
                                nlLockVariable(2*v->u.id);
                                nlLockVariable(2*v->u.id + 1);
 
-                               nlSetVariable(2*v->u.id, v->uv[0]);
-                               nlSetVariable(2*v->u.id + 1, v->uv[1]);
+                               nlSetVariable(0, 2*v->u.id, v->uv[0]);
+                               nlSetVariable(0, 2*v->u.id + 1, v->uv[1]);
                        }
                }
        }
@@ -2841,6 +2850,7 @@ static PBool p_chart_lscm_solve(PChart *chart)
 
        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;
@@ -2879,10 +2889,11 @@ static PBool p_chart_lscm_solve(PChart *chart)
                }
 
                /* angle based lscm formulation */
-               ratio = (sina3 == 0.0f)? 0.0f: sina2/sina3;
+               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);
@@ -2898,6 +2909,21 @@ static PBool p_chart_lscm_solve(PChart *chart)
                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);
@@ -2908,6 +2934,12 @@ static PBool p_chart_lscm_solve(PChart *chart)
                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;
 }
@@ -3073,144 +3105,6 @@ static void p_chart_stretch_minimize(PChart *chart, RNG *rng)
        }
 }
 
-/* Packing */
-
-static int p_compare_chart_area(const void *a, const void *b)
-{
-       PChart *ca = *((PChart**)a);
-       PChart *cb = *((PChart**)b);
-
-    if (ca->u.pack.area > cb->u.pack.area)
-               return -1;
-       else if (ca->u.pack.area == cb->u.pack.area)
-               return 0;
-       else
-               return 1;
-}
-
-static PBool p_pack_try(PHandle *handle, float side)
-{
-       PChart *chart;
-       float packx, packy, rowh, groupw, w, h;
-       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];
-
-               if (chart->flag & PCHART_NOPACK)
-                       continue;
-
-               w = chart->u.pack.size[0];
-               h = chart->u.pack.size[1];
-
-               if(w <= (side-packx)) {
-                       chart->u.pack.trans[0] = packx;
-                       chart->u.pack.trans[1] = packy;
-
-                       packx += w;
-                       rowh= MAX2(rowh, h);
-               }
-               else {
-                       packy += rowh;
-                       packx = w;
-                       rowh = h;
-
-                       chart->u.pack.trans[0] = 0.0;
-                       chart->u.pack.trans[1] = packy;
-               }
-
-               if (packy+rowh > side)
-                       return P_FALSE;
-       }
-
-       return P_TRUE;
-}
-
-#define PACK_SEARCH_DEPTH 15
-
-void p_charts_pack(PHandle *handle)
-{
-       PChart *chart;
-       float uv_area, area, trans[2], minside, maxside, totarea, side;
-       int i;
-
-       /* very simple rectangle packing */
-
-       if (handle->ncharts == 0)
-               return;
-
-       totarea = 0.0f;
-       maxside = 0.0f;
-
-       for (i = 0; i < handle->ncharts; i++) {
-               chart = handle->charts[i];
-
-               if (chart->flag & PCHART_NOPACK) {
-                       chart->u.pack.area = 0.0f;
-                       continue;
-               }
-
-               p_chart_area(chart, &uv_area, &area);
-               p_chart_uv_bbox(chart, trans, chart->u.pack.size);
-
-               /* translate to origin and make area equal to 3d area */
-               chart->u.pack.rescale = (uv_area > 0.0f)? sqrt(area)/sqrt(uv_area): 0.0f;
-               chart->u.pack.area = area;
-               totarea += area;
-
-               trans[0] = -trans[0];
-               trans[1] = -trans[1];
-               p_chart_uv_translate(chart, trans);
-               p_chart_uv_scale(chart, chart->u.pack.rescale);
-
-               /* compute new dimensions for packing */
-               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]);
-       }
-
-       /* sort by chart area, largest first */
-       qsort(handle->charts, handle->ncharts, sizeof(PChart*), p_compare_chart_area);
-
-       /* binary search over pack region size */
-       minside = MAX2(sqrt(totarea), maxside);
-       maxside = (((int)sqrt(handle->ncharts-1))+1)*maxside;
-
-       if (minside < maxside) { /* should always be true */
-
-               for (i = 0; i < PACK_SEARCH_DEPTH; i++) {
-                       if (p_pack_try(handle, (minside+maxside)*0.5f + 1e-5))
-                               maxside = (minside+maxside)*0.5f;
-                       else
-                               minside = (minside+maxside)*0.5f;
-               }
-       }
-
-       /* do the actual packing */
-       side = maxside + 1e-5;
-       if (!p_pack_try(handle, side))
-               param_warning("packing failed.\n");
-
-       for (i = 0; i < handle->ncharts; i++) {
-               chart = handle->charts[i];
-
-               if (chart->flag & PCHART_NOPACK)
-                       continue;
-
-               p_chart_uv_scale(chart, 1.0f/side);
-               trans[0] = chart->u.pack.trans[0]/side;
-               trans[1] = chart->u.pack.trans[1]/side;
-               p_chart_uv_translate(chart, trans);
-       }
-}
-
 /* Minimum area enclosing rectangle for packing */
 
 static int p_compare_geometric_uv(const void *a, const void *b)
@@ -4004,7 +3898,7 @@ void param_face_add(ParamHandle *handle, ParamKey key, int nverts,
        param_assert((nverts == 3) || (nverts == 4));
 
        if (nverts == 4) {
-               if (!p_quad_split_direction(co)) {
+               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);
                }
@@ -4106,8 +4000,8 @@ void param_lscm_solve(ParamHandle *handle)
                if (chart->u.lscm.context) {
                        result = p_chart_lscm_solve(chart);
 
-                       /*if (result)
-                               p_chart_rotate_minimum_area(chart);*/
+                       if (result && !(chart->flag & PCHART_NOPACK))
+                               p_chart_rotate_minimum_area(chart);
 
                        if (!result || (chart->u.lscm.pin1))
                                p_chart_lscm_end(chart);
@@ -4211,10 +4105,66 @@ void param_smooth_area(ParamHandle *handle)
                p_smooth(chart);
        }
 }
-
 void param_pack(ParamHandle *handle)
 {
-       p_charts_pack((PHandle*)handle);
+       /* box packing variables */
+       boxPack *boxarray, *box;
+       float tot_width, tot_height, scale;
+        
+       PChart *chart;
+       int i, unpacked=0;
+       float trans[2];
+       
+       PHandle *phandle = (PHandle*)handle;
+       
+       
+       if (phandle->ncharts == 0)
+               return;
+       
+       /* 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];
+               
+               
+               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);
+               
+               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);
 }
 
 void param_flush(ParamHandle *handle)