2.5: Space Image ported back

Organized as follows:

uvedit/
	uv editing related code
	uvedit_draw.c: drawing code
	uvedit_ops.c: operators, just a few done
	uvedit_unwrap_ops.c: will be operators for unwrapping
	uvedit_paramatrizer.c: lscm/abf/stretch/pack

space_image/
	space_image.c: registration and common getter/setters
	image_draw.c: drawing code, mostly functional
	image_panels.c: panels, all commented out
	image_render.c: render callbacks, non functional
	image_ops.c: operators, only view navigation done
	image_header.c: header, menus mostly done but missing buttons

Notes:
* Header menus consist only of Operator and RNA buttons, if they
  are not implemented they're displayed grayed out. Ideally the full
  header could work like this, but std_libbuttons looks problematic.

* Started using view2d code more than the old code, but for now it
  still does own view2d management due to some very specific
  requirements that the image window has. The drawing code however
  is more clear hopefully, it only uses view2d, and there is no
  switching between 'p' and 'f' view2d's anymore, it is always 'f'.

* In order to make uvedit operators more independent I move some
  image space settings to scene toolsettings, and the current image
  and its buffer is in the context. Especially sync selection and
  select mode belonged there anyway as this cannot work correct with
  different spaces having different settings anyway.

* Image paint is not back yet, did not want to put that together with
  uvedit because there's really no code sharing.. perhaps vertex paint,
  image paint and sculpt would be good to have in one module to share
  brush code, partial redraw, etc better.

1 files changed, 4481 insertions, 0 deletions

diff --git a/source/blender/editors/uvedit/uvedit_parametrizer.c b/source/blender/editors/uvedit/uvedit_parametrizer.c
new file mode 100644
index 00000000000..635a7903d6f
--- /dev/null
+++ b/source/blender/editors/uvedit/uvedit_parametrizer.c
@@ -0,0 +1,4481 @@
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_memarena.h"
+#include "BLI_arithb.h"
+#include "BLI_rand.h"
+#include "BLI_heap.h"
+#include "BLI_boxpack2d.h"
+
+#include "BKE_utildefines.h"
+
+#include "ONL_opennl.h"
+
+#include "uvedit_intern.h"
+#include "uvedit_parametrizer.h"
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "BLO_sys_types.h" // for intptr_t support
+
+#if defined(_WIN32)
+#define M_PI 3.14159265358979323846
+#endif
+
+/* Utils */
+
+#if 0
+	#define param_assert(condition);
+	#define param_warning(message);
+	#define param_test_equals_ptr(condition);
+	#define param_test_equals_int(condition);
+#else
+	#define param_assert(condition) \
+		if (!(condition)) \
+			{ /*printf("Assertion %s:%d\n", __FILE__, __LINE__); abort();*/ }
+	#define param_warning(message) \
+		{ /*printf("Warning %s:%d: %s\n", __FILE__, __LINE__, message);*/ }
+	#define param_test_equals_ptr(str, a, b) \
+		if (a != b) \
+			{ /*printf("Equals %s => %p != %p\n", str, a, b);*/ };
+	#define param_test_equals_int(str, a, b) \
+		if (a != b) \
+			{ /*printf("Equals %s => %d != %d\n", str, a, b);*/ };
+#endif
+
+typedef enum PBool {
+	P_TRUE = 1,
+	P_FALSE = 0
+} PBool;
+
+/* Special Purpose Hash */
+
+typedef intptr_t PHashKey;
+
+typedef struct PHashLink {
+	struct PHashLink *next;
+	PHashKey key;
+} PHashLink;
+
+typedef struct PHash {
+	PHashLink **list;
+	PHashLink **buckets;
+	int size, cursize, cursize_id;
+} PHash;
+
+
+
+struct PVert;
+struct PEdge;
+struct PFace;
+struct PChart;
+struct PHandle;
+
+/* Simplices */
+
+typedef struct PVert {
+	struct PVert *nextlink;
+
+	union PVertUnion {
+		PHashKey key;			/* construct */
+		int id;					/* abf/lscm matrix index */
+		float distortion;		/* area smoothing */
+		HeapNode *heaplink;		/* edge collapsing */
+	} u;
+
+	struct PEdge *edge;
+	float *co;
+	float uv[2];
+	unsigned char flag;
+
+} PVert; 
+
+typedef struct PEdge {
+	struct PEdge *nextlink;
+
+	union PEdgeUnion {
+		PHashKey key;					/* construct */
+		int id;							/* abf matrix index */
+		HeapNode *heaplink;				/* fill holes */
+		struct PEdge *nextcollapse;		/* simplification */
+	} u;
+
+	struct PVert *vert;
+	struct PEdge *pair;
+	struct PEdge *next;
+	struct PFace *face;
+	float *orig_uv, old_uv[2];
+	unsigned short flag;
+
+} PEdge;
+
+typedef struct PFace {
+	struct PFace *nextlink;
+
+	union PFaceUnion {
+		PHashKey key;			/* construct */
+		int chart;				/* construct splitting*/
+		float area3d;			/* stretch */
+		int id;					/* abf matrix index */
+	} u;
+
+	struct PEdge *edge;
+	unsigned char flag;
+
+} PFace;
+
+enum PVertFlag {
+	PVERT_PIN = 1,
+	PVERT_SELECT = 2,
+	PVERT_INTERIOR = 4,
+	PVERT_COLLAPSE = 8,
+	PVERT_SPLIT = 16
+};
+
+enum PEdgeFlag {
+	PEDGE_SEAM = 1,
+	PEDGE_VERTEX_SPLIT = 2,
+	PEDGE_PIN = 4,
+	PEDGE_SELECT = 8,
+	PEDGE_DONE = 16,
+	PEDGE_FILLED = 32,
+	PEDGE_COLLAPSE = 64,
+	PEDGE_COLLAPSE_EDGE = 128,
+	PEDGE_COLLAPSE_PAIR = 256
+};
+
+/* for flipping faces */
+#define PEDGE_VERTEX_FLAGS (PEDGE_PIN)
+
+enum PFaceFlag {
+	PFACE_CONNECTED = 1,
+	PFACE_FILLED = 2,
+	PFACE_COLLAPSE = 4
+};
+
+/* Chart */
+
+typedef struct PChart {
+	PVert *verts;
+	PEdge *edges;
+	PFace *faces;
+	int nverts, nedges, nfaces;
+
+	PVert *collapsed_verts;
+	PEdge *collapsed_edges;
+	PFace *collapsed_faces;
+
+	union PChartUnion {
+		struct PChartLscm {
+			NLContext context;
+			float *abf_alpha;
+			PVert *pin1, *pin2;
+		} lscm;
+		struct PChartPack {
+			float rescale, area;
+			float size[2], trans[2];
+		} pack;
+	} u;
+
+	unsigned char flag;
+	struct PHandle *handle;
+} PChart;
+
+enum PChartFlag {
+	PCHART_NOPACK = 1
+};
+
+enum PHandleState {
+	PHANDLE_STATE_ALLOCATED,
+	PHANDLE_STATE_CONSTRUCTED,
+	PHANDLE_STATE_LSCM,
+	PHANDLE_STATE_STRETCH
+};
+
+typedef struct PHandle {
+	enum PHandleState state;
+	MemArena *arena;
+
+	PChart *construction_chart;
+	PHash *hash_verts;
+	PHash *hash_edges;
+	PHash *hash_faces;
+
+	PChart **charts;
+	int ncharts;
+
+	float aspx, aspy;
+
+	RNG *rng;
+	float blend;
+} PHandle;
+
+
+/* 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, 
+	16411, 32771, 65537, 131101, 262147, 524309, 1048583, 2097169, 
+	4194319, 8388617, 16777259, 33554467, 67108879, 134217757, 268435459
+};
+
+#define PHASH_hash(ph, item) (((uintptr_t) (item))%((unsigned int) (ph)->cursize))
+#define PHASH_edge(v1, v2)	 ((v1)^(v2))
+
+static PHash *phash_new(PHashLink **list, int sizehint)
+{
+	PHash *ph = (PHash*)MEM_callocN(sizeof(PHash), "PHash");
+	ph->size = 0;
+	ph->cursize_id = 0;
+	ph->list = list;
+
+	while (PHashSizes[ph->cursize_id] < sizehint)
+		ph->cursize_id++;
+
+	ph->cursize = PHashSizes[ph->cursize_id];
+	ph->buckets = (PHashLink**)MEM_callocN(ph->cursize*sizeof(*ph->buckets), "PHashBuckets");
+
+	return ph;
+}
+
+static void phash_delete(PHash *ph)
+{
+	MEM_freeN(ph->buckets);
+	MEM_freeN(ph);
+}
+
+static int phash_size(PHash *ph)
+{
+	return ph->size;
+}
+
+static void phash_insert(PHash *ph, PHashLink *link)
+{
+	int size = ph->cursize;
+	int hash = PHASH_hash(ph, link->key);
+	PHashLink *lookup = ph->buckets[hash];
+
+	if (lookup == NULL) {
+		/* insert in front of the list */
+		ph->buckets[hash] = link;
+		link->next = *(ph->list);
+		*(ph->list) = link;
+	}
+	else {
+		/* insert after existing element */
+		link->next = lookup->next;
+		lookup->next = link;
+	}
+		
+	ph->size++;
+
+	if (ph->size > (size*3)) {
+		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->list) = NULL;
+
+		for (link = first; link; link = next) {
+			next = link->next;
+			phash_insert(ph, link);
+		}
+	}
+}
+
+static PHashLink *phash_lookup(PHash *ph, PHashKey key)
+{
+	PHashLink *link;
+	int hash = PHASH_hash(ph, key);
+
+	for (link = ph->buckets[hash]; link; link = link->next)
+		if (link->key == key)
+			return link;
+		else if (PHASH_hash(ph, link->key) != hash)
+			return NULL;
+	
+	return link;
+}
+
+static PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
+{
+	int hash = PHASH_hash(ph, key);
+
+	for (link = link->next; link; link = link->next)
+		if (link->key == key)
+			return link;
+		else if (PHASH_hash(ph, link->key) != hash)
+			return NULL;
+	
+	return link;
+}
+
+/* Geometry */
+
+static float p_vec_angle_cos(float *v1, float *v2, float *v3)
+{
+	float d1[3], d2[3];
+
+	d1[0] = v1[0] - v2[0];
+	d1[1] = v1[1] - v2[1];
+	d1[2] = v1[2] - v2[2];
+
+	d2[0] = v3[0] - v2[0];
+	d2[1] = v3[1] - v2[1];
+	d2[2] = v3[2] - v2[2];
+
+	Normalize(d1);
+	Normalize(d2);
+
+	return d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2];
+}
+
+static float p_vec_angle(float *v1, float *v2, float *v3)
+{
+	float dot = p_vec_angle_cos(v1, v2, v3);
+
+	if (dot <= -1.0f)
+		return (float)M_PI;
+	else if (dot >= 1.0f)
+		return 0.0f;
+	else
+		return (float)acos(dot);
+}
+
+static float p_vec2_angle(float *v1, float *v2, float *v3)
+{
+	float u1[3], u2[3], u3[3];
+
+	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 p_vec_angle(u1, u2, u3);
+}
+
+static void p_triangle_angles(float *v1, float *v2, float *v3, float *a1, float *a2, float *a3)
+{
+	*a1 = p_vec_angle(v3, v1, v2);
+	*a2 = p_vec_angle(v1, v2, v3);
+	*a3 = M_PI - *a2 - *a1;
+}
+
+static void p_face_angles(PFace *f, float *a1, float *a2, float *a3)
+{
+	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);
+}
+
+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;
+
+	return AreaT3Dfl(v1->co, v2->co, v3->co);
+}
+
+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])));
+}
+
+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;
+
+	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])));
+}
+
+static float p_edge_length(PEdge *e)
+{
+    PVert *v1 = e->vert, *v2 = e->next->vert;
+    float d[3];
+
+    d[0] = v2->co[0] - v1->co[0];
+    d[1] = v2->co[1] - v1->co[1];
+    d[2] = v2->co[2] - v1->co[2];
+
+    return sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
+}
+
+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]);
+}
+
+static void p_chart_uv_bbox(PChart *chart, float *minv, float *maxv)
+{
+	PVert *v;
+
+	INIT_MINMAX2(minv, maxv);
+
+	for (v=chart->verts; v; v=v->nextlink) {
+		DO_MINMAX2(v->uv, minv, maxv);
+	}
+}
+
+static void p_chart_uv_scale(PChart *chart, float scale)
+{
+	PVert *v;
+
+	for (v=chart->verts; v; v=v->nextlink) {
+		v->uv[0] *= scale;
+		v->uv[1] *= scale;
+	}
+}
+
+static void p_chart_uv_scale_xy(PChart *chart, float x, float y)
+{
+	PVert *v;
+
+	for (v=chart->verts; v; v=v->nextlink) {
+		v->uv[0] *= x;
+		v->uv[1] *= y;
+	}
+}
+
+static void p_chart_uv_translate(PChart *chart, float trans[2])
+{
+	PVert *v;
+
+	for (v=chart->verts; v; v=v->nextlink) {
+		v->uv[0] += trans[0];
+		v->uv[1] += trans[1];
+	}
+}
+
+static PBool p_intersect_line_2d_dir(float *v1, float *dir1, float *v2, float *dir2, float *isect)
+{
+	float lmbda, div;
+
+	div= dir2[0]*dir1[1] - dir2[1]*dir1[0];
+
+	if (div == 0.0f)
+		return P_FALSE;
+
+	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 P_TRUE;
+}
+
+#if 0
+static PBool p_intersect_line_2d(float *v1, float *v2, float *v3, float *v4, float *isect)
+{
+	float dir1[2], dir2[2];
+
+	dir1[0] = v4[0] - v3[0];
+	dir1[1] = v4[1] - v3[1];
+
+	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;
+	}
+
+	return P_TRUE;
+}
+#endif
+
+/* Topological Utilities */
+
+static PEdge *p_wheel_edge_next(PEdge *e)
+{
+	return e->next->next->pair;
+}
+
+static PEdge *p_wheel_edge_prev(PEdge *e)
+{   
+	return (e->pair)? e->pair->next: NULL;
+}
+
+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 PBool p_vert_interior(PVert *v)
+{
+	return (v->edge->pair != NULL);
+}
+
+static void p_face_flip(PFace *f)
+{
+	PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+	PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+	int f1 = e1->flag, f2 = e2->flag, f3 = e3->flag;
+
+	e1->vert = v2;
+	e1->next = e3;
+	e1->flag = (f1 & ~PEDGE_VERTEX_FLAGS) | (f2 & PEDGE_VERTEX_FLAGS);
+
+	e2->vert = v3;
+	e2->next = e1;
+	e2->flag = (f2 & ~PEDGE_VERTEX_FLAGS) | (f3 & PEDGE_VERTEX_FLAGS);
+
+	e3->vert = v1;
+	e3->next = e2;
+	e3->flag = (f3 & ~PEDGE_VERTEX_FLAGS) | (f1 & PEDGE_VERTEX_FLAGS);
+}
+
+#if 0
+static void p_chart_topological_sanity_check(PChart *chart)
+{
+	PVert *v;
+	PEdge *e;
+
+	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);
+		}
+	}
+}
+#endif
+
+/* Loading / Flushing */
+
+static void p_vert_load_pin_select_uvs(PHandle *handle, PVert *v)
+{
+	PEdge *e;
+	int nedges = 0, npins = 0;
+	float pinuv[2];
+
+	v->uv[0] = v->uv[1] = 0.0f;
+	pinuv[0] = pinuv[1] = 0.0f;
+	e = v->edge;
+	do {
+		if (e->orig_uv) {
+			if (e->flag & PEDGE_SELECT)
+				v->flag |= PVERT_SELECT;
+
+ 			if (e->flag & PEDGE_PIN) {
+				pinuv[0] += e->orig_uv[0]*handle->aspx;
+				pinuv[1] += e->orig_uv[1]*handle->aspy;
+				npins++;
+			}
+			else {
+				v->uv[0] += e->orig_uv[0]*handle->aspx;
+				v->uv[1] += e->orig_uv[1]*handle->aspy;
+			}
+
+			nedges++;
+		}
+
+		e = p_wheel_edge_next(e);
+	} while (e && e != (v->edge));
+
+	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_flush_uvs(PHandle *handle, PChart *chart)
+{
+	PEdge *e;
+
+	for (e=chart->edges; e; e=e->nextlink) {
+		if (e->orig_uv) {
+			e->orig_uv[0] = e->vert->uv[0]/handle->aspx;
+			e->orig_uv[1] = e->vert->uv[1]/handle->aspy;
+		}
+	}
+}
+
+static void p_flush_uvs_blend(PHandle *handle, PChart *chart, float blend)
+{
+	PEdge *e;
+	float invblend = 1.0f - blend;
+
+	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]/handle->aspx;
+			e->orig_uv[1] = blend*e->old_uv[1] + invblend*e->vert->uv[1]/handle->aspy;
+		}
+	}
+}
+
+static void p_face_backup_uvs(PFace *f)
+{
+	PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+
+	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];
+	}
+}
+
+static void p_face_restore_uvs(PFace *f)
+{
+	PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+
+	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];
+	}
+}
+
+/* Construction (use only during construction, relies on u.key being set */
+
+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;
+
+	phash_insert(handle->hash_verts, (PHashLink*)v);
+
+	return v;
+}
+
+static PVert *p_vert_lookup(PHandle *handle, PHashKey key, float *co, PEdge *e)
+{
+	PVert *v = (PVert*)phash_lookup(handle->hash_verts, key);
+
+	if (v)
+		return v;
+	else
+		return p_vert_add(handle, key, co, e);
+}
+
+static PVert *p_vert_copy(PChart *chart, PVert *v)
+{
+	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->u.key = v->u.key;
+	nv->edge = v->edge;
+	nv->flag = v->flag;
+
+	return nv;
+}
+
+static PEdge *p_edge_lookup(PHandle *handle, PHashKey *vkeys)
+{
+	PHashKey key = PHASH_edge(vkeys[0], vkeys[1]);
+	PEdge *e = (PEdge*)phash_lookup(handle->hash_edges, key);
+
+	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;
+
+		e = (PEdge*)phash_next(handle->hash_edges, key, (PHashLink*)e);
+	}
+
+	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");
+	chart->handle = handle;
+
+	return chart;
+}
+
+static void p_chart_delete(PChart *chart)
+{
+	/* the actual links are free by memarena */
+	MEM_freeN(chart);
+}
+
+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->u.key == ep->vert->u.key) {
+		uvp1 = ep->orig_uv;
+		uvp2 = ep->next->orig_uv;
+	}
+	else {
+		uvp1 = ep->next->orig_uv;
+		uvp2 = ep->orig_uv;
+	}
+
+	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])) {
+		e->flag |= PEDGE_SEAM;
+		ep->flag |= PEDGE_SEAM;
+		return P_TRUE;
+	}
+	
+	return P_FALSE;
+}
+
+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->u.key;
+	PHashKey key2 = e->next->vert->u.key;
+
+	if (e->flag & PEDGE_SEAM)
+		return P_FALSE;
+	
+	key = PHASH_edge(key1, key2);
+	pe = (PEdge*)phash_lookup(handle->hash_edges, key);
+	*pair = NULL;
+
+	while (pe) {
+		if (pe != e) {
+			v1 = pe->vert;
+			v2 = pe->next->vert;
+
+			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 ||
+				    (impl && p_edge_implicit_seam(e, pe))) {
+					*pair = NULL;
+					return P_FALSE;
+				}
+
+				*pair = pe;
+			}
+		}
+
+		pe = (PEdge*)phash_next(handle->hash_edges, key, (PHashLink*)pe);
+	}
+
+	if (*pair && (e->vert == (*pair)->vert)) {
+		if ((*pair)->next->pair || (*pair)->next->next->pair) {
+			/* non unfoldable, maybe mobius ring or klein bottle */
+			*pair = NULL;
+			return P_FALSE;
+		}
+	}
+
+	return (*pair != NULL);
+}
+
+static PBool p_edge_connect_pair(PHandle *handle, PEdge *e, PEdge ***stack, PBool impl)
+{
+	PEdge *pair = NULL;
+
+	if(!e->pair && p_edge_has_pair(handle, e, &pair, impl)) {
+		if (e->vert == pair->vert)
+			p_face_flip(pair->face);
+
+		e->pair = pair;
+		pair->pair = e;
+
+		if (!(pair->face->flag & PFACE_CONNECTED)) {
+			**stack = pair;
+			(*stack)++;
+		}
+	}
+
+	return (e->pair != NULL);
+}
+
+static int p_connect_pairs(PHandle *handle, PBool impl)
+{
+	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=chart->faces; first; first=first->nextlink) {
+		if (first->flag & PFACE_CONNECTED)
+			continue;
+
+		*stack = first->edge;
+		stack++;
+
+		while (stack != stackbase) {
+			stack--;
+			e = *stack;
+			e1 = e->next;
+			e2 = e1->next;
+
+			f = e->face;
+			f->flag |= PFACE_CONNECTED;
+
+			/* assign verts to charts so we can sort them later */
+			f->u.chart = ncharts;
+
+			if (!p_edge_connect_pair(handle, e, &stack, impl))
+				e->vert->edge = e;
+			if (!p_edge_connect_pair(handle, e1, &stack, impl))
+				e1->vert->edge = e1;
+			if (!p_edge_connect_pair(handle, e2, &stack, impl))
+				e2->vert->edge = e2;
+		}
+
+		ncharts++;
+	}
+
+	MEM_freeN(stackbase);
+
+	return ncharts;
+}
+
+static void p_split_vert(PChart *chart, PEdge *e)
+{
+	PEdge *we, *lastwe = NULL;
+	PVert *v = e->vert;
+	PBool copy = P_TRUE;
+
+	if (e->flag & PEDGE_VERTEX_SPLIT)
+		return;
+
+	/* rewind to start */
+	lastwe = e;
+	for (we = p_wheel_edge_prev(e); we && (we != e); we = p_wheel_edge_prev(we))
+		lastwe = we;
+	
+	/* go over all edges in wheel */
+	for (we = lastwe; we; we = p_wheel_edge_next(we)) {
+		if (we->flag & PEDGE_VERTEX_SPLIT)
+			break;
+
+		we->flag |= PEDGE_VERTEX_SPLIT;
+
+		if (we == v->edge) {
+			/* found it, no need to copy */
+			copy = P_FALSE;
+			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;
+		do {
+			we->vert = v;
+			we = p_wheel_edge_next(we);
+		} while (we && (we != lastwe));
+	}
+}
+
+static PChart **p_split_charts(PHandle *handle, PChart *chart, int ncharts)
+{
+	PChart **charts = MEM_mallocN(sizeof*charts * ncharts, "PCharts"), *nchart;
+	PFace *f, *nextf;
+	int i;
+
+	for (i = 0; i < ncharts; i++)
+		charts[i] = p_chart_new(handle);
+
+	f = chart->faces;
+	while (f) {
+		PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+		nextf = f->nextlink;
+
+		nchart = charts[f->u.chart];
+
+		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);
+		p_split_vert(nchart, e3);
+
+		f = nextf;
+	}
+
+	return charts;
+}
+
+static PFace *p_face_add(PHandle *handle)
+{
+	PFace *f;
+	PEdge *e1, *e2, *e3;
+
+	/* allocate */
+	f = (PFace*)BLI_memarena_alloc(handle->arena, sizeof *f);
+	f->flag=0; // init !
+
+	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;
+	e1->face = e2->face = e3->face = f;
+
+	e1->next = e2;
+	e2->next = e3;
+	e3->next = e1;
+
+	e1->pair = NULL;
+	e2->pair = NULL;
+	e3->pair = NULL;
+   
+	e1->flag =0;
+	e2->flag =0;
+	e3->flag =0;
+
+	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;
+		if (pin[i2]) e2->flag |= PEDGE_PIN;
+		if (pin[i3]) e3->flag |= PEDGE_PIN;
+	}
+
+	if (select) {
+		if (select[i1]) e1->flag |= PEDGE_SELECT;
+		if (select[i2]) e2->flag |= PEDGE_SELECT;
+		if (select[i3]) e3->flag |= PEDGE_SELECT;
+	}
+
+	/* insert into hash */
+	f->u.key = key;
+	phash_insert(handle->hash_faces, (PHashLink*)f);
+
+	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(handle->hash_edges, (PHashLink*)e1);
+	phash_insert(handle->hash_edges, (PHashLink*)e2);
+	phash_insert(handle->hash_edges, (PHashLink*)e3);
+
+	return f;
+}
+
+static PFace *p_face_add_fill(PChart *chart, PVert *v1, PVert *v2, PVert *v3)
+{
+	PFace *f = p_face_add(chart->handle);
+	PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+
+	e1->vert = v1;
+	e2->vert = v2;
+	e3->vert = v3;
+
+	e1->orig_uv = e2->orig_uv = e3->orig_uv = NULL;
+
+	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;
+
+	chart->nfaces++;
+	chart->nedges += 3;
+
+	return f;
+}
+
+static PBool p_quad_split_direction(PHandle *handle, float **co, PHashKey *vkeys)
+{
+	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;
+	}
+	else {
+		if (p_face_exists(handle,vkeys,0,1,3) || p_face_exists(handle,vkeys,1,2,3))
+			return !dir;
+	}
+
+	return dir;
+}
+
+/* Construction: boundary filling */
+
+static void p_chart_boundaries(PChart *chart, int *nboundaries, PEdge **outer)
+{   
+    PEdge *e, *be;
+    float len, maxlen = -1.0;
+
+	if (nboundaries)
+		*nboundaries = 0;
+	if (outer)
+		*outer = NULL;
+
+	for (e=chart->edges; e; e=e->nextlink) {
+        if (e->pair || (e->flag & PEDGE_DONE))
+            continue;
+
+		if (nboundaries)
+			(*nboundaries)++;
+
+        len = 0.0f;
+    
+		be = e;
+		do {
+            be->flag |= PEDGE_DONE;
+            len += p_edge_length(be);
+			be = be->next->vert->edge;
+        } while(be != e);
+
+        if (outer && (len > maxlen)) {
+			*outer = e;
+            maxlen = len;
+        }
+    }
+
+	for (e=chart->edges; e; e=e->nextlink)
+        e->flag &= ~PEDGE_DONE;
+}
+
+static float p_edge_boundary_angle(PEdge *e)
+{
+	PEdge *we;
+	PVert *v, *v1, *v2;
+	float angle;
+	int n = 0;
+
+	v = e->vert;
+
+	/* concave angle check -- could be better */
+	angle = M_PI;
+
+	we = v->edge;
+	do {
+		v1 = we->next->vert;
+		v2 = we->next->next->vert;
+		angle -= p_vec_angle(v1->co, v->co, v2->co);
+
+		we = we->next->next->pair;
+		n++;
+	} while (we && (we != v->edge));
+
+	return angle;
+}
+
+static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
+{
+	PEdge *e, *e1, *e2;
+
+	PFace *f;
+	struct Heap *heap = BLI_heap_new();
+	float angle;
+
+	e = be;
+	do {
+		angle = p_edge_boundary_angle(e);
+		e->u.heaplink = BLI_heap_insert(heap, angle, e);
+
+		e = p_boundary_edge_next(e);
+	} while(e != be);
+
+	if (nedges == 2) {
+		/* no real boundary, but an isolated seam */
+		e = be->next->vert->edge;
+		e->pair = be;
+		be->pair = e;
+
+		BLI_heap_remove(heap, e->u.heaplink);
+		BLI_heap_remove(heap, be->u.heaplink);
+	}
+	else {
+		while (nedges > 2) {
+			PEdge *ne, *ne1, *ne2;
+
+			e = (PEdge*)BLI_heap_popmin(heap);
+
+			e1 = p_boundary_edge_prev(e);
+			e2 = p_boundary_edge_next(e);
+
+			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;
+
+
+
+
+
+			f = p_face_add_fill(chart, e->vert, e1->vert, e2->vert);
+			f->flag |= PFACE_FILLED;
+
+			ne = f->edge->next->next;
+			ne1 = f->edge;
+			ne2 = f->edge->next;
+
+			ne->flag = ne1->flag = ne2->flag = PEDGE_FILLED;
+
+			e->pair = ne;
+			ne->pair = e;
+			e1->pair = ne1;
+			ne1->pair = e1;
+
+			ne->vert = e2->vert;
+			ne1->vert = e->vert;
+			ne2->vert = e1->vert;
+
+			if (nedges == 3) {
+				e2->pair = ne2;
+				ne2->pair = e2;
+			}
+			else {
+				ne2->vert->edge = ne2;
+				
+				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--;
+		}
+	}
+
+	BLI_heap_free(heap, NULL);
+}
+
+static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
+{
+	PEdge *e, *be; /* *enext - as yet unused */
+	int nedges;
+
+	for (e=chart->edges; e; e=e->nextlink) {
+		/* enext = e->nextlink; - as yet unused */
+
+        if (e->pair || (e->flag & PEDGE_FILLED))
+            continue;
+
+		nedges = 0;
+		be = e;
+		do {
+			be->flag |= PEDGE_FILLED;
+			be = be->next->vert->edge;
+			nedges++;
+		} while(be != e);
+
+		if (e != outer)
+			p_chart_fill_boundary(chart, e, nedges);
+    }
+}
+
+#if 0
+/* Polygon kernel for inserting uv's non overlapping */
+
+static int p_polygon_point_in(float *cp1, float *cp2, float *p)
+{
+	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);
+}
+
+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;
+	}
+}
+
+static void p_polygon_kernel_center(float (*points)[2], int npoints, float *center)
+{
+	int i, size, nnewpoints = npoints;
+	float (*oldpoints)[2], (*newpoints)[2], *p1, *p2;
+	
+	size = npoints*3;
+	oldpoints = MEM_mallocN(sizeof(float)*2*size, "PPolygonOldPoints");
+	newpoints = MEM_mallocN(sizeof(float)*2*size, "PPolygonNewPoints");
+
+	memcpy(oldpoints, points, sizeof(float)*2*npoints);
+
+	for (i = 0; i < npoints; i++) {
+		p1 = points[i];
+		p2 = points[(i+1)%npoints];
+		p_polygon_kernel_clip(oldpoints, nnewpoints, newpoints, &nnewpoints, p1, p2);
+
+		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];
+
+			MEM_freeN(oldpoints);
+			MEM_freeN(newpoints);
+
+			return;
+		}
+
+		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 {
+			float (*sw_points)[2] = oldpoints;
+			oldpoints = newpoints;
+			newpoints = sw_points;
+		}
+	}
+
+	center[0] = center[1] = 0.0f;
+
+	for (i = 0; i < nnewpoints; i++) {
+		center[0] += oldpoints[i][0];
+		center[1] += oldpoints[i][1];
+	}
+
+	center[0] /= nnewpoints;
+	center[1] /= nnewpoints;
+
+	MEM_freeN(oldpoints);
+	MEM_freeN(newpoints);
+}
+#endif
+
+#if 0
+/* Edge Collapser */
+
+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;
+
+	VecSubf(a, v2, v1);
+	VecSubf(b, v3, v1);
+	Crossf(c, a, b);
+
+	clen = VecLength(c);
+
+	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;
+
+	do {
+		if (p_face_uv_area_signed(e->face) < 0.0f)
+			return P_TRUE;
+
+		e = p_wheel_edge_next(e);
+	} while (e && (e != v->edge));
+
+	return P_FALSE;
+}
+
+static PBool p_vert_map_harmonic_weights(PVert *v)
+{
+	float weightsum, positionsum[2], olduv[2];
+
+	weightsum = 0.0f;
+	positionsum[0] = positionsum[1] = 0.0f;
+
+	if (p_vert_interior(v)) {
+		PEdge *e = v->edge;
+
+		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;
+
+		do {
+			float t1, t2;
+			PVert *v1, *v2;
+
+			v2 = e->next->vert;
+			v1 = e->next->next->vert;
+
+			t1 = p_vert_cotan(v1->co, v->co, v2->co);
+			t2 = p_vert_cotan(v2->co, v->co, v1->co);
+
+			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));
+	}
+
+	if (weightsum != 0.0f) {
+		weightsum = 1.0f/weightsum;
+		positionsum[0] *= weightsum;
+		positionsum[1] *= weightsum;
+	}
+
+	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];
+
+		return P_FALSE;
+	}
+
+	return P_TRUE;
+}
+
+static void p_vert_harmonic_insert(PVert *v)
+{
+	PEdge *e;
+
+	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 */
+
+		int npoints = 0, i;
+		float (*points)[2];
+
+		e = v->edge;
+		do {
+			npoints++;	
+			e = p_wheel_edge_next(e);
+		} while (e && (e != v->edge));
+
+		if (e == NULL)
+			npoints++;
+
+		points = MEM_mallocN(sizeof(float)*2*npoints, "PHarmonicPoints");
+
+		e = v->edge;
+		i = 0;
+		do {
+			PEdge *nexte = p_wheel_edge_next(e);
+
+			points[i][0] = e->next->vert->uv[0]; 
+			points[i][1] = e->next->vert->uv[1]; 
+
+			if (nexte == NULL) {
+				i++;
+				points[i][0] = e->next->next->vert->uv[0]; 
+				points[i][1] = e->next->next->vert->uv[1]; 
+				break;
+			}
+
+			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(PHandle *handle, 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(handle, 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 */
+		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;
+}
+
+static 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;
+}
+
+static 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;
+
+		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];
+	}
+}
+
+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;
+}
+
+static 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);
+
+		// XXX 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];
+				}
+			}
+		}
+	}
+
+	/* 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);
+
+		// XXX 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->aspx = 1.0f;
+	handle->aspy = 1.0f;
+
+	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_aspect_ratio(ParamHandle *handle, float aspx, float aspy)
+{
+	PHandle *phandle = (PHandle*)handle;
+
+	phandle->aspx = aspx;
+	phandle->aspy = aspy;
+}
+
+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;
+
+	for (i = j = 0; i < phandle->ncharts; i++) {
+		PVert *v;
+		chart = phandle->charts[i];
+
+		p_chart_boundaries(chart, &nboundaries, &outer);
+
+		if (!impl && nboundaries == 0) {
+			p_chart_delete(chart);
+			continue;
+		}
+
+		phandle->charts[j] = chart;
+		j++;
+
+		if (fill && (nboundaries > 1))
+			p_chart_fill_boundaries(chart, outer);
+
+		for (v=chart->verts; v; v=v->nextlink)
+			p_vert_load_pin_select_uvs(handle, v);
+	}
+
+	phandle->ncharts = j;
+
+	phandle->state = PHANDLE_STATE_CONSTRUCTED;
+}
+
+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++) {
+		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)
+{
+	PHandle *phandle = (PHandle*)handle;
+	PChart *chart;
+	int i;
+	PBool result;
+
+	param_assert(phandle->state == PHANDLE_STATE_LSCM);
+
+	for (i = 0; i < phandle->ncharts; i++) {
+		chart = phandle->charts[i];
+
+		if (chart->u.lscm.context) {
+			result = p_chart_lscm_solve(phandle, chart);
+
+			if (result && !(chart->flag & PCHART_NOPACK))
+				p_chart_rotate_minimum_area(chart);
+
+			if (!result || (chart->u.lscm.pin1))
+				p_chart_lscm_end(chart);
+		}
+	}
+}
+
+void param_lscm_end(ParamHandle *handle)
+{
+	PHandle *phandle = (PHandle*)handle;
+	int i;
+
+	param_assert(phandle->state == PHANDLE_STATE_LSCM);
+
+	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;
+}
+
+void param_stretch_begin(ParamHandle *handle)
+{
+	PHandle *phandle = (PHandle*)handle;
+	PChart *chart;
+	PVert *v;
+	PFace *f;
+	int i;
+
+	param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
+	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=chart->verts; v; v=v->nextlink)
+			v->flag &= ~PVERT_PIN; /* don't use user-defined pins */
+
+		p_stretch_pin_boundary(chart);
+
+		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;
+	PChart *chart;
+	int i;
+
+	param_assert(phandle->state == PHANDLE_STATE_STRETCH);
+
+	for (i = 0; i < phandle->ncharts; i++) {
+		chart = phandle->charts[i];
+		p_chart_stretch_minimize(chart, phandle->rng);
+	}
+}
+
+void param_stretch_end(ParamHandle *handle)
+{
+	PHandle *phandle = (PHandle*)handle;
+
+	param_assert(phandle->state == PHANDLE_STATE_STRETCH);
+	phandle->state = PHANDLE_STATE_CONSTRUCTED;
+
+	rng_free(phandle->rng);
+	phandle->rng = NULL;
+}
+
+void param_smooth_area(ParamHandle *handle)
+{
+	PHandle *phandle = (PHandle*)handle;
+	int i;
+
+	param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
+
+	for (i = 0; i < phandle->ncharts; i++) {
+		PChart *chart = phandle->charts[i];
+		PVert *v;
+
+		for (v=chart->verts; v; v=v->nextlink)
+			v->flag &= ~PVERT_PIN;
+
+		p_smooth(chart);
+	}
+}
+ 
+void param_pack(ParamHandle *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;
+	
+	if(phandle->aspx != phandle->aspy)
+		param_scale(handle, 1.0f/phandle->aspx, 1.0f/phandle->aspy);
+	
+	/* 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);
+
+	if(phandle->aspx != phandle->aspy)
+		param_scale(handle, phandle->aspx, phandle->aspy);
+}
+
+void param_average(ParamHandle *handle)
+{
+	PChart *chart;
+	int i;
+	float tot_uvarea = 0.0f, tot_facearea = 0.0f;
+	float tot_fac, fac;
+	float minv[2], maxv[2], trans[2];
+	PHandle *phandle = (PHandle*)handle;
+	
+	if (phandle->ncharts == 0)
+		return;
+	
+	for (i = 0; i < phandle->ncharts; i++) {
+		PFace *f;
+		chart = phandle->charts[i];
+		
+		chart->u.pack.area = 0.0f; /* 3d area */
+		chart->u.pack.rescale = 0.0f; /* UV area, abusing rescale for tmp storage, oh well :/ */
+		
+		for (f=chart->faces; f; f=f->nextlink) {
+			chart->u.pack.area += p_face_area(f);
+			chart->u.pack.rescale += fabs(p_face_uv_area_signed(f));
+		}
+		
+		tot_facearea += chart->u.pack.area;
+		tot_uvarea += chart->u.pack.rescale;
+	}
+	
+	if (tot_facearea == tot_uvarea || tot_facearea==0.0f || tot_uvarea==0.0f) {
+		/* nothing to do */
+		return;
+	}
+	
+	tot_fac = tot_facearea/tot_uvarea;
+	
+	for (i = 0; i < phandle->ncharts; i++) {
+		chart = phandle->charts[i];
+		if (chart->u.pack.area != 0.0f && chart->u.pack.rescale != 0.0f) {
+			fac = chart->u.pack.area / chart->u.pack.rescale;
+			
+			/* Get the island center */
+			p_chart_uv_bbox(chart, minv, maxv);
+			trans[0] = (minv[0] + maxv[0]) /-2.0f;
+			trans[1] = (minv[1] + maxv[1]) /-2.0f;
+			
+			/* Move center to 0,0 */
+			p_chart_uv_translate(chart, trans);
+			p_chart_uv_scale(chart, sqrt(fac / tot_fac));
+			
+			/* Move to original center */
+			trans[0] = -trans[0];
+			trans[1] = -trans[1];
+			p_chart_uv_translate(chart, trans);
+		}
+	}
+}
+
+void param_scale(ParamHandle *handle, float x, float y)
+{
+	PHandle *phandle = (PHandle*)handle;
+	PChart *chart;
+	int i;
+
+	for (i = 0; i < phandle->ncharts; i++) {
+		chart = phandle->charts[i];
+		p_chart_uv_scale_xy(chart, x, y);
+	}
+}
+
+void param_flush(ParamHandle *handle)
+{
+	PHandle *phandle = (PHandle*)handle;
+	PChart *chart;
+	int i;
+
+	for (i = 0; i < phandle->ncharts; i++) {
+		chart = phandle->charts[i];
+
+		if ((phandle->state == PHANDLE_STATE_LSCM) && !chart->u.lscm.context)
+			continue;
+
+		if (phandle->blend == 0.0f)
+			p_flush_uvs(phandle, chart);
+		else
+			p_flush_uvs_blend(phandle, chart, phandle->blend);
+	}
+}
+
+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];
+
+		for (f=chart->faces; f; f=f->nextlink)
+			p_face_restore_uvs(f);
+	}
+}
+