Orange branch commit.

This commit adds new underlying uv unwrapper code, intended to be
more extensible. At the moment this has a re-implementation of LSCM.
This has not been activated yet, since it doesn't add anything new.

What's new is the stretch minimize tool from tuhopuu. It works by
selecting some some uv's in the uv editor window, and then pressing
ctrl+V. The uv's on the boundary stay fixed.

More stuff will follow as I port it over & fix it.

9 files changed, 3031 insertions, 832 deletions

diff --git a/intern/opennl/extern/ONL_opennl.h b/intern/opennl/extern/ONL_opennl.h
index 128f1b137bc..345cf0dc717 100644
--- a/intern/opennl/extern/ONL_opennl.h
+++ b/intern/opennl/extern/ONL_opennl.h
@@ -89,10 +89,6 @@ typedef void* NLContext;
 #define NL_SYMMETRIC        0x106
 #define NL_ERROR            0x108
 
-/* Enable / Disable */
-
-#define NL_NORMALIZE_ROWS  0x400
-
 /* Row parameters */
 
 #define NL_RIGHT_HAND_SIDE 0x500
@@ -142,7 +138,9 @@ void nlRightHandSideAdd(NLuint index, NLfloat value);
 
 /* Solve */
 
-NLboolean nlSolve(void);
+void nlPrintMatrix(void);
+NLboolean nlSolve();
+NLboolean nlSolveAdvanced(NLint *permutation, NLboolean solveAgain);
 
 #ifdef __cplusplus
 }
diff --git a/intern/opennl/intern/opennl.c b/intern/opennl/intern/opennl.c
index 7773582e027..ad40f45c73a 100644
--- a/intern/opennl/intern/opennl.c
+++ b/intern/opennl/intern/opennl.c
@@ -24,13 +24,13 @@
  *
  *  Contact: Bruno Levy
  *
- *     levy@loria.fr
+ *	 levy@loria.fr
  *
- *     ISA Project
- *     LORIA, INRIA Lorraine, 
- *     Campus Scientifique, BP 239
- *     54506 VANDOEUVRE LES NANCY CEDEX 
- *     FRANCE
+ *	 ISA Project
+ *	 LORIA, INRIA Lorraine, 
+ *	 Campus Scientifique, BP 239
+ *	 54506 VANDOEUVRE LES NANCY CEDEX 
+ *	 FRANCE
  *
  *  Note that the GNU General Public License does not permit incorporating
  *  the Software into proprietary programs. 
@@ -58,51 +58,51 @@
 
 
 static void __nl_assertion_failed(char* cond, char* file, int line) {
-    fprintf(
-        stderr, 
-        "OpenNL assertion failed: %s, file:%s, line:%d\n",
-        cond,file,line
-    );
-    abort();
+	fprintf(
+		stderr, 
+		"OpenNL assertion failed: %s, file:%s, line:%d\n",
+		cond,file,line
+	);
+	abort();
 }
 
 static void __nl_range_assertion_failed(
-    float x, float min_val, float max_val, char* file, int line
+	float x, float min_val, float max_val, char* file, int line
 ) {
-    fprintf(
-        stderr, 
-        "OpenNL range assertion failed: %f in [ %f ... %f ], file:%s, line:%d\n",
-        x, min_val, max_val, file,line
-    );
-    abort();
+	fprintf(
+		stderr, 
+		"OpenNL range assertion failed: %f in [ %f ... %f ], file:%s, line:%d\n",
+		x, min_val, max_val, file,line
+	);
+	abort();
 }
 
 static void __nl_should_not_have_reached(char* file, int line) {
-    fprintf(
-        stderr, 
-        "OpenNL should not have reached this point: file:%s, line:%d\n",
-        file,line
-    );
-    abort();
+	fprintf(
+		stderr, 
+		"OpenNL should not have reached this point: file:%s, line:%d\n",
+		file,line
+	);
+	abort();
 }
 
 
-#define __nl_assert(x) {                                        \
-    if(!(x)) {                                                  \
-        __nl_assertion_failed(#x,__FILE__, __LINE__);          \
-    }                                                           \
+#define __nl_assert(x) {										\
+	if(!(x)) {												  \
+		__nl_assertion_failed(#x,__FILE__, __LINE__);		  \
+	}														   \
 } 
 
-#define __nl_range_assert(x,min_val,max_val) {                  \
-    if(((x) < (min_val)) || ((x) > (max_val))) {                \
-        __nl_range_assertion_failed(x, min_val, max_val,        \
-            __FILE__, __LINE__                                  \
-        );                                                     \
-    }                                                           \
+#define __nl_range_assert(x,min_val,max_val) {				  \
+	if(((x) < (min_val)) || ((x) > (max_val))) {				\
+		__nl_range_assertion_failed(x, min_val, max_val,		\
+			__FILE__, __LINE__								  \
+		);													 \
+	}														   \
 }
 
-#define __nl_assert_not_reached {                               \
-    __nl_should_not_have_reached(__FILE__, __LINE__);          \
+#define __nl_assert_not_reached {							   \
+	__nl_should_not_have_reached(__FILE__, __LINE__);		  \
 }
 
 #ifdef NL_DEBUG
@@ -135,301 +135,301 @@ static void __nl_should_not_have_reached(char* file, int line) {
 /************************************************************************************/
 /* memory management */
 
-#define __NL_NEW(T)                (T*)(calloc(1, sizeof(T))) 
-#define __NL_NEW_ARRAY(T,NB)       (T*)(calloc((NB),sizeof(T))) 
+#define __NL_NEW(T)				(T*)(calloc(1, sizeof(T))) 
+#define __NL_NEW_ARRAY(T,NB)	   (T*)(calloc((NB),sizeof(T))) 
 #define __NL_RENEW_ARRAY(T,x,NB)   (T*)(realloc(x,(NB)*sizeof(T))) 
-#define __NL_DELETE(x)             free(x); x = NULL 
-#define __NL_DELETE_ARRAY(x)       free(x); x = NULL
+#define __NL_DELETE(x)			 free(x); x = NULL 
+#define __NL_DELETE_ARRAY(x)	   free(x); x = NULL
 
-#define __NL_CLEAR(T, x)           memset(x, 0, sizeof(T)) 
+#define __NL_CLEAR(T, x)		   memset(x, 0, sizeof(T)) 
 #define __NL_CLEAR_ARRAY(T,x,NB)   memset(x, 0, (NB)*sizeof(T)) 
 
 /************************************************************************************/
 /* Dynamic arrays for sparse row/columns */
 
 typedef struct {
-    NLuint   index;
-    NLfloat value;
+	NLuint   index;
+	NLfloat value;
 } __NLCoeff;
 
 typedef struct {
-    NLuint size;
-    NLuint capacity;
-    __NLCoeff* coeff;
+	NLuint size;
+	NLuint capacity;
+	__NLCoeff* coeff;
 } __NLRowColumn;
 
 static void __nlRowColumnConstruct(__NLRowColumn* c) {
-    c->size     = 0;
-    c->capacity = 0;
-    c->coeff    = NULL;
+	c->size	 = 0;
+	c->capacity = 0;
+	c->coeff	= NULL;
 }
 
 static void __nlRowColumnDestroy(__NLRowColumn* c) {
-    __NL_DELETE_ARRAY(c->coeff);
+	__NL_DELETE_ARRAY(c->coeff);
 #ifdef NL_PARANOID
-    __NL_CLEAR(__NLRowColumn, c); 
+	__NL_CLEAR(__NLRowColumn, c); 
 #endif
 }
 
 static void __nlRowColumnGrow(__NLRowColumn* c) {
-    if(c->capacity != 0) {
-        c->capacity = 2 * c->capacity;
-        c->coeff = __NL_RENEW_ARRAY(__NLCoeff, c->coeff, c->capacity);
-    } else {
-        c->capacity = 4;
-        c->coeff = __NL_NEW_ARRAY(__NLCoeff, c->capacity);
-    }
+	if(c->capacity != 0) {
+		c->capacity = 2 * c->capacity;
+		c->coeff = __NL_RENEW_ARRAY(__NLCoeff, c->coeff, c->capacity);
+	} else {
+		c->capacity = 4;
+		c->coeff = __NL_NEW_ARRAY(__NLCoeff, c->capacity);
+	}
 }
 
 static void __nlRowColumnAdd(__NLRowColumn* c, NLint index, NLfloat value) {
-    NLuint i;
-    for(i=0; i<c->size; i++) {
-        if(c->coeff[i].index == (NLuint)index) {
-            c->coeff[i].value += value;
-            return;
-        }
-    }
-    if(c->size == c->capacity) {
-        __nlRowColumnGrow(c);
-    }
-    c->coeff[c->size].index = index;
-    c->coeff[c->size].value = value;
-    c->size++;
+	NLuint i;
+	for(i=0; i<c->size; i++) {
+		if(c->coeff[i].index == (NLuint)index) {
+			c->coeff[i].value += value;
+			return;
+		}
+	}
+	if(c->size == c->capacity) {
+		__nlRowColumnGrow(c);
+	}
+	c->coeff[c->size].index = index;
+	c->coeff[c->size].value = value;
+	c->size++;
 }
 
 /* Does not check whether the index already exists */
 static void __nlRowColumnAppend(__NLRowColumn* c, NLint index, NLfloat value) {
-    if(c->size == c->capacity) {
-        __nlRowColumnGrow(c);
-    }
-    c->coeff[c->size].index = index;
-    c->coeff[c->size].value = value;
-    c->size++;
+	if(c->size == c->capacity) {
+		__nlRowColumnGrow(c);
+	}
+	c->coeff[c->size].index = index;
+	c->coeff[c->size].value = value;
+	c->size++;
 }
 
 static void __nlRowColumnZero(__NLRowColumn* c) {
-    c->size = 0;
+	c->size = 0;
 }
 
 static void __nlRowColumnClear(__NLRowColumn* c) {
-    c->size     = 0;
-    c->capacity = 0;
-    __NL_DELETE_ARRAY(c->coeff);
+	c->size	 = 0;
+	c->capacity = 0;
+	__NL_DELETE_ARRAY(c->coeff);
 }
 
 /************************************************************************************/
 /* SparseMatrix data structure */
 
-#define __NL_ROWS      1
+#define __NL_ROWS	  1
 #define __NL_COLUMNS   2
 #define __NL_SYMMETRIC 4
 
 typedef struct {
-    NLuint m;
-    NLuint n;
-    NLuint diag_size;
-    NLenum storage;
-    __NLRowColumn* row;
-    __NLRowColumn* column;
-    NLfloat*      diag;
+	NLuint m;
+	NLuint n;
+	NLuint diag_size;
+	NLenum storage;
+	__NLRowColumn* row;
+	__NLRowColumn* column;
+	NLfloat*	  diag;
 } __NLSparseMatrix;
 
 
 static void __nlSparseMatrixConstruct(
-    __NLSparseMatrix* M, NLuint m, NLuint n, NLenum storage
+	__NLSparseMatrix* M, NLuint m, NLuint n, NLenum storage
 ) {
-    NLuint i;
-    M->m = m;
-    M->n = n;
-    M->storage = storage;
-    if(storage & __NL_ROWS) {
-        M->row = __NL_NEW_ARRAY(__NLRowColumn, m);
-        for(i=0; i<n; i++) {
-            __nlRowColumnConstruct(&(M->row[i]));
-        }
-    } else {
-        M->row = NULL;
-    }
-
-    if(storage & __NL_COLUMNS) {
-        M->column = __NL_NEW_ARRAY(__NLRowColumn, n);
-        for(i=0; i<n; i++) {
-            __nlRowColumnConstruct(&(M->column[i]));
-        }
-    } else {
-        M->column = NULL;
-    }
-
-    M->diag_size = MIN(m,n);
-    M->diag = __NL_NEW_ARRAY(NLfloat, M->diag_size);
+	NLuint i;
+	M->m = m;
+	M->n = n;
+	M->storage = storage;
+	if(storage & __NL_ROWS) {
+		M->row = __NL_NEW_ARRAY(__NLRowColumn, m);
+		for(i=0; i<n; i++) {
+			__nlRowColumnConstruct(&(M->row[i]));
+		}
+	} else {
+		M->row = NULL;
+	}
+
+	if(storage & __NL_COLUMNS) {
+		M->column = __NL_NEW_ARRAY(__NLRowColumn, n);
+		for(i=0; i<n; i++) {
+			__nlRowColumnConstruct(&(M->column[i]));
+		}
+	} else {
+		M->column = NULL;
+	}
+
+	M->diag_size = MIN(m,n);
+	M->diag = __NL_NEW_ARRAY(NLfloat, M->diag_size);
 }
 
 static void __nlSparseMatrixDestroy(__NLSparseMatrix* M) {
-    NLuint i;
-    __NL_DELETE_ARRAY(M->diag);
-    if(M->storage & __NL_ROWS) {
-        for(i=0; i<M->m; i++) {
-            __nlRowColumnDestroy(&(M->row[i]));
-        }
-        __NL_DELETE_ARRAY(M->row);
-    }
-    if(M->storage & __NL_COLUMNS) {
-        for(i=0; i<M->n; i++) {
-            __nlRowColumnDestroy(&(M->column[i]));
-        }
-        __NL_DELETE_ARRAY(M->column);
-    }
+	NLuint i;
+	__NL_DELETE_ARRAY(M->diag);
+	if(M->storage & __NL_ROWS) {
+		for(i=0; i<M->m; i++) {
+			__nlRowColumnDestroy(&(M->row[i]));
+		}
+		__NL_DELETE_ARRAY(M->row);
+	}
+	if(M->storage & __NL_COLUMNS) {
+		for(i=0; i<M->n; i++) {
+			__nlRowColumnDestroy(&(M->column[i]));
+		}
+		__NL_DELETE_ARRAY(M->column);
+	}
 #ifdef NL_PARANOID
-    __NL_CLEAR(__NLSparseMatrix,M);
+	__NL_CLEAR(__NLSparseMatrix,M);
 #endif
 }
 
 static void __nlSparseMatrixAdd(
-    __NLSparseMatrix* M, NLuint i, NLuint j, NLfloat value
+	__NLSparseMatrix* M, NLuint i, NLuint j, NLfloat value
 ) {
-    __nl_parano_range_assert(i, 0, M->m - 1);
-    __nl_parano_range_assert(j, 0, M->n - 1);
-    if((M->storage & __NL_SYMMETRIC) && (j > i)) {
-        return;
-    }
-    if(i == j) {
-        M->diag[i] += value;
-    }
-    if(M->storage & __NL_ROWS) {
-        __nlRowColumnAdd(&(M->row[i]), j, value);
-    }
-    if(M->storage & __NL_COLUMNS) {
-        __nlRowColumnAdd(&(M->column[j]), i, value);
-    }
+	__nl_parano_range_assert(i, 0, M->m - 1);
+	__nl_parano_range_assert(j, 0, M->n - 1);
+	if((M->storage & __NL_SYMMETRIC) && (j > i)) {
+		return;
+	}
+	if(i == j) {
+		M->diag[i] += value;
+	}
+	if(M->storage & __NL_ROWS) {
+		__nlRowColumnAdd(&(M->row[i]), j, value);
+	}
+	if(M->storage & __NL_COLUMNS) {
+		__nlRowColumnAdd(&(M->column[j]), i, value);
+	}
 }
 
 static void __nlSparseMatrixClear( __NLSparseMatrix* M) {
-    NLuint i;
-    if(M->storage & __NL_ROWS) {
-        for(i=0; i<M->m; i++) {
-            __nlRowColumnClear(&(M->row[i]));
-        }
-    }
-    if(M->storage & __NL_COLUMNS) {
-        for(i=0; i<M->n; i++) {
-            __nlRowColumnClear(&(M->column[i]));
-        }
-    }
-    __NL_CLEAR_ARRAY(NLfloat, M->diag, M->diag_size);    
+	NLuint i;
+	if(M->storage & __NL_ROWS) {
+		for(i=0; i<M->m; i++) {
+			__nlRowColumnClear(&(M->row[i]));
+		}
+	}
+	if(M->storage & __NL_COLUMNS) {
+		for(i=0; i<M->n; i++) {
+			__nlRowColumnClear(&(M->column[i]));
+		}
+	}
+	__NL_CLEAR_ARRAY(NLfloat, M->diag, M->diag_size);	
 }
 
 /* Returns the number of non-zero coefficients */
 static NLuint __nlSparseMatrixNNZ( __NLSparseMatrix* M) {
-    NLuint nnz = 0;
-    NLuint i;
-    if(M->storage & __NL_ROWS) {
-        for(i = 0; i<M->m; i++) {
-            nnz += M->row[i].size;
-        }
-    } else if (M->storage & __NL_COLUMNS) {
-        for(i = 0; i<M->n; i++) {
-            nnz += M->column[i].size;
-        }
-    } else {
-        __nl_assert_not_reached;
-    }
-    return nnz;
+	NLuint nnz = 0;
+	NLuint i;
+	if(M->storage & __NL_ROWS) {
+		for(i = 0; i<M->m; i++) {
+			nnz += M->row[i].size;
+		}
+	} else if (M->storage & __NL_COLUMNS) {
+		for(i = 0; i<M->n; i++) {
+			nnz += M->column[i].size;
+		}
+	} else {
+		__nl_assert_not_reached;
+	}
+	return nnz;
 }
 
 /************************************************************************************/
 /* SparseMatrix x Vector routines, internal helper routines */
 
 static void __nlSparseMatrix_mult_rows_symmetric(
-    __NLSparseMatrix* A, NLfloat* x, NLfloat* y
+	__NLSparseMatrix* A, NLfloat* x, NLfloat* y
 ) {
-    NLuint m = A->m;
-    NLuint i,ij;
-    __NLRowColumn* Ri = NULL;
-    __NLCoeff* c = NULL;
-    for(i=0; i<m; i++) {
-        y[i] = 0;
-        Ri = &(A->row[i]);
-        for(ij=0; ij<Ri->size; ij++) {
-            c = &(Ri->coeff[ij]);
-            y[i] += c->value * x[c->index];
-            if(i != c->index) {
-                y[c->index] += c->value * x[i];
-            }
-        }
-    }
+	NLuint m = A->m;
+	NLuint i,ij;
+	__NLRowColumn* Ri = NULL;
+	__NLCoeff* c = NULL;
+	for(i=0; i<m; i++) {
+		y[i] = 0;
+		Ri = &(A->row[i]);
+		for(ij=0; ij<Ri->size; ij++) {
+			c = &(Ri->coeff[ij]);
+			y[i] += c->value * x[c->index];
+			if(i != c->index) {
+				y[c->index] += c->value * x[i];
+			}
+		}
+	}
 }
 
 static void __nlSparseMatrix_mult_rows(
-    __NLSparseMatrix* A, NLfloat* x, NLfloat* y
+	__NLSparseMatrix* A, NLfloat* x, NLfloat* y
 ) {
-    NLuint m = A->m;
-    NLuint i,ij;
-    __NLRowColumn* Ri = NULL;
-    __NLCoeff* c = NULL;
-    for(i=0; i<m; i++) {
-        y[i] = 0;
-        Ri = &(A->row[i]);
-        for(ij=0; ij<Ri->size; ij++) {
-            c = &(Ri->coeff[ij]);
-            y[i] += c->value * x[c->index];
-        }
-    }
+	NLuint m = A->m;
+	NLuint i,ij;
+	__NLRowColumn* Ri = NULL;
+	__NLCoeff* c = NULL;
+	for(i=0; i<m; i++) {
+		y[i] = 0;
+		Ri = &(A->row[i]);
+		for(ij=0; ij<Ri->size; ij++) {
+			c = &(Ri->coeff[ij]);
+			y[i] += c->value * x[c->index];
+		}
+	}
 }
 
 static void __nlSparseMatrix_mult_cols_symmetric(
-    __NLSparseMatrix* A, NLfloat* x, NLfloat* y
+	__NLSparseMatrix* A, NLfloat* x, NLfloat* y
 ) {
-    NLuint n = A->n;
-    NLuint j,ii;
-    __NLRowColumn* Cj = NULL;
-    __NLCoeff* c = NULL;
-    for(j=0; j<n; j++) {
-        y[j] = 0;
-        Cj = &(A->column[j]);
-        for(ii=0; ii<Cj->size; ii++) {
-            c = &(Cj->coeff[ii]);
-            y[c->index] += c->value * x[j];
-            if(j != c->index) {
-                y[j] += c->value * x[c->index];
-            }
-        }
-    }
+	NLuint n = A->n;
+	NLuint j,ii;
+	__NLRowColumn* Cj = NULL;
+	__NLCoeff* c = NULL;
+	for(j=0; j<n; j++) {
+		y[j] = 0;
+		Cj = &(A->column[j]);
+		for(ii=0; ii<Cj->size; ii++) {
+			c = &(Cj->coeff[ii]);
+			y[c->index] += c->value * x[j];
+			if(j != c->index) {
+				y[j] += c->value * x[c->index];
+			}
+		}
+	}
 }
 
 static void __nlSparseMatrix_mult_cols(
-    __NLSparseMatrix* A, NLfloat* x, NLfloat* y
+	__NLSparseMatrix* A, NLfloat* x, NLfloat* y
 ) {
-    NLuint n = A->n;
-    NLuint j,ii; 
-    __NLRowColumn* Cj = NULL;
-    __NLCoeff* c = NULL;
-    __NL_CLEAR_ARRAY(NLfloat, y, A->m);
-    for(j=0; j<n; j++) {
-        Cj = &(A->column[j]);
-        for(ii=0; ii<Cj->size; ii++) {
-            c = &(Cj->coeff[ii]);
-            y[c->index] += c->value * x[j];
-        }
-    }
+	NLuint n = A->n;
+	NLuint j,ii; 
+	__NLRowColumn* Cj = NULL;
+	__NLCoeff* c = NULL;
+	__NL_CLEAR_ARRAY(NLfloat, y, A->m);
+	for(j=0; j<n; j++) {
+		Cj = &(A->column[j]);
+		for(ii=0; ii<Cj->size; ii++) {
+			c = &(Cj->coeff[ii]);
+			y[c->index] += c->value * x[j];
+		}
+	}
 }
 
 /************************************************************************************/
 /* SparseMatrix x Vector routines, main driver routine */
 
 static void __nlSparseMatrixMult(__NLSparseMatrix* A, NLfloat* x, NLfloat* y) {
-    if(A->storage & __NL_ROWS) {
-        if(A->storage & __NL_SYMMETRIC) {
-            __nlSparseMatrix_mult_rows_symmetric(A, x, y);
-        } else {
-            __nlSparseMatrix_mult_rows(A, x, y);
-        }
-    } else {
-        if(A->storage & __NL_SYMMETRIC) {
-            __nlSparseMatrix_mult_cols_symmetric(A, x, y);
-        } else {
-            __nlSparseMatrix_mult_cols(A, x, y);
-        }
-    }
+	if(A->storage & __NL_ROWS) {
+		if(A->storage & __NL_SYMMETRIC) {
+			__nlSparseMatrix_mult_rows_symmetric(A, x, y);
+		} else {
+			__nlSparseMatrix_mult_rows(A, x, y);
+		}
+	} else {
+		if(A->storage & __NL_SYMMETRIC) {
+			__nlSparseMatrix_mult_cols_symmetric(A, x, y);
+		} else {
+			__nlSparseMatrix_mult_cols(A, x, y);
+		}
+	}
 }
 
 /************************************************************************************/
@@ -438,513 +438,468 @@ static void __nlSparseMatrixMult(__NLSparseMatrix* A, NLfloat* x, NLfloat* y) {
 typedef void(*__NLMatrixFunc)(float* x, float* y);
 
 typedef struct {
-    NLfloat  value;
-    NLboolean locked;
-    NLuint    index;
+	NLfloat  value;
+	NLboolean locked;
+	NLuint	index;
 } __NLVariable;
 
-#define __NL_STATE_INITIAL            0
-#define __NL_STATE_SYSTEM             1
-#define __NL_STATE_MATRIX             2
-#define __NL_STATE_ROW                3
-#define __NL_STATE_MATRIX_CONSTRUCTED 4
-#define __NL_STATE_SYSTEM_CONSTRUCTED 5
-#define __NL_STATE_SOLVED             6
+#define __NL_STATE_INITIAL				0
+#define __NL_STATE_SYSTEM				1
+#define __NL_STATE_MATRIX				2
+#define __NL_STATE_ROW					3
+#define __NL_STATE_MATRIX_CONSTRUCTED	4
+#define __NL_STATE_SYSTEM_CONSTRUCTED	5
+#define __NL_STATE_SYSTEM_SOLVED		7
 
 typedef struct {
-    NLenum           state;
-    __NLVariable*    variable;
-    NLuint           n;
-    __NLSparseMatrix M;
-    __NLRowColumn    af;
-    __NLRowColumn    al;
-    __NLRowColumn    xl;
-    NLfloat*        x;
-    NLfloat*        b;
-    NLfloat         right_hand_side;
-    NLfloat         row_scaling;
-    NLuint           nb_variables;
-    NLuint           current_row;
-    NLboolean        least_squares;
-    NLboolean        symmetric;
-    NLboolean        normalize_rows;
-    NLboolean        alloc_M;
-    NLboolean        alloc_af;
-    NLboolean        alloc_al;
-    NLboolean        alloc_xl;
-    NLboolean        alloc_variable;
-    NLboolean        alloc_x;
-    NLboolean        alloc_b;
-    NLfloat         error;
-    __NLMatrixFunc   matrix_vector_prod;
+	NLenum		   state;
+	__NLVariable*	variable;
+	NLuint		   n;
+	__NLSparseMatrix M;
+	__NLRowColumn	af;
+	__NLRowColumn	al;
+	NLfloat*		x;
+	NLfloat*		b;
+	NLfloat		 right_hand_side;
+	NLuint		   nb_variables;
+	NLuint		   current_row;
+	NLboolean		least_squares;
+	NLboolean		symmetric;
+	NLboolean		solve_again;
+	NLboolean		alloc_M;
+	NLboolean		alloc_af;
+	NLboolean		alloc_al;
+	NLboolean		alloc_variable;
+	NLboolean		alloc_x;
+	NLboolean		alloc_b;
+	NLfloat		 error;
+	__NLMatrixFunc   matrix_vector_prod;
+
+	struct __NLSuperLUContext {
+		NLboolean alloc_slu;
+		SuperMatrix L, U;
+		NLint *perm_c, *perm_r;
+		SuperLUStat_t stat;
+	} slu;
 } __NLContext;
 
 static __NLContext* __nlCurrentContext = NULL;
 
 static void __nlMatrixVectorProd_default(NLfloat* x, NLfloat* y) {
-    __nlSparseMatrixMult(&(__nlCurrentContext->M), x, y);
+	__nlSparseMatrixMult(&(__nlCurrentContext->M), x, y);
 }
 
 
 NLContext nlNewContext(void) {
-    __NLContext* result      = __NL_NEW(__NLContext);
-    result->state            = __NL_STATE_INITIAL;
-    result->row_scaling      = 1.0;
-    result->right_hand_side  = 0.0;
-    result->matrix_vector_prod = __nlMatrixVectorProd_default;
-    nlMakeCurrent(result);
-    return result;
+	__NLContext* result	  = __NL_NEW(__NLContext);
+	result->state			= __NL_STATE_INITIAL;
+	result->right_hand_side  = 0.0;
+	result->matrix_vector_prod = __nlMatrixVectorProd_default;
+	nlMakeCurrent(result);
+	return result;
 }
 
+static void __nlFree_SUPERLU(__NLContext *context);
+
 void nlDeleteContext(NLContext context_in) {
-    __NLContext* context = (__NLContext*)(context_in);
-    if(__nlCurrentContext == context) {
-        __nlCurrentContext = NULL;
-    }
-    if(context->alloc_M) {
-        __nlSparseMatrixDestroy(&context->M);
-    }
-    if(context->alloc_af) {
-        __nlRowColumnDestroy(&context->af);
-    }
-    if(context->alloc_al) {
-        __nlRowColumnDestroy(&context->al);
-    }
-    if(context->alloc_xl) {
-        __nlRowColumnDestroy(&context->xl);
-    }
-    if(context->alloc_variable) {
-        __NL_DELETE_ARRAY(context->variable);
-    }
-    if(context->alloc_x) {
-        __NL_DELETE_ARRAY(context->x);
-    }
-    if(context->alloc_b) {
-        __NL_DELETE_ARRAY(context->b);
-    }
+	__NLContext* context = (__NLContext*)(context_in);
+	if(__nlCurrentContext == context) {
+		__nlCurrentContext = NULL;
+	}
+	if(context->alloc_M) {
+		__nlSparseMatrixDestroy(&context->M);
+	}
+	if(context->alloc_af) {
+		__nlRowColumnDestroy(&context->af);
+	}
+	if(context->alloc_al) {
+		__nlRowColumnDestroy(&context->al);
+	}
+	if(context->alloc_variable) {
+		__NL_DELETE_ARRAY(context->variable);
+	}
+	if(context->alloc_x) {
+		__NL_DELETE_ARRAY(context->x);
+	}
+	if(context->alloc_b) {
+		__NL_DELETE_ARRAY(context->b);
+	}
+	if (context->slu.alloc_slu) {
+		__nlFree_SUPERLU(context);
+	}
 
 #ifdef NL_PARANOID
-    __NL_CLEAR(__NLContext, context);
+	__NL_CLEAR(__NLContext, context);
 #endif
-    __NL_DELETE(context);
+	__NL_DELETE(context);
 }
 
 void nlMakeCurrent(NLContext context) {
-    __nlCurrentContext = (__NLContext*)(context);
+	__nlCurrentContext = (__NLContext*)(context);
 }
 
 NLContext nlGetCurrent(void) {
-    return __nlCurrentContext;
+	return __nlCurrentContext;
 }
 
 static void __nlCheckState(NLenum state) {
-    __nl_assert(__nlCurrentContext->state == state);
+	__nl_assert(__nlCurrentContext->state == state);
 }
 
 static void __nlTransition(NLenum from_state, NLenum to_state) {
-    __nlCheckState(from_state);
-    __nlCurrentContext->state = to_state;
+	__nlCheckState(from_state);
+	__nlCurrentContext->state = to_state;
 }
 
 /************************************************************************************/
 /* Get/Set parameters */
 
 void nlSolverParameterf(NLenum pname, NLfloat param) {
-    __nlCheckState(__NL_STATE_INITIAL);
-    switch(pname) {
-    case NL_NB_VARIABLES: {
-        __nl_assert(param > 0);
-        __nlCurrentContext->nb_variables = (NLuint)param;
-    } break;
-    case NL_LEAST_SQUARES: {
-        __nlCurrentContext->least_squares = (NLboolean)param;
-    } break;
-    case NL_SYMMETRIC: {
-        __nlCurrentContext->symmetric = (NLboolean)param;        
-    }
-    default: {
-        __nl_assert_not_reached;
-    } break;
-    }
+	__nlCheckState(__NL_STATE_INITIAL);
+	switch(pname) {
+	case NL_NB_VARIABLES: {
+		__nl_assert(param > 0);
+		__nlCurrentContext->nb_variables = (NLuint)param;
+	} break;
+	case NL_LEAST_SQUARES: {
+		__nlCurrentContext->least_squares = (NLboolean)param;
+	} break;
+	case NL_SYMMETRIC: {
+		__nlCurrentContext->symmetric = (NLboolean)param;		
+	}
+	default: {
+		__nl_assert_not_reached;
+	} break;
+	}
 }
 
 void nlSolverParameteri(NLenum pname, NLint param) {
-    __nlCheckState(__NL_STATE_INITIAL);
-    switch(pname) {
-    case NL_NB_VARIABLES: {
-        __nl_assert(param > 0);
-        __nlCurrentContext->nb_variables = (NLuint)param;
-    } break;
-    case NL_LEAST_SQUARES: {
-        __nlCurrentContext->least_squares = (NLboolean)param;
-    } break;
-    case NL_SYMMETRIC: {
-        __nlCurrentContext->symmetric = (NLboolean)param;        
-    }
-    default: {
-        __nl_assert_not_reached;
-    } break;
-    }
+	__nlCheckState(__NL_STATE_INITIAL);
+	switch(pname) {
+	case NL_NB_VARIABLES: {
+		__nl_assert(param > 0);
+		__nlCurrentContext->nb_variables = (NLuint)param;
+	} break;
+	case NL_LEAST_SQUARES: {
+		__nlCurrentContext->least_squares = (NLboolean)param;
+	} break;
+	case NL_SYMMETRIC: {
+		__nlCurrentContext->symmetric = (NLboolean)param;		
+	}
+	default: {
+		__nl_assert_not_reached;
+	} break;
+	}
 }
 
 void nlRowParameterf(NLenum pname, NLfloat param) {
-    __nlCheckState(__NL_STATE_MATRIX);
-    switch(pname) {
-    case NL_RIGHT_HAND_SIDE: {
-        __nlCurrentContext->right_hand_side = param;
-    } break;
-    case NL_ROW_SCALING: {
-        __nlCurrentContext->row_scaling = param;
-    } break;
-    }
+	__nlCheckState(__NL_STATE_MATRIX);
+	switch(pname) {
+	case NL_RIGHT_HAND_SIDE: {
+		__nlCurrentContext->right_hand_side = param;
+	} break;
+	}
 }
 
 void nlRowParameteri(NLenum pname, NLint param) {
-    __nlCheckState(__NL_STATE_MATRIX);
-    switch(pname) {
-    case NL_RIGHT_HAND_SIDE: {
-        __nlCurrentContext->right_hand_side = (NLfloat)param;
-    } break;
-    case NL_ROW_SCALING: {
-        __nlCurrentContext->row_scaling = (NLfloat)param;
-    } break;
-    }
+	__nlCheckState(__NL_STATE_MATRIX);
+	switch(pname) {
+	case NL_RIGHT_HAND_SIDE: {
+		__nlCurrentContext->right_hand_side = (NLfloat)param;
+	} break;
+	}
 }
 
 void nlGetBooleanv(NLenum pname, NLboolean* params) {
-    switch(pname) {
-    case NL_LEAST_SQUARES: {
-        *params = __nlCurrentContext->least_squares;
-    } break;
-    case NL_SYMMETRIC: {
-        *params = __nlCurrentContext->symmetric;
-    } break;
-    default: {
-        __nl_assert_not_reached;
-    } break;
-    }
+	switch(pname) {
+	case NL_LEAST_SQUARES: {
+		*params = __nlCurrentContext->least_squares;
+	} break;
+	case NL_SYMMETRIC: {
+		*params = __nlCurrentContext->symmetric;
+	} break;
+	default: {
+		__nl_assert_not_reached;
+	} break;
+	}
 }
 
 void nlGetFloatv(NLenum pname, NLfloat* params) {
-    switch(pname) {
-    case NL_NB_VARIABLES: {
-        *params = (NLfloat)(__nlCurrentContext->nb_variables);
-    } break;
-    case NL_LEAST_SQUARES: {
-        *params = (NLfloat)(__nlCurrentContext->least_squares);
-    } break;
-    case NL_SYMMETRIC: {
-        *params = (NLfloat)(__nlCurrentContext->symmetric);
-    } break;
-    case NL_ERROR: {
-        *params = (NLfloat)(__nlCurrentContext->error);
-    } break;
-    default: {
-        __nl_assert_not_reached;
-    } break;
-    }
+	switch(pname) {
+	case NL_NB_VARIABLES: {
+		*params = (NLfloat)(__nlCurrentContext->nb_variables);
+	} break;
+	case NL_LEAST_SQUARES: {
+		*params = (NLfloat)(__nlCurrentContext->least_squares);
+	} break;
+	case NL_SYMMETRIC: {
+		*params = (NLfloat)(__nlCurrentContext->symmetric);
+	} break;
+	case NL_ERROR: {
+		*params = (NLfloat)(__nlCurrentContext->error);
+	} break;
+	default: {
+		__nl_assert_not_reached;
+	} break;
+	}
 }
 
 void nlGetIntergerv(NLenum pname, NLint* params) {
-    switch(pname) {
-    case NL_NB_VARIABLES: {
-        *params = (NLint)(__nlCurrentContext->nb_variables);
-    } break;
-    case NL_LEAST_SQUARES: {
-        *params = (NLint)(__nlCurrentContext->least_squares);
-    } break;
-    case NL_SYMMETRIC: {
-        *params = (NLint)(__nlCurrentContext->symmetric);
-    } break;
-    default: {
-        __nl_assert_not_reached;
-    } break;
-    }
+	switch(pname) {
+	case NL_NB_VARIABLES: {
+		*params = (NLint)(__nlCurrentContext->nb_variables);
+	} break;
+	case NL_LEAST_SQUARES: {
+		*params = (NLint)(__nlCurrentContext->least_squares);
+	} break;
+	case NL_SYMMETRIC: {
+		*params = (NLint)(__nlCurrentContext->symmetric);
+	} break;
+	default: {
+		__nl_assert_not_reached;
+	} break;
+	}
 }
 
 /************************************************************************************/
 /* Enable / Disable */
 
 void nlEnable(NLenum pname) {
-    switch(pname) {
-    case NL_NORMALIZE_ROWS: {
-        __nl_assert(__nlCurrentContext->state != __NL_STATE_ROW);
-        __nlCurrentContext->normalize_rows = NL_TRUE;
-    } break;
-    default: {
-        __nl_assert_not_reached;
-    }
-    }
+	switch(pname) {
+	default: {
+		__nl_assert_not_reached;
+	}
+	}
 }
 
 void nlDisable(NLenum pname) {
-    switch(pname) {
-    case NL_NORMALIZE_ROWS: {
-        __nl_assert(__nlCurrentContext->state != __NL_STATE_ROW);
-        __nlCurrentContext->normalize_rows = NL_FALSE;
-    } break;
-    default: {
-        __nl_assert_not_reached;
-    }
-    }
+	switch(pname) {
+	default: {
+		__nl_assert_not_reached;
+	}
+	}
 }
 
 NLboolean nlIsEnabled(NLenum pname) {
-    switch(pname) {
-    case NL_NORMALIZE_ROWS: {
-        return __nlCurrentContext->normalize_rows;
-    } break;
-    default: {
-        __nl_assert_not_reached;
-    }
-    }
-    return NL_FALSE;
+	switch(pname) {
+	default: {
+		__nl_assert_not_reached;
+	}
+	}
+	return NL_FALSE;
 }
 
 /************************************************************************************/
 /* Get/Set Lock/Unlock variables */
 
 void nlSetVariable(NLuint index, NLfloat value) {
-    __nlCheckState(__NL_STATE_SYSTEM);
-    __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
-    __nlCurrentContext->variable[index].value = value;    
+	__nlCheckState(__NL_STATE_SYSTEM);
+	__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
+	__nlCurrentContext->variable[index].value = value;	
 }
 
 NLfloat nlGetVariable(NLuint index) {
-    __nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL);
-    __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
-    return __nlCurrentContext->variable[index].value;
+	__nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL);
+	__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
+	return __nlCurrentContext->variable[index].value;
 }
 
 void nlLockVariable(NLuint index) {
-    __nlCheckState(__NL_STATE_SYSTEM);
-    __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
-    __nlCurrentContext->variable[index].locked = NL_TRUE;
+	__nlCheckState(__NL_STATE_SYSTEM);
+	__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
+	__nlCurrentContext->variable[index].locked = NL_TRUE;
 }
 
 void nlUnlockVariable(NLuint index) {
-    __nlCheckState(__NL_STATE_SYSTEM);
-    __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
-    __nlCurrentContext->variable[index].locked = NL_FALSE;
+	__nlCheckState(__NL_STATE_SYSTEM);
+	__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
+	__nlCurrentContext->variable[index].locked = NL_FALSE;
 }
 
 NLboolean nlVariableIsLocked(NLuint index) {
-    __nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL);
-    __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
-    return __nlCurrentContext->variable[index].locked;
+	__nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL);
+	__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
+	return __nlCurrentContext->variable[index].locked;
 }
 
 /************************************************************************************/
 /* System construction */
 
 static void __nlVariablesToVector() {
-    NLuint i;
-    __nl_assert(__nlCurrentContext->alloc_x);
-    __nl_assert(__nlCurrentContext->alloc_variable);
-    for(i=0; i<__nlCurrentContext->nb_variables; i++) {
-        __NLVariable* v = &(__nlCurrentContext->variable[i]);
-        if(!v->locked) {
-            __nl_assert(v->index < __nlCurrentContext->n);
-            __nlCurrentContext->x[v->index] = v->value;
-        }
-    }
+	NLuint i;
+	__nl_assert(__nlCurrentContext->alloc_x);
+	__nl_assert(__nlCurrentContext->alloc_variable);
+	for(i=0; i<__nlCurrentContext->nb_variables; i++) {
+		__NLVariable* v = &(__nlCurrentContext->variable[i]);
+		if(!v->locked) {
+			__nl_assert(v->index < __nlCurrentContext->n);
+			__nlCurrentContext->x[v->index] = v->value;
+		}
+	}
 }
 
 static void __nlVectorToVariables() {
-    NLuint i;
-    __nl_assert(__nlCurrentContext->alloc_x);
-    __nl_assert(__nlCurrentContext->alloc_variable);
-    for(i=0; i<__nlCurrentContext->nb_variables; i++) {
-        __NLVariable* v = &(__nlCurrentContext->variable[i]);
-        if(!v->locked) {
-            __nl_assert(v->index < __nlCurrentContext->n);
-            v->value = __nlCurrentContext->x[v->index];
-        }
-    }
+	NLuint i;
+	__nl_assert(__nlCurrentContext->alloc_x);
+	__nl_assert(__nlCurrentContext->alloc_variable);
+	for(i=0; i<__nlCurrentContext->nb_variables; i++) {
+		__NLVariable* v = &(__nlCurrentContext->variable[i]);
+		if(!v->locked) {
+			__nl_assert(v->index < __nlCurrentContext->n);
+			v->value = __nlCurrentContext->x[v->index];
+		}
+	}
 }
 
-
 static void __nlBeginSystem() {
-    __nlTransition(__NL_STATE_INITIAL, __NL_STATE_SYSTEM);
-    __nl_assert(__nlCurrentContext->nb_variables > 0);
-    __nlCurrentContext->variable = __NL_NEW_ARRAY(
-        __NLVariable, __nlCurrentContext->nb_variables
-    );
-    __nlCurrentContext->alloc_variable = NL_TRUE;
+	__nl_assert(__nlCurrentContext->nb_variables > 0);
+
+	if (__nlCurrentContext->solve_again)
+		__nlTransition(__NL_STATE_SYSTEM_SOLVED, __NL_STATE_SYSTEM);
+	else {
+		__nlTransition(__NL_STATE_INITIAL, __NL_STATE_SYSTEM);
+
+		__nlCurrentContext->variable = __NL_NEW_ARRAY(
+			__NLVariable, __nlCurrentContext->nb_variables
+		);
+		__nlCurrentContext->alloc_variable = NL_TRUE;
+	}
 }
 
 static void __nlEndSystem() {
-    __nlTransition(__NL_STATE_MATRIX_CONSTRUCTED, __NL_STATE_SYSTEM_CONSTRUCTED);    
+	__nlTransition(__NL_STATE_MATRIX_CONSTRUCTED, __NL_STATE_SYSTEM_CONSTRUCTED);	
 }
 
 static void __nlBeginMatrix() {
-    NLuint i;
-    NLuint n = 0;
-    NLenum storage = __NL_ROWS;
-
-    __nlTransition(__NL_STATE_SYSTEM, __NL_STATE_MATRIX);
-
-    for(i=0; i<__nlCurrentContext->nb_variables; i++) {
-        if(!__nlCurrentContext->variable[i].locked) {
-            __nlCurrentContext->variable[i].index = n;
-            n++;
-        } else {
-            __nlCurrentContext->variable[i].index = ~0;
-        }
-    }
-
-    __nlCurrentContext->n = n;
-
-    /* a least squares problem results in a symmetric matrix */
-    if(__nlCurrentContext->least_squares) {
-        __nlCurrentContext->symmetric = NL_TRUE;
-    }
-
-    if(__nlCurrentContext->symmetric) {
-        storage = (storage | __NL_SYMMETRIC);
-    }
-
-    /* SuperLU storage does not support symmetric storage */
-    storage = (storage & ~__NL_SYMMETRIC);
-
-    __nlSparseMatrixConstruct(&__nlCurrentContext->M, n, n, storage);
-    __nlCurrentContext->alloc_M = NL_TRUE;
-
-    __nlCurrentContext->x = __NL_NEW_ARRAY(NLfloat, n);
-    __nlCurrentContext->alloc_x = NL_TRUE;
-    
-    __nlCurrentContext->b = __NL_NEW_ARRAY(NLfloat, n);
-    __nlCurrentContext->alloc_b = NL_TRUE;
-
-    __nlVariablesToVector();
-
-    __nlRowColumnConstruct(&__nlCurrentContext->af);
-    __nlCurrentContext->alloc_af = NL_TRUE;
-    __nlRowColumnConstruct(&__nlCurrentContext->al);
-    __nlCurrentContext->alloc_al = NL_TRUE;
-    __nlRowColumnConstruct(&__nlCurrentContext->xl);
-    __nlCurrentContext->alloc_xl = NL_TRUE;
-
-    __nlCurrentContext->current_row = 0;
+	NLuint i;
+	NLuint n = 0;
+	NLenum storage = __NL_ROWS;
+
+	__nlTransition(__NL_STATE_SYSTEM, __NL_STATE_MATRIX);
+
+	if (!__nlCurrentContext->solve_again) {
+		for(i=0; i<__nlCurrentContext->nb_variables; i++) {
+			if(!__nlCurrentContext->variable[i].locked)
+				__nlCurrentContext->variable[i].index = n++;
+			else
+				__nlCurrentContext->variable[i].index = ~0;
+		}
+
+		__nlCurrentContext->n = n;
+
+		/* a least squares problem results in a symmetric matrix */
+		if(__nlCurrentContext->least_squares)
+			__nlCurrentContext->symmetric = NL_TRUE;
+
+		if(__nlCurrentContext->symmetric)
+			storage = (storage | __NL_SYMMETRIC);
+
+		/* SuperLU storage does not support symmetric storage */
+		storage = (storage & ~__NL_SYMMETRIC);
+
+		__nlSparseMatrixConstruct(&__nlCurrentContext->M, n, n, storage);
+		__nlCurrentContext->alloc_M = NL_TRUE;
+
+		__nlCurrentContext->x = __NL_NEW_ARRAY(NLfloat, n);
+		__nlCurrentContext->alloc_x = NL_TRUE;
+		
+		__nlCurrentContext->b = __NL_NEW_ARRAY(NLfloat, n);
+		__nlCurrentContext->alloc_b = NL_TRUE;
+	}
+	else {
+		/* need to recompute b only, A is not constructed anymore */
+		__NL_CLEAR_ARRAY(NLfloat, __nlCurrentContext->b, n);
+	}
+
+	__nlVariablesToVector();
+
+	__nlRowColumnConstruct(&__nlCurrentContext->af);
+	__nlCurrentContext->alloc_af = NL_TRUE;
+	__nlRowColumnConstruct(&__nlCurrentContext->al);
+	__nlCurrentContext->alloc_al = NL_TRUE;
+
+	__nlCurrentContext->current_row = 0;
 }
 
 static void __nlEndMatrix() {
-    __nlTransition(__NL_STATE_MATRIX, __NL_STATE_MATRIX_CONSTRUCTED);    
-    
-    __nlRowColumnDestroy(&__nlCurrentContext->af);
-    __nlCurrentContext->alloc_af = NL_FALSE;
-    __nlRowColumnDestroy(&__nlCurrentContext->al);
-    __nlCurrentContext->alloc_al = NL_FALSE;
-    __nlRowColumnDestroy(&__nlCurrentContext->xl);
-    __nlCurrentContext->alloc_al = NL_FALSE;
-    
+	__nlTransition(__NL_STATE_MATRIX, __NL_STATE_MATRIX_CONSTRUCTED);	
+	
+	__nlRowColumnDestroy(&__nlCurrentContext->af);
+	__nlCurrentContext->alloc_af = NL_FALSE;
+	__nlRowColumnDestroy(&__nlCurrentContext->al);
+	__nlCurrentContext->alloc_al = NL_FALSE;
+	
 #if 0
-    if(!__nlCurrentContext->least_squares) {
-        __nl_assert(
-            __nlCurrentContext->current_row == 
-            __nlCurrentContext->n
-        );
-    }
+	if(!__nlCurrentContext->least_squares) {
+		__nl_assert(
+			__nlCurrentContext->current_row == 
+			__nlCurrentContext->n
+		);
+	}
 #endif
 }
 
 static void __nlBeginRow() {
-    __nlTransition(__NL_STATE_MATRIX, __NL_STATE_ROW);
-    __nlRowColumnZero(&__nlCurrentContext->af);
-    __nlRowColumnZero(&__nlCurrentContext->al);
-    __nlRowColumnZero(&__nlCurrentContext->xl);
-}
-
-static void __nlScaleRow(NLfloat s) {
-    __NLRowColumn*    af = &__nlCurrentContext->af;
-    __NLRowColumn*    al = &__nlCurrentContext->al;
-    NLuint nf            = af->size;
-    NLuint nl            = al->size;
-    NLuint i;
-    for(i=0; i<nf; i++) {
-        af->coeff[i].value *= s;
-    }
-    for(i=0; i<nl; i++) {
-        al->coeff[i].value *= s;
-    }
-    __nlCurrentContext->right_hand_side *= s;
-}
-
-static void __nlNormalizeRow(NLfloat weight) {
-    __NLRowColumn*    af = &__nlCurrentContext->af;
-    __NLRowColumn*    al = &__nlCurrentContext->al;
-    NLuint nf            = af->size;
-    NLuint nl            = al->size;
-    NLuint i;
-    NLfloat norm = 0.0;
-    for(i=0; i<nf; i++) {
-        norm += af->coeff[i].value * af->coeff[i].value;
-    }
-    for(i=0; i<nl; i++) {
-        norm += al->coeff[i].value * al->coeff[i].value;
-    }
-    norm = sqrt(norm);
-    __nlScaleRow(weight / norm);
+	__nlTransition(__NL_STATE_MATRIX, __NL_STATE_ROW);
+	__nlRowColumnZero(&__nlCurrentContext->af);
+	__nlRowColumnZero(&__nlCurrentContext->al);
 }
 
 static void __nlEndRow() {
-    __NLRowColumn*    af = &__nlCurrentContext->af;
-    __NLRowColumn*    al = &__nlCurrentContext->al;
-    __NLRowColumn*    xl = &__nlCurrentContext->xl;
-    __NLSparseMatrix* M  = &__nlCurrentContext->M;
-    NLfloat* b        = __nlCurrentContext->b;
-    NLuint nf          = af->size;
-    NLuint nl          = al->size;
-    NLuint current_row = __nlCurrentContext->current_row;
-    NLuint i;
-    NLuint j;
-    NLfloat S;
-    __nlTransition(__NL_STATE_ROW, __NL_STATE_MATRIX);
-
-    if(__nlCurrentContext->normalize_rows) {
-        __nlNormalizeRow(__nlCurrentContext->row_scaling);
-    } else {
-        __nlScaleRow(__nlCurrentContext->row_scaling);
-    }
-
-    if(__nlCurrentContext->least_squares) {
-        for(i=0; i<nf; i++) {
-            for(j=0; j<nf; j++) {
-                __nlSparseMatrixAdd(
-                    M, af->coeff[i].index, af->coeff[j].index,
-                    af->coeff[i].value * af->coeff[j].value
-                );
-            }
-        }
-        S = -__nlCurrentContext->right_hand_side;
-        for(j=0; j<nl; j++) {
-            S += al->coeff[j].value * xl->coeff[j].value;
-        }
-        for(i=0; i<nf; i++) {
-            b[ af->coeff[i].index ] -= af->coeff[i].value * S;
-        }
-    } else {
-        for(i=0; i<nf; i++) {
-            __nlSparseMatrixAdd(
-                M, current_row, af->coeff[i].index, af->coeff[i].value
-            );
-        }
-        b[current_row] = -__nlCurrentContext->right_hand_side;
-        for(i=0; i<nl; i++) {
-            b[current_row] -= al->coeff[i].value * xl->coeff[i].value;
-        }
-    }
-    __nlCurrentContext->current_row++;
-    __nlCurrentContext->right_hand_side = 0.0;    
-    __nlCurrentContext->row_scaling     = 1.0;
+	__NLRowColumn*	af = &__nlCurrentContext->af;
+	__NLRowColumn*	al = &__nlCurrentContext->al;
+	__NLSparseMatrix* M  = &__nlCurrentContext->M;
+	NLfloat* b		= __nlCurrentContext->b;
+	NLuint nf		  = af->size;
+	NLuint nl		  = al->size;
+	NLuint current_row = __nlCurrentContext->current_row;
+	NLuint i;
+	NLuint j;
+	NLfloat S;
+	__nlTransition(__NL_STATE_ROW, __NL_STATE_MATRIX);
+
+	if(__nlCurrentContext->least_squares) {
+		if (!__nlCurrentContext->solve_again) {
+			for(i=0; i<nf; i++) {
+				for(j=0; j<nf; j++) {
+					__nlSparseMatrixAdd(
+						M, af->coeff[i].index, af->coeff[j].index,
+						af->coeff[i].value * af->coeff[j].value
+					);
+				}
+			}
+		}
+
+		S = -__nlCurrentContext->right_hand_side;
+		for(j=0; j<nl; j++) {
+			S += al->coeff[j].value;
+		}
+		for(i=0; i<nf; i++) {
+			b[ af->coeff[i].index ] -= af->coeff[i].value * S;
+		}
+	} else {
+		if (!__nlCurrentContext->solve_again) {
+			for(i=0; i<nf; i++) {
+				__nlSparseMatrixAdd(
+					M, current_row, af->coeff[i].index, af->coeff[i].value
+				);
+			}
+		}
+		b[current_row] = -__nlCurrentContext->right_hand_side;
+		for(i=0; i<nl; i++) {
+			b[current_row] -= al->coeff[i].value;
+		}
+	}
+	__nlCurrentContext->current_row++;
+	__nlCurrentContext->right_hand_side = 0.0;	
 }
 
 void nlMatrixAdd(NLuint row, NLuint col, NLfloat value)
 {
-    __NLSparseMatrix* M  = &__nlCurrentContext->M;
-    __nlCheckState(__NL_STATE_MATRIX);
-    __nl_range_assert(row, 0, __nlCurrentContext->n - 1);
-    __nl_range_assert(col, 0, __nlCurrentContext->nb_variables - 1);
+	__NLSparseMatrix* M  = &__nlCurrentContext->M;
+	__nlCheckState(__NL_STATE_MATRIX);
+	__nl_range_assert(row, 0, __nlCurrentContext->n - 1);
+	__nl_range_assert(col, 0, __nlCurrentContext->nb_variables - 1);
 	__nl_assert(!__nlCurrentContext->least_squares);
 
 	__nlSparseMatrixAdd(M, row, col, value);
@@ -952,226 +907,256 @@ void nlMatrixAdd(NLuint row, NLuint col, NLfloat value)
 
 void nlRightHandSideAdd(NLuint index, NLfloat value)
 {
-    NLfloat* b = __nlCurrentContext->b;
+	NLfloat* b = __nlCurrentContext->b;
 
-    __nlCheckState(__NL_STATE_MATRIX);
-    __nl_range_assert(index, 0, __nlCurrentContext->n - 1);
+	__nlCheckState(__NL_STATE_MATRIX);
+	__nl_range_assert(index, 0, __nlCurrentContext->n - 1);
 	__nl_assert(!__nlCurrentContext->least_squares);
 
 	b[index] += value;
 }
 
 void nlCoefficient(NLuint index, NLfloat value) {
-    __NLVariable* v;
+	__NLVariable* v;
 	unsigned int zero= 0;
-    __nlCheckState(__NL_STATE_ROW);
-    __nl_range_assert(index, zero, __nlCurrentContext->nb_variables - 1);
-    v = &(__nlCurrentContext->variable[index]);
-    if(v->locked) {
-        __nlRowColumnAppend(&(__nlCurrentContext->al), 0, value);
-        __nlRowColumnAppend(&(__nlCurrentContext->xl), 0, v->value);
-    } else {
-        __nlRowColumnAppend(&(__nlCurrentContext->af), v->index, value);
-    }
+	__nlCheckState(__NL_STATE_ROW);
+	__nl_range_assert(index, zero, __nlCurrentContext->nb_variables - 1);
+	v = &(__nlCurrentContext->variable[index]);
+	if(v->locked)
+		__nlRowColumnAppend(&(__nlCurrentContext->al), 0, value*v->value);
+	else
+		__nlRowColumnAppend(&(__nlCurrentContext->af), v->index, value);
 }
 
 void nlBegin(NLenum prim) {
-    switch(prim) {
-    case NL_SYSTEM: {
-        __nlBeginSystem();
-    } break;
-    case NL_MATRIX: {
-        __nlBeginMatrix();
-    } break;
-    case NL_ROW: {
-        __nlBeginRow();
-    } break;
-    default: {
-        __nl_assert_not_reached;
-    }
-    }
+	switch(prim) {
+	case NL_SYSTEM: {
+		__nlBeginSystem();
+	} break;
+	case NL_MATRIX: {
+		__nlBeginMatrix();
+	} break;
+	case NL_ROW: {
+		__nlBeginRow();
+	} break;
+	default: {
+		__nl_assert_not_reached;
+	}
+	}
 }
 
 void nlEnd(NLenum prim) {
-    switch(prim) {
-    case NL_SYSTEM: {
-        __nlEndSystem();
-    } break;
-    case NL_MATRIX: {
-        __nlEndMatrix();
-    } break;
-    case NL_ROW: {
-        __nlEndRow();
-    } break;
-    default: {
-        __nl_assert_not_reached;
-    }
-    }
+	switch(prim) {
+	case NL_SYSTEM: {
+		__nlEndSystem();
+	} break;
+	case NL_MATRIX: {
+		__nlEndMatrix();
+	} break;
+	case NL_ROW: {
+		__nlEndRow();
+	} break;
+	default: {
+		__nl_assert_not_reached;
+	}
+	}
 }
 
 /************************************************************************/
 /* SuperLU wrapper */
 
-/* Note: SuperLU is difficult to call, but it is worth it.    */
+/* Note: SuperLU is difficult to call, but it is worth it.	*/
 /* Here is a driver inspired by A. Sheffer's "cow flattener". */
-static NLboolean __nlSolve_SUPERLU( NLboolean do_perm) {
-
-    /* OpenNL Context */
-    __NLSparseMatrix* M  = &(__nlCurrentContext->M);
-    NLfloat* b          = __nlCurrentContext->b;
-    NLfloat* x          = __nlCurrentContext->x;
-
-    /* Compressed Row Storage matrix representation */
-    NLuint    n      = __nlCurrentContext->n;
-    NLuint    nnz    = __nlSparseMatrixNNZ(M); /* Number of Non-Zero coeffs */
-    NLint*    xa     = __NL_NEW_ARRAY(NLint, n+1);
-    NLfloat* rhs    = __NL_NEW_ARRAY(NLfloat, n);
-    NLfloat* a      = __NL_NEW_ARRAY(NLfloat, nnz);
-    NLint*    asub   = __NL_NEW_ARRAY(NLint, nnz);
-
-    /* Permutation vector */
-    NLint*    perm_r  = __NL_NEW_ARRAY(NLint, n);
-    NLint*    perm    = __NL_NEW_ARRAY(NLint, n);
-
-    /* SuperLU variables */
-    SuperMatrix A, B; /* System       */
-    SuperMatrix L, U; /* Inverse of A */
-    NLint info;       /* status code  */
-    DNformat *vals = NULL; /* access to result */
-    float *rvals  = NULL; /* access to result */
-
-    /* SuperLU options and stats */
-    superlu_options_t options;
-    SuperLUStat_t     stat;
-
-
-    /* Temporary variables */
-    __NLRowColumn* Ri = NULL;
-    NLuint         i,jj,count;
-    
-    __nl_assert(!(M->storage & __NL_SYMMETRIC));
-    __nl_assert(M->storage & __NL_ROWS);
-    __nl_assert(M->m == M->n);
-    
-    
-    /*
-     * Step 1: convert matrix M into SuperLU compressed column 
-     *   representation.
-     * -------------------------------------------------------
-     */
-
-    count = 0;
-    for(i=0; i<n; i++) {
-        Ri = &(M->row[i]);
-        xa[i] = count;
-        for(jj=0; jj<Ri->size; jj++) {
-            a[count]    = Ri->coeff[jj].value;
-            asub[count] = Ri->coeff[jj].index;
-            count++;
-        }
-    }
-    xa[n] = nnz;
-
-    /* Save memory for SuperLU */
-    __nlSparseMatrixClear(M);
-
-
-    /*
-     * Rem: symmetric storage does not seem to work with
-     * SuperLU ... (->deactivated in main SLS::Solver driver)
-     */
-    sCreate_CompCol_Matrix(
-        &A, n, n, nnz, a, asub, xa, 
-        SLU_NR,              /* Row_wise, no supernode */
-        SLU_S,               /* floats                */ 
-        SLU_GE               /* general storage        */
-    );
-
-    /* Step 2: create vector */
-    sCreate_Dense_Matrix(
-        &B, n, 1, b, n, 
-        SLU_DN, /* Fortran-type column-wise storage */
-        SLU_S,  /* floats                          */
-        SLU_GE  /* general                          */
-    );
-            
-
-    /* Step 3: get permutation matrix 
-     * ------------------------------
-     * com_perm: 0 -> no re-ordering
-     *           1 -> re-ordering for A^t.A
-     *           2 -> re-ordering for A^t+A
-     *           3 -> approximate minimum degree ordering
-     */
-    get_perm_c(do_perm ? 3 : 0, &A, perm);
-
-    /* Step 4: call SuperLU main routine
-     * ---------------------------------
-     */
-
-    set_default_options(&options);
-    options.ColPerm = MY_PERMC;
-    StatInit(&stat);
-
-    sgssv(&options, &A, perm, perm_r, &L, &U, &B, &stat, &info);
-
-    /* Step 5: get the solution
-     * ------------------------
-     * Fortran-type column-wise storage
-     */
-    vals = (DNformat*)B.Store;
-    rvals = (float*)(vals->nzval);
-    if(info == 0) {
-        for(i = 0; i <  n; i++){
-            x[i] = rvals[i];
-        }
-    }
-
-    /* Step 6: cleanup
-     * ---------------
-     */
-
-    /*
-     *  For these two ones, only the "store" structure
-     * needs to be deallocated (the arrays have been allocated
-     * by us).
-     */
-    Destroy_SuperMatrix_Store(&A);
-    Destroy_SuperMatrix_Store(&B);
-
-    
-    /*
-     *   These ones need to be fully deallocated (they have been
-     * allocated by SuperLU).
-     */
-    Destroy_SuperNode_Matrix(&L);
-    Destroy_CompCol_Matrix(&U);
-
-    StatFree(&stat);
-
-    __NL_DELETE_ARRAY(xa);
-    __NL_DELETE_ARRAY(rhs);
-    __NL_DELETE_ARRAY(a);
-    __NL_DELETE_ARRAY(asub);
-    __NL_DELETE_ARRAY(perm_r);
-    __NL_DELETE_ARRAY(perm);
-
-    return (info == 0);
+static NLboolean __nlFactorize_SUPERLU(__NLContext *context, NLint *permutation) {
+
+	/* OpenNL Context */
+	__NLSparseMatrix* M = &(context->M);
+	NLuint n = context->n;
+	NLuint nnz = __nlSparseMatrixNNZ(M); /* number of non-zero coeffs */
+
+	/* Compressed Row Storage matrix representation */
+	NLint	*xa		= __NL_NEW_ARRAY(NLint, n+1);
+	NLfloat	*rhs	= __NL_NEW_ARRAY(NLfloat, n);
+	NLfloat	*a		= __NL_NEW_ARRAY(NLfloat, nnz);
+	NLint	*asub	= __NL_NEW_ARRAY(NLint, nnz);
+	NLint	*etree	= __NL_NEW_ARRAY(NLint, n);
+
+	/* SuperLU variables */
+	SuperMatrix At, AtP;
+	NLint info, panel_size, relax;
+	superlu_options_t options;
+
+	/* Temporary variables */
+	__NLRowColumn* Ri = NULL;
+	NLuint i, jj, count;
+	
+	__nl_assert(!(M->storage & __NL_SYMMETRIC));
+	__nl_assert(M->storage & __NL_ROWS);
+	__nl_assert(M->m == M->n);
+	
+	/* Convert M to compressed column format */
+	for(i=0, count=0; i<n; i++) {
+		__NLRowColumn *Ri = M->row + i;
+		xa[i] = count;
+
+		for(jj=0; jj<Ri->size; jj++, count++) {
+			a[count] = Ri->coeff[jj].value;
+			asub[count] = Ri->coeff[jj].index;
+		}
+	}
+	xa[n] = nnz;
+
+	/* Free M, don't need it anymore at this point */
+	__nlSparseMatrixClear(M);
+
+	/* Create superlu A matrix transposed */
+	sCreate_CompCol_Matrix(
+		&At, n, n, nnz, a, asub, xa, 
+		SLU_NC,		/* Colum wise, no supernode */
+		SLU_S,		/* floats */ 
+		SLU_GE		/* general storage */
+	);
+
+	/* Set superlu options */
+	set_default_options(&options);
+	options.ColPerm = MY_PERMC;
+	options.Fact = DOFACT;
+
+	StatInit(&(context->slu.stat));
+
+	panel_size = sp_ienv(1); /* sp_ienv give us the defaults */
+	relax = sp_ienv(2);
+
+	/* Compute permutation and permuted matrix */
+	context->slu.perm_r = __NL_NEW_ARRAY(NLint, n);
+	context->slu.perm_c = __NL_NEW_ARRAY(NLint, n);
+
+	if ((permutation == NULL) || (*permutation == -1)) {
+		get_perm_c(3, &At, context->slu.perm_c);
+
+		if (permutation)
+			memcpy(permutation, context->slu.perm_c, sizeof(NLint)*n);
+	}
+	else
+		memcpy(context->slu.perm_c, permutation, sizeof(NLint)*n);
+
+	sp_preorder(&options, &At, context->slu.perm_c, etree, &AtP);
+
+	/* Decompose into L and U */
+	sgstrf(&options, &AtP, relax, panel_size,
+		etree, NULL, 0, context->slu.perm_c, context->slu.perm_r,
+		&(context->slu.L), &(context->slu.U), &(context->slu.stat), &info);
+
+	/* Cleanup */
+
+	Destroy_SuperMatrix_Store(&At);
+	Destroy_SuperMatrix_Store(&AtP);
+
+	__NL_DELETE_ARRAY(etree);
+	__NL_DELETE_ARRAY(xa);
+	__NL_DELETE_ARRAY(rhs);
+	__NL_DELETE_ARRAY(a);
+	__NL_DELETE_ARRAY(asub);
+
+	context->slu.alloc_slu = NL_TRUE;
+
+	return (info == 0);
+}
+
+static NLboolean __nlInvert_SUPERLU(__NLContext *context) {
+
+	/* OpenNL Context */
+	NLfloat* b = context->b;
+	NLfloat* x = context->x;
+	NLuint n = context->n;
+
+	/* SuperLU variables */
+	SuperMatrix B;
+	NLint info;
+
+	/* Create superlu array for B */
+	sCreate_Dense_Matrix(
+		&B, n, 1, b, n, 
+		SLU_DN, /* Fortran-type column-wise storage */
+		SLU_S,  /* floats						  */
+		SLU_GE  /* general						  */
+	);
+
+	/* Forward/Back substitution to compute x */
+	sgstrs(TRANS, &(context->slu.L), &(context->slu.U),
+		context->slu.perm_c, context->slu.perm_r, &B,
+		&(context->slu.stat), &info);
+
+	if(info == 0)
+		memcpy(x, ((DNformat*)B.Store)->nzval, sizeof(*x)*n);
+
+	Destroy_SuperMatrix_Store(&B);
+
+	return (info == 0);
+}
+
+static void __nlFree_SUPERLU(__NLContext *context) {
+
+	Destroy_SuperNode_Matrix(&(context->slu.L));
+	Destroy_CompCol_Matrix(&(context->slu.U));
+
+	StatFree(&(context->slu.stat));
+
+	__NL_DELETE_ARRAY(context->slu.perm_r);
+	__NL_DELETE_ARRAY(context->slu.perm_c);
+
+	context->slu.alloc_slu = NL_FALSE;
 }
 
+void nlPrintMatrix(void) {
+	__NLSparseMatrix* M  = &(__nlCurrentContext->M);
+	NLuint i, jj, k;
+	NLuint n = __nlCurrentContext->n;
+	__NLRowColumn* Ri = NULL;
+	float *value = malloc(sizeof(*value)*n);
+
+	printf("M:\n");
+	for(i=0; i<n; i++) {
+		Ri = &(M->row[i]);
+
+		memset(value, 0.0, sizeof(*value)*n);
+		for(jj=0; jj<Ri->size; jj++)
+			value[Ri->coeff[jj].index] = Ri->coeff[jj].value;
+
+		for (k = 0; k<n; k++)
+			printf("%.3f ", value[k]);
+		printf("\n");
+	}
+}
 
 /************************************************************************/
 /* nlSolve() driver routine */
 
-NLboolean nlSolve(void) {
-    NLboolean result = NL_TRUE;
+NLboolean nlSolveAdvanced(NLint *permutation, NLboolean solveAgain) {
+	NLboolean result = NL_TRUE;
+
+	__nlCheckState(__NL_STATE_SYSTEM_CONSTRUCTED);
+
+	if (!__nlCurrentContext->solve_again)
+		result = __nlFactorize_SUPERLU(__nlCurrentContext, permutation);
+
+	if (result) {
+		result = __nlInvert_SUPERLU(__nlCurrentContext);
 
-    __nlCheckState(__NL_STATE_SYSTEM_CONSTRUCTED);
-    result = __nlSolve_SUPERLU(NL_TRUE);
+		if (result) {
+			__nlVectorToVariables();
 
-    __nlVectorToVariables();
-    __nlTransition(__NL_STATE_SYSTEM_CONSTRUCTED, __NL_STATE_SOLVED);
+			if (solveAgain)
+				__nlCurrentContext->solve_again = NL_TRUE;
+
+			__nlTransition(__NL_STATE_SYSTEM_CONSTRUCTED, __NL_STATE_SYSTEM_SOLVED);
+		}
+	}
+
+	return result;
+}
 
-    return result;
+NLboolean nlSolve() {
+	return nlSolveAdvanced(NULL, NL_FALSE);
 }
 
diff --git a/source/blender/include/BDR_unwrapper.h b/source/blender/include/BDR_unwrapper.h
index 176bf8c8e68..e853b48bc9d 100644
--- a/source/blender/include/BDR_unwrapper.h
+++ b/source/blender/include/BDR_unwrapper.h
@@ -37,6 +37,7 @@ void set_seamtface(void); /* set TF_SEAM flags in tfaces */
 void unwrap_lscm(void); /* unwrap selected tfaces */
 void unwrap_lscm_live(void); /* unwrap selected tfaces (for live mode, with no undo pushes) */
 void select_linked_tfaces_with_seams(int mode, Mesh *me, unsigned int index);
+void minimize_stretch_tface_uv(void);
 
 #endif /* BDR_UNWRAPPER_H */
 
diff --git a/source/blender/src/header_image.c b/source/blender/src/header_image.c
index e0571e94e45..9d12e417d02 100644
--- a/source/blender/src/header_image.c
+++ b/source/blender/src/header_image.c
@@ -1013,6 +1013,9 @@ static void do_image_uvsmenu(void *arg, int event)
 		if(G.sima->flag & SI_LSCM_LIVE) G.sima->flag &= ~SI_LSCM_LIVE;
 		else G.sima->flag |= SI_LSCM_LIVE;
 		break;
+	case 12:
+		minimize_stretch_tface_uv();
+		break;
 	}
 }
 
@@ -1047,6 +1050,7 @@ static uiBlock *image_uvsmenu(void *arg_unused)
 
 	uiDefBut(block, SEPR, 0, "", 0, yco-=6, menuwidth, 6, NULL, 0.0, 0.0, 0, 0, "");	
 
+	uiDefIconTextBut(block, BUTM, 1, ICON_BLANK1, "Minimize Stretch|Ctrl V", 0, yco-=20, menuwidth, 19, NULL, 0.0, 0.0, 0, 3, "");
 	uiDefIconTextBut(block, BUTM, 1, ICON_BLANK1, "Limit Stitch...|Shift V", 0, yco-=20, menuwidth, 19, NULL, 0.0, 0.0, 0, 3, "");
 	uiDefIconTextBut(block, BUTM, 1, ICON_BLANK1, "Stitch|V", 0, yco-=20, menuwidth, 19, NULL, 0.0, 0.0, 0, 4, "");
 	uiDefIconTextBlockBut(block, image_uvs_transformmenu, NULL, ICON_RIGHTARROW_THIN, "Transform", 0, yco-=20, 120, 19, "");
diff --git a/source/blender/src/parametrizer.c b/source/blender/src/parametrizer.c
new file mode 100644
index 00000000000..ff5d92e4604
--- /dev/null
+++ b/source/blender/src/parametrizer.c
@@ -0,0 +1,1785 @@
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_memarena.h"
+#include "BLI_arithb.h"
+#include "BLI_rand.h"
+
+#include "BKE_utildefines.h"
+
+#include "BIF_editsima.h"
+#include "BIF_toolbox.h"
+
+#include "ONL_opennl.h"
+
+#include "parametrizer.h"
+#include "parametrizer_intern.h"
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+#if defined(_WIN32)
+#define M_PI 3.14159265358979323846
+#endif
+
+/* Hash */
+
+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) (((unsigned long) (item))%((unsigned int) (ph)->cursize))
+
+PHash *phash_new(int sizehint)
+{
+	PHash *ph = (PHash*)MEM_callocN(sizeof(PHash), "PHash");
+	ph->size = 0;
+	ph->cursize_id = 0;
+	ph->first = NULL;
+
+	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;
+}
+
+void phash_delete(PHash *ph)
+{
+	MEM_freeN(ph->buckets);
+	MEM_freeN(ph);
+}
+
+void phash_delete_with_links(PHash *ph)
+{
+	PHashLink *link, *next=NULL;
+
+	for (link = ph->first; link; link = next) {
+		next = link->next;
+		MEM_freeN(link);
+	}
+
+	phash_delete(ph);
+}
+
+int phash_size(PHash *ph)
+{
+	return ph->size;
+}
+
+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->first;
+		ph->first = 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->first;
+
+		ph->cursize = PHashSizes[++ph->cursize_id];
+		MEM_freeN(ph->buckets);
+		ph->buckets = (PHashLink**)MEM_callocN(ph->cursize*sizeof(*ph->buckets), "PHashBuckets");
+		ph->size = 0;
+		ph->first = NULL;
+
+		for (link = first; link; link = next) {
+			next = link->next;
+			phash_insert(ph, link);
+		}
+	}
+}
+
+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;
+}
+
+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;
+}
+
+/* Heap */
+
+#define PHEAP_PARENT(i) ((i-1)>>1)
+#define PHEAP_LEFT(i)   ((i<<1)+1)
+#define PHEAP_RIGHT(i)  ((i<<1)+2)
+#define PHEAP_COMPARE(a, b) (a->value < b->value)
+#define PHEAP_EQUALS(a, b) (a->value == b->value)
+#define PHEAP_SWAP(heap, i, j) \
+	{ SWAP(int, heap->tree[i]->index, heap->tree[j]->index); \
+	  SWAP(PHeapLink*, heap->tree[i], heap->tree[j]);  }
+
+static void pheap_down(PHeap *heap, int i)
+{
+	while (P_TRUE) {
+		int size = heap->size, smallest;
+		int l = PHEAP_LEFT(i);
+		int r = PHEAP_RIGHT(i);
+
+		smallest = ((l < size) && PHEAP_COMPARE(heap->tree[l], heap->tree[i]))? l: i;
+
+		if ((r < size) && PHEAP_COMPARE(heap->tree[r], heap->tree[smallest]))
+			smallest = r;
+		
+		if (smallest == i)
+			break;
+
+		PHEAP_SWAP(heap, i, smallest);
+		i = smallest;
+	}
+}
+
+static void pheap_up(PHeap *heap, int i)
+{
+	while (i > 0) {
+		int p = PHEAP_PARENT(i);
+
+		if (PHEAP_COMPARE(heap->tree[p], heap->tree[i]))
+			break;
+
+		PHEAP_SWAP(heap, p, i);
+		i = p;
+	}
+}
+
+PHeap *pheap_new()
+{
+	/* TODO: replace mallocN with something faster */
+
+	PHeap *heap = (PHeap*)MEM_callocN(sizeof(PHeap), "PHeap");
+	heap->bufsize = 1;
+	heap->tree = (PHeapLink**)MEM_mallocN(sizeof(PHeapLink*), "PHeapTree");
+
+	return heap;
+}
+
+void pheap_delete(PHeap *heap)
+{
+	MEM_freeN(heap->tree);
+	MEM_freeN(heap);
+}
+
+PHeapLink *pheap_insert(PHeap *heap, float value, void *ptr)
+{
+	PHeapLink *link;
+
+	if ((heap->size + 1) > heap->bufsize) {
+		int newsize = heap->bufsize*2;
+
+		PHeapLink **ntree = (PHeapLink**)MEM_mallocN(newsize*sizeof(PHeapLink*), "PHeapTree");
+		memcpy(ntree, heap->tree, sizeof(PHeapLink*)*heap->size);
+		MEM_freeN(heap->tree);
+
+		heap->tree = ntree;
+		heap->bufsize = newsize;
+	}
+
+	param_assert(heap->size < heap->bufsize);
+
+	link = MEM_mallocN(sizeof *link, "PHeapLink");
+	link->value = value;
+	link->ptr = ptr;
+	link->index = heap->size;
+
+	heap->tree[link->index] = link;
+
+	heap->size++;
+
+	pheap_up(heap, heap->size-1);
+
+	return link;
+}
+
+int pheap_empty(PHeap *heap)
+{
+	return (heap->size == 0);
+}
+
+int pheap_size(PHeap *heap)
+{
+	return heap->size;
+}
+
+void *pheap_min(PHeap *heap)
+{
+	return heap->tree[0]->ptr;
+}
+
+void *pheap_popmin(PHeap *heap)
+{
+	void *ptr = heap->tree[0]->ptr;
+
+	MEM_freeN(heap->tree[0]);
+
+	if (heap->size == 1)
+		heap->size--;
+	else {
+		PHEAP_SWAP(heap, 0, heap->size-1);
+		heap->size--;
+
+		pheap_down(heap, 0);
+	}
+
+	return ptr;
+}
+
+static void pheap_remove(PHeap *heap, PHeapLink *link)
+{
+	int i = link->index;
+
+	while (i > 0) {
+		int p = PHEAP_PARENT(i);
+
+		PHEAP_SWAP(heap, p, i);
+		i = p;
+	}
+
+	pheap_popmin(heap);
+}
+
+/* Construction */
+
+PEdge *p_wheel_edge_next(PEdge *e)
+{
+	return e->next->next->pair;
+}
+
+PEdge *p_wheel_edge_prev(PEdge *e)
+{   
+	return (e->pair)? e->pair->next: NULL;
+}
+
+static PVert *p_vert_add(PChart *chart, PHashKey key, float *co, PEdge *e)
+{
+	PVert *v = (PVert*)BLI_memarena_alloc(chart->handle->arena, sizeof *v);
+	v->co = co;
+	v->link.key = key;
+	v->edge = e;
+
+	phash_insert(chart->verts, (PHashLink*)v);
+
+	return v;
+}
+
+static PVert *p_vert_lookup(PChart *chart, PHashKey key, float *co, PEdge *e)
+{
+	PVert *v = (PVert*)phash_lookup(chart->verts, key);
+
+	if (v)
+		return v;
+	else
+		return p_vert_add(chart, 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->link.key = v->link.key;
+	nv->edge = v->edge;
+
+	phash_insert(chart->verts, (PHashLink*)nv);
+
+	return nv;
+}
+
+static PEdge *p_edge_lookup(PChart *chart, PHashKey *vkeys)
+{
+	PHashKey key = vkeys[0]^vkeys[1];
+	PEdge *e = (PEdge*)phash_lookup(chart->edges, key);
+
+	while (e) {
+		if ((e->vert->link.key == vkeys[0]) && (e->next->vert->link.key == vkeys[1]))
+			return e;
+		else if ((e->vert->link.key == vkeys[1]) && (e->next->vert->link.key == vkeys[0]))
+			return e;
+
+		e = (PEdge*)phash_next(chart->edges, key, (PHashLink*)e);
+	}
+
+	return NULL;
+}
+
+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);
+}
+
+static void p_vert_load_pin_uvs(PVert *v)
+{
+	PEdge *e;
+	int nedges = 0;
+
+	v->uv[0] = v->uv[1] = 0.0f;
+	nedges = 0;
+	e = v->edge;
+	do {
+		if (e->orig_uv && (e->flag & PEDGE_PIN)) {
+			v->flag |= PVERT_PIN;
+			v->uv[0] += e->orig_uv[0];
+			v->uv[1] += e->orig_uv[1];
+			nedges++;
+		}
+
+		e = p_wheel_edge_next(e);
+	} while (e && e != (v->edge));
+
+	if (nedges > 0) {
+		v->uv[0] /= nedges;
+		v->uv[1] /= nedges;
+	}
+}
+
+static void p_vert_load_select_uvs(PVert *v)
+{
+	PEdge *e;
+	int nedges = 0;
+
+	v->uv[0] = v->uv[1] = 0.0f;
+	nedges = 0;
+	e = v->edge;
+	do {
+		if (e->orig_uv && (e->flag & PEDGE_SELECT))
+			v->flag |= PVERT_SELECT;
+
+		v->uv[0] += e->orig_uv[0];
+		v->uv[1] += e->orig_uv[1];
+		nedges++;
+
+		e = p_wheel_edge_next(e);
+	} while (e && e != (v->edge));
+
+	if (nedges > 0) {
+		v->uv[0] /= nedges;
+		v->uv[1] /= nedges;
+	}
+}
+
+static void p_extrema_verts(PChart *chart, PVert **v1, PVert **v2)
+{
+	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 = (PVert*)chart->verts->first; v; v=v->link.next) {
+		for (i = 0; i < 3; i++) {
+			if (v->co[i] < minv[i]) {
+				minv[i] = v->co[i];
+				minvert[i] = v;
+			}
+			if (v->co[i] > maxv[i]) {
+				maxv[i] = v->co[i];
+				maxvert[i] = v;
+			}
+		}
+	}
+
+	/* 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];
+		}
+	}
+
+	if (minvert[dir] == maxvert[dir]) {
+		/* degenerate case */
+		PFace *f = (PFace*)chart->faces->first;
+		*v1 = f->edge->vert;
+		*v2 = f->edge->next->vert;
+	}
+	else {
+		*v1 = minvert[dir];
+		*v2 = maxvert[dir];
+	}
+}
+
+static float p_vec_normalise(float *v)
+{
+   float d;
+   
+    d = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
+
+    if(d != 0.0f) {
+		d = 1.0f/d;
+
+		v[0] *= d;
+		v[1] *= d;
+		v[2] *= d;
+	}
+
+	return d;
+}
+
+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];
+
+	p_vec_normalise(d1);
+	p_vec_normalise(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 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;
+
+	*a1 = p_vec_angle(v3->co, v1->co, v2->co);
+	*a2 = p_vec_angle(v1->co, v2->co, v3->co);
+	*a3 = M_PI - *a2 - *a1;
+}
+
+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_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_face_uv_area(PFace *f)
+{
+	return fabs(p_face_uv_area_signed(f));
+}
+
+static void p_chart_area(PChart *chart, float *uv_area, float *area)
+{
+	PFace *f;
+
+	*uv_area = *area = 0.0f;
+
+	for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
+		*uv_area += p_face_uv_area(f);
+		*area += p_face_area(f);
+	}
+}
+
+static PChart *p_chart_new(PHandle *handle)
+{
+	PChart *chart = (PChart*)MEM_callocN(sizeof*chart, "PChart");
+	chart->verts = phash_new(1);
+	chart->edges = phash_new(1);
+	chart->faces = phash_new(1);
+	chart->handle = handle;
+
+	return chart;
+}
+
+static void p_chart_delete(PChart *chart)
+{
+	/* the actual links are free by memarena */
+	phash_delete(chart->verts);
+	phash_delete(chart->edges);
+	phash_delete(chart->faces);
+
+	MEM_freeN(chart);
+}
+
+static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
+{
+	float *uv1, *uv2, *uvp1, *uvp2;
+	float limit[2];
+
+	uv1 = e->orig_uv;
+	uv2 = e->next->orig_uv;
+
+	if (e->vert->link.key == ep->vert->link.key) {
+		uvp1 = ep->orig_uv;
+		uvp2 = ep->next->orig_uv;
+	}
+	else {
+		uvp1 = ep->next->orig_uv;
+		uvp2 = ep->orig_uv;
+	}
+
+	get_connected_limit_tface_uv(limit);
+
+	if((fabs(uv1[0]-uvp1[0]) > limit[0]) && (fabs(uv1[1]-uvp1[1]) > limit[1])) {
+		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(PChart *chart, PEdge *e, PEdge **pair, PBool impl)
+{
+	PHashKey key;
+	PEdge *pe;
+	PVert *v1, *v2;
+	PHashKey key1 = e->vert->link.key;
+	PHashKey key2 = e->next->vert->link.key;
+
+	if (e->flag & PEDGE_SEAM)
+		return P_FALSE;
+	
+	key = key1 ^ key2;
+	pe = (PEdge*)phash_lookup(chart->edges, key);
+	*pair = NULL;
+
+	while (pe) {
+		if (pe != e) {
+			v1 = pe->vert;
+			v2 = pe->next->vert;
+
+			if (((v1->link.key == key1) && (v2->link.key == key2)) ||
+				((v1->link.key == key2) && (v2->link.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(chart->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(PChart *chart, PEdge *e, PEdge ***stack, PBool impl)
+{
+	PEdge *pair;
+
+	if(!e->pair && p_edge_has_pair(chart, 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(PChart *chart, PBool impl)
+{
+	PEdge **stackbase = MEM_mallocN(sizeof*stackbase * phash_size(chart->faces), "Pstackbase");
+	PEdge **stack = stackbase;
+	PFace *f, *first;
+	PEdge *e, *e1, *e2;
+	int ncharts = 0;
+
+	/* connect pairs, count edges, set vertex-edge pointer to a pairless edge */
+	for (first=(PFace*)chart->faces->first; first; first=first->link.next) {
+		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(chart, e, &stack, impl))
+				e->vert->edge = e;
+			if (!p_edge_connect_pair(chart, e1, &stack, impl))
+				e1->vert->edge = e1;
+			if (!p_edge_connect_pair(chart, 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;
+			phash_insert(chart->verts, (PHashLink*)v);
+		}
+	}
+
+	if (copy) {
+		/* not found, copying */
+		v = p_vert_copy(chart, v);
+		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 = (PFace*)chart->faces->first;
+	while (f) {
+		PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+		nextf = f->link.next;
+
+		nchart = charts[f->u.chart];
+
+		phash_insert(nchart->faces, (PHashLink*)f);
+		phash_insert(nchart->edges, (PHashLink*)e1);
+		phash_insert(nchart->edges, (PHashLink*)e2);
+		phash_insert(nchart->edges, (PHashLink*)e3);
+
+		p_split_vert(nchart, e1);
+		p_split_vert(nchart, e2);
+		p_split_vert(nchart, e3);
+
+		f = nextf;
+	}
+
+	return charts;
+}
+
+static void p_face_backup_uvs(PFace *f)
+{
+	PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+
+	e1->old_uv[0] = e1->orig_uv[0];
+	e1->old_uv[1] = e1->orig_uv[1];
+	e2->old_uv[0] = e2->orig_uv[0];
+	e2->old_uv[1] = e2->orig_uv[1];
+	e3->old_uv[0] = e3->orig_uv[0];
+	e3->old_uv[1] = e3->orig_uv[1];
+}
+
+static void p_face_restore_uvs(PFace *f)
+{
+	PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+
+	e1->orig_uv[0] = e1->old_uv[0];
+	e1->orig_uv[1] = e1->old_uv[1];
+	e2->orig_uv[0] = e2->old_uv[0];
+	e2->orig_uv[1] = e2->old_uv[1];
+	e3->orig_uv[0] = e3->old_uv[0];
+	e3->orig_uv[1] = e3->old_uv[1];
+}
+
+static PFace *p_face_add(PChart *chart, ParamKey key, ParamKey *vkeys,
+                         float *co[3], float *uv[3], int i1, int i2, int i3,
+                         ParamBool *pin, ParamBool *select)
+{
+	PFace *f;
+	PEdge *e1, *e2, *e3;
+
+	/* allocate */
+	f = (PFace*)BLI_memarena_alloc(chart->handle->arena, sizeof *f);
+
+	e1 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e1);
+	e2 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e2);
+	e3 = (PEdge*)BLI_memarena_alloc(chart->handle->arena, sizeof *e3);
+
+	/* set up edges */
+	f->edge = e1;
+	e1->face = e2->face = e3->face = f;
+
+	e1->next = e2;
+	e2->next = e3;
+	e3->next = e1;
+
+	if (co && uv) {
+		e1->vert = p_vert_lookup(chart, vkeys[i1], co[i1], e1);
+		e2->vert = p_vert_lookup(chart, vkeys[i2], co[i2], e2);
+		e3->vert = p_vert_lookup(chart, vkeys[i3], co[i3], e3);
+
+		e1->orig_uv = uv[i1];
+		e2->orig_uv = uv[i2];
+		e3->orig_uv = uv[i3];
+	}
+	else {
+		/* internal call to add face */
+		e1->vert = e2->vert = e3->vert = NULL;
+		e1->orig_uv = e2->orig_uv = e3->orig_uv = NULL;
+	}
+
+	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->link.key = key;
+	phash_insert(chart->faces, (PHashLink*)f);
+
+	e1->link.key = vkeys[i1]^vkeys[i2];
+	e2->link.key = vkeys[i2]^vkeys[i3];
+	e3->link.key = vkeys[i3]^vkeys[i1];
+
+	phash_insert(chart->edges, (PHashLink*)e1);
+	phash_insert(chart->edges, (PHashLink*)e2);
+	phash_insert(chart->edges, (PHashLink*)e3);
+
+	return f;
+}
+
+static PBool p_quad_split_direction(float **co)
+{
+    float a1, a2;
+	
+	a1 = p_vec_angle_cos(co[0], co[1], co[2]);
+	a1 += p_vec_angle_cos(co[1], co[0], co[2]);
+	a1 += p_vec_angle_cos(co[2], co[0], co[1]);
+
+	a2 = p_vec_angle_cos(co[0], co[1], co[3]);
+	a2 += p_vec_angle_cos(co[1], co[0], co[3]);
+	a2 += p_vec_angle_cos(co[3], co[0], co[1]);
+
+	return (a1 > a2);
+}
+
+static float p_edge_length(PEdge *e)
+{
+    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]);
+}
+
+void p_chart_uv_bbox(PChart *chart, float *minv, float *maxv)
+{
+	PVert *v;
+
+	INIT_MINMAX2(minv, maxv);
+
+	for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
+		DO_MINMAX2(v->uv, minv, maxv);
+	}
+}
+
+static void p_chart_uv_scale(PChart *chart, float scale)
+{
+	PVert *v;
+
+	for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
+		v->uv[0] *= scale;
+		v->uv[1] *= scale;
+	}
+}
+
+static void p_chart_uv_translate(PChart *chart, float trans[2])
+{
+	PVert *v;
+
+	for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
+		v->uv[0] += trans[0];
+		v->uv[1] += trans[1];
+	}
+}
+
+static void p_chart_boundaries(PChart *chart, int *nboundaries, PEdge **outer)
+{   
+    PEdge *e, *be;
+    float len, maxlen = -1.0;
+
+	*nboundaries = 0;
+	*outer = NULL;
+
+	for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
+        if (e->pair || (e->flag & PEDGE_DONE))
+            continue;
+
+		(*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 (len > maxlen) {
+			*outer = e;
+            maxlen = len;
+        }
+    }
+
+	for (e=(PEdge*)chart->edges->first; e; e=e->link.next)
+        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 PEdge *p_boundary_edge_next(PEdge *e)
+{
+	return e->next->vert->edge;
+}
+
+static PEdge *p_boundary_edge_prev(PEdge *e)
+{
+    PEdge *we = e, *last;
+
+	do {
+		last = we;
+		we = p_wheel_edge_next(we);
+	} while (we && (we != e));
+
+	return last->next->next;
+}
+
+static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
+{
+	PEdge *e, *e1, *e2;
+	PHashKey vkeys[3];
+	PFace *f;
+	struct PHeap *heap = pheap_new(nedges);
+	float angle;
+
+	e = be;
+	do {
+		angle = p_edge_boundary_angle(e);
+		e->u.heaplink = pheap_insert(heap, angle, e);
+
+		e = e->next->vert->edge;
+	} while(e != be);
+
+	if (nedges == 2) {
+		/* no real boundary, but an isolated seam */
+		e = be->next->vert->edge;
+		e->pair = be;
+		be->pair = e;
+
+		pheap_remove(heap, e->u.heaplink);
+		pheap_remove(heap, be->u.heaplink);
+	}
+	else {
+		while (nedges > 2) {
+			PEdge *ne, *ne1, *ne2;
+
+			e = pheap_popmin(heap);
+
+			e1 = p_boundary_edge_prev(e);
+			e2 = p_boundary_edge_next(e);
+
+			pheap_remove(heap, e1->u.heaplink);
+			pheap_remove(heap, e2->u.heaplink);
+			e->u.heaplink = e1->u.heaplink = e2->u.heaplink = NULL;
+
+			e->flag |= PEDGE_FILLED;
+			e1->flag |= PEDGE_FILLED;
+
+			vkeys[0] = e->vert->link.key;
+			vkeys[1] = e1->vert->link.key;
+			vkeys[2] = e2->vert->link.key;
+
+			f = p_face_add(chart, -1, vkeys, NULL, NULL, 0, 1, 2, NULL, NULL);
+			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 = pheap_insert(heap, p_edge_boundary_angle(ne2), ne2);
+				e2->u.heaplink = pheap_insert(heap, p_edge_boundary_angle(e2), e2);
+			}
+
+			nedges--;
+		}
+	}
+
+	pheap_delete(heap);
+}
+
+static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
+{
+	PEdge *e, *enext, *be;
+	int nedges;
+
+	for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
+		enext = e->link.next;
+
+        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);
+    }
+}
+
+static void p_flush_uvs(PChart *chart)
+{
+	PEdge *e;
+
+	for (e=(PEdge*)chart->edges->first; e; e=e->link.next) {
+		if (e->orig_uv) {
+			e->orig_uv[0] = e->vert->uv[0];
+			e->orig_uv[1] = e->vert->uv[1];
+		}
+	}
+}
+
+/* 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));
+	
+	return (ParamHandle*)handle;
+}
+
+void param_delete(ParamHandle *handle)
+{
+	PHandle *phandle = (PHandle*)handle;
+	int i;
+
+	param_assert((phandle->state == PHANDLE_STATE_ALLOCATED) ||
+	             (phandle->state == PHANDLE_STATE_CONSTRUCTED));
+
+	for (i = 0; i < phandle->ncharts; i++)
+		p_chart_delete(phandle->charts[i]);
+	
+	if (phandle->charts)
+		MEM_freeN(phandle->charts);
+
+	if (phandle->construction_chart)
+		p_chart_delete(phandle->construction_chart);
+
+	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;
+	PChart *chart = phandle->construction_chart;
+
+	param_assert(phash_lookup(chart->faces, key) == NULL);
+	param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
+	param_assert((nverts == 3) || (nverts == 4));
+
+	if (nverts == 4) {
+		if (p_quad_split_direction(co)) {
+			p_face_add(chart, key, vkeys, co, uv, 0, 1, 2, pin, select);
+			p_face_add(chart, key, vkeys, co, uv, 0, 2, 3, pin, select);
+		}
+		else {
+			p_face_add(chart, key, vkeys, co, uv, 0, 1, 3, pin, select);
+			p_face_add(chart, key, vkeys, co, uv, 1, 2, 3, pin, select);
+		}
+	}
+	else
+		p_face_add(chart, key, vkeys, co, uv, 0, 1, 2, pin, select);
+}
+
+void param_edge_set_seam(ParamHandle *handle, ParamKey *vkeys)
+{
+	PHandle *phandle = (PHandle*)handle;
+	PChart *chart = phandle->construction_chart;
+	PEdge *e;
+
+	param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
+
+	e = p_edge_lookup(chart, 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(chart, impl);
+	phandle->charts = p_split_charts(phandle, chart, phandle->ncharts);
+
+	p_chart_delete(chart);
+	phandle->construction_chart = NULL;
+
+	for (i = j = 0; i < phandle->ncharts; i++) {
+		p_chart_boundaries(phandle->charts[i], &nboundaries, &outer);
+
+		if (nboundaries > 0) {
+			phandle->charts[j] = phandle->charts[i];
+			j++;
+
+			if (fill && (nboundaries > 1))
+				p_chart_fill_boundaries(phandle->charts[i], outer);
+		}
+		else
+			p_chart_delete(phandle->charts[i]);
+	}
+
+	phandle->ncharts = j;
+
+	phandle->state = PHANDLE_STATE_CONSTRUCTED;
+}
+
+/* Least Squares Conformal Maps */
+
+static void p_chart_lscm_load_solution(PChart *chart)
+{
+	PVert *pin = chart->u.lscm.singlepin, *v;
+	float translation[2];
+
+	if (pin) {
+		translation[0] = pin->uv[0] - nlGetVariable(2*pin->u.index);
+		translation[1] = pin->uv[1] - nlGetVariable(2*pin->u.index + 1);
+	}
+	else
+		translation[0] = translation[1] = 0.0f;
+
+	for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
+		v->uv[0] = nlGetVariable(2*v->u.index) + translation[0];
+		v->uv[1] = nlGetVariable(2*v->u.index + 1) + translation[1];
+	}
+}
+
+static void p_chart_lscm_begin(PChart *chart)
+{
+	PVert *v, *pin1, *pin2;
+	int npins = 0, id = 0;
+
+	nlNewContext();
+	nlSolverParameteri(NL_NB_VARIABLES, 2*phash_size(chart->verts));
+	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
+
+	/* give vertices matrix indices and count pins */
+	for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
+		p_vert_load_pin_uvs(v);
+
+		if (v->flag & PVERT_PIN) {
+			npins++;
+			chart->u.lscm.singlepin = v;
+		}
+
+		v->u.index = id++;
+	}
+
+	if (npins <= 1) {
+		/* not enough pins, lets find some ourself */
+		p_extrema_verts(chart, &pin1, &pin2);
+
+		chart->u.lscm.pin1 = pin1;
+		chart->u.lscm.pin2 = pin2;
+	}
+	else
+		chart->u.lscm.singlepin = NULL;
+
+	chart->u.lscm.context = nlGetCurrent();
+}
+
+static PBool p_chart_lscm_solve(PChart *chart)
+{
+	PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
+	PFace *f;
+
+	nlMakeCurrent(chart->u.lscm.context);
+
+	nlBegin(NL_SYSTEM);
+
+	if (chart->u.lscm.pin1) {
+		nlLockVariable(2*pin1->u.index);
+		nlLockVariable(2*pin1->u.index + 1);
+		nlLockVariable(2*pin2->u.index);
+		nlLockVariable(2*pin2->u.index + 1);
+	
+		nlSetVariable(2*pin1->u.index, 0.0f);
+		nlSetVariable(2*pin1->u.index + 1, 0.5f);
+		nlSetVariable(2*pin2->u.index, 1.0f);
+		nlSetVariable(2*pin2->u.index + 1, 0.5f);
+	}
+	else {
+		/* set and lock the pins */
+		for (v=(PVert*)chart->verts->first; v; v=v->link.next) {
+			if (v->flag & PVERT_PIN) {
+				nlLockVariable(2*v->u.index);
+				nlLockVariable(2*v->u.index + 1);
+
+				nlSetVariable(2*v->u.index, v->uv[0]);
+				nlSetVariable(2*v->u.index + 1, v->uv[1]);
+			}
+		}
+	}
+
+	/* construct matrix */
+
+	nlBegin(NL_MATRIX);
+
+	for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
+		PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+		PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+		float a1, a2, a3, ratio, cosine, sine;
+		float sina1, sina2, sina3, sinmax;
+
+		if (chart->u.lscm.abf_alpha) {
+			/* use abf angles if passed on */
+			a1 = *(chart->u.lscm.abf_alpha++);
+			a2 = *(chart->u.lscm.abf_alpha++);
+			a3 = *(chart->u.lscm.abf_alpha++);
+		}
+		else
+			p_face_angles(f, &a1, &a2, &a3);
+
+		sina1 = sin(a1);
+		sina2 = sin(a2);
+		sina3 = sin(a3);
+
+		sinmax = MAX3(sina1, sina2, sina3);
+
+		/* shift vertices to find most stable order */
+		#define SHIFT3(type, a, b, c) \
+			{ type tmp; tmp = a; a = c; c = b; b = tmp; }
+
+		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 = sina2/sina3;
+		cosine = cos(a1)*ratio;
+		sine = sina1*ratio;
+
+		nlBegin(NL_ROW);
+		nlCoefficient(2*v1->u.index,   cosine - 1.0);
+		nlCoefficient(2*v1->u.index+1, -sine);
+		nlCoefficient(2*v2->u.index,   -cosine);
+		nlCoefficient(2*v2->u.index+1, sine);
+		nlCoefficient(2*v3->u.index,   1.0);
+		nlEnd(NL_ROW);
+
+		nlBegin(NL_ROW);
+		nlCoefficient(2*v1->u.index,   sine);
+		nlCoefficient(2*v1->u.index+1, cosine - 1.0);
+		nlCoefficient(2*v2->u.index,   -sine);
+		nlCoefficient(2*v2->u.index+1, -cosine);
+		nlCoefficient(2*v3->u.index+1, 1.0);
+		nlEnd(NL_ROW);
+	}
+
+	nlEnd(NL_MATRIX);
+
+	nlEnd(NL_SYSTEM);
+
+	if (nlSolveAdvanced(NULL, NL_TRUE)) {
+		p_chart_lscm_load_solution(chart);
+		p_flush_uvs(chart);
+
+		return P_TRUE;
+	}
+
+	return P_FALSE;
+}
+
+static void p_chart_lscm_end(PChart *chart)
+{
+	if (chart->u.lscm.context)
+		nlDeleteContext(chart->u.lscm.context);
+
+	chart->u.lscm.context = NULL;
+	chart->u.lscm.singlepin = NULL;
+	chart->u.lscm.pin1 = NULL;
+	chart->u.lscm.pin2 = NULL;
+}
+
+void param_lscm_begin(ParamHandle *handle)
+{
+	PHandle *phandle = (PHandle*)handle;
+	int i;
+
+	param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
+	phandle->state = PHANDLE_STATE_LSCM;
+
+	for (i = 0; i < phandle->ncharts; i++)
+		p_chart_lscm_begin(phandle->charts[i]);
+}
+
+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(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]);
+
+	phandle->state = PHANDLE_STATE_CONSTRUCTED;
+}
+
+/* Stretch */
+
+#define P_STRETCH_ITER 20
+
+static void p_stretch_pin_boundary(PChart *chart)
+{
+	PVert *v;
+
+	for(v=(PVert*)chart->verts->first; v; v=v->link.next)
+		if (v->edge->pair == NULL)
+			v->flag |= PVERT_PIN;
+		else
+			v->flag &= ~PVERT_PIN;
+}
+
+static float p_face_stretch(PFace *f)
+{
+	float T, w, tmp[3];
+	float Ps[3], Pt[3];
+	float a, c, area;
+	PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
+	PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
+
+	area = p_face_uv_area_signed(f);
+
+	if (area <= 0.0f) /* flipped face -> infinite stretch */
+		return 1e10f;
+	
+	if (f->flag & PFACE_FILLED)
+		return 0.0f;
+
+	w= 1.0f/(2.0f*area);
+
+	/* compute derivatives */
+	VecCopyf(Ps, v1->co);
+	VecMulf(Ps, (v2->uv[1] - v3->uv[1]));
+
+	VecCopyf(tmp, v2->co);
+	VecMulf(tmp, (v3->uv[1] - v1->uv[1]));
+	VecAddf(Ps, Ps, tmp);
+
+	VecCopyf(tmp, v3->co);
+	VecMulf(tmp, (v1->uv[1] - v2->uv[1]));
+	VecAddf(Ps, Ps, tmp);
+
+	VecMulf(Ps, w);
+
+	VecCopyf(Pt, v1->co);
+	VecMulf(Pt, (v3->uv[0] - v2->uv[0]));
+
+	VecCopyf(tmp, v2->co);
+	VecMulf(tmp, (v1->uv[0] - v3->uv[0]));
+	VecAddf(Pt, Pt, tmp);
+
+	VecCopyf(tmp, v3->co);
+	VecMulf(tmp, (v2->uv[0] - v1->uv[0]));
+	VecAddf(Pt, Pt, tmp);
+
+	VecMulf(Pt, w);
+
+	/* Sander Tensor */
+	a= Inpf(Ps, Ps);
+	c= Inpf(Pt, Pt);
+
+	T = sqrt(0.5f*(a + c)*f->u.area3d);
+
+	return T;
+}
+
+static float p_stretch_compute_vertex(PVert *v)
+{
+	PEdge *e = v->edge;
+	float sum = 0.0f;
+
+	do {
+		sum += p_face_stretch(e->face);
+		e = p_wheel_edge_next(e);
+	} while (e && e != (v->edge));
+
+	return sum;
+}
+
+static void p_chart_stretch_minimize(PChart *chart, RNG *rng)
+{
+	PVert *v;
+	PEdge *e;
+	int j, nedges;
+	float orig_stretch, low, stretch_low, high, stretch_high, mid, stretch;
+	float orig_uv[2], dir[2], random_angle, trusted_radius;
+
+	for(v=(PVert*)chart->verts->first; v; v=v->link.next) {
+		if((v->flag & PVERT_PIN) || !(v->flag & PVERT_SELECT))
+			continue;
+
+		orig_stretch = p_stretch_compute_vertex(v);
+		orig_uv[0] = v->uv[0];
+		orig_uv[1] = v->uv[1];
+
+		/* move vertex in a random direction */
+		trusted_radius = 0.0f;
+		nedges = 0;
+		e = v->edge;
+
+		do {
+			trusted_radius += p_edge_uv_length(e);
+			nedges++;
+
+			e = p_wheel_edge_next(e);
+		} while (e && e != (v->edge));
+
+		trusted_radius /= 2 * nedges;
+
+		random_angle = rng_getFloat(rng) * 2.0 * M_PI;
+		dir[0] = trusted_radius * cos(random_angle);
+		dir[1] = trusted_radius * sin(random_angle);
+
+		/* calculate old and new stretch */
+		low = 0;
+		stretch_low = orig_stretch;
+
+		Vec2Addf(v->uv, orig_uv, dir);
+		high = 1;
+		stretch = stretch_high = p_stretch_compute_vertex(v);
+
+		/* binary search for lowest stretch position */
+		for (j = 0; j < P_STRETCH_ITER; j++) {
+			mid = 0.5 * (low + high);
+			v->uv[0]= orig_uv[0] + mid*dir[0];
+			v->uv[1]= orig_uv[1] + mid*dir[1];
+			stretch = p_stretch_compute_vertex(v);
+
+			if (stretch_low < stretch_high) {
+				high = mid;
+				stretch_high = stretch;
+			}
+			else {
+				low = mid;
+				stretch_low = stretch;
+			}
+		}
+
+		/* no luck, stretch has increased, reset to old values */
+		if(stretch >= orig_stretch)
+			Vec2Copyf(v->uv, orig_uv);
+	}
+}
+
+void param_stretch_begin(ParamHandle *handle)
+{
+	PHandle *phandle = (PHandle*)handle;
+	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);
+
+	for (i = 0; i < phandle->ncharts; i++) {
+		chart = phandle->charts[i];
+
+		for (v=(PVert*)chart->verts->first; v; v=v->link.next)
+			p_vert_load_select_uvs(v);
+
+		p_stretch_pin_boundary(chart);
+
+		for (f=(PFace*)chart->faces->first; f; f=f->link.next) {
+			p_face_backup_uvs(f);
+			f->u.area3d = p_face_area(f);
+		}
+	}
+}
+
+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);
+		p_flush_uvs(chart);
+	}
+}
+
+void param_stretch_end(ParamHandle *handle, ParamBool restore)
+{
+	PHandle *phandle = (PHandle*)handle;
+	PChart *chart;
+	PFace *f;
+	int i;
+
+	param_assert(phandle->state == PHANDLE_STATE_STRETCH);
+	phandle->state = PHANDLE_STATE_CONSTRUCTED;
+
+	rng_free(phandle->rng);
+	phandle->rng = NULL;
+
+	if (restore) {
+		for (i = 0; i < phandle->ncharts; i++) {
+			chart = phandle->charts[i];
+
+			for (f=(PFace*)chart->faces->first; f; f=f->link.next)
+				p_face_restore_uvs(f);
+		}
+	}
+}
+
+/* Packing */
+
+static PBool p_pack_try(PHandle *handle, float side)
+{
+	PChart *chart;
+	float packx, packy, rowh, groupw, w, h;
+	int i;
+
+	packx= packy= 0.0;
+	rowh= 0.0;
+	groupw= 1.0/sqrt(handle->ncharts);
+
+	for (i = 0; i < handle->ncharts; i++) {
+		chart = handle->charts[i];
+		w = chart->u.pack.size[0];
+		h = chart->u.pack.size[1];
+
+		if(w <= (1.0-packx)) {
+			chart->u.pack.trans[0] = packx;
+			chart->u.pack.trans[1] = packy;
+
+			packx += w;
+			rowh= MAX2(rowh, h);
+		}
+		else {
+			packy += rowh;
+			packx = w;
+			rowh = h;
+
+			chart->u.pack.trans[0] = 0.0;
+			chart->u.pack.trans[1] = packy;
+		}
+
+		if (rowh > side)
+			return P_FALSE;
+	}
+
+	return P_TRUE;
+}
+
+#define PACK_SEARCH_DEPTH 7
+
+void param_pack(ParamHandle *handle)
+{
+	PHandle *phandle = (PHandle*)handle;
+	PChart *chart;
+	float uv_area, area, trans[2], minside, maxside, totarea, side;
+	int i;
+
+	/* very simple rectangle packing */
+
+	if (phandle->ncharts == 0)
+		return;
+
+	totarea = 0.0f;
+	maxside = 0.0f;
+
+	for (i = 0; i < phandle->ncharts; i++) {
+		chart = phandle->charts[i];
+
+		p_chart_area(chart, &uv_area, &area);
+		chart->u.pack.rescale = (uv_area > 0.0f)? area/uv_area: 0.0f;
+		totarea += uv_area*chart->u.pack.rescale;
+
+		p_chart_uv_bbox(chart, trans, chart->u.pack.size);
+		trans[0] = -trans[0];
+		trans[1] = -trans[1];
+		p_chart_uv_translate(chart, trans);
+		p_chart_uv_scale(chart, chart->u.pack.rescale);
+
+		chart->u.pack.size[0] -= trans[0];
+		chart->u.pack.size[1] -= trans[1];
+		chart->u.pack.size[0] *= chart->u.pack.rescale;
+		chart->u.pack.size[1] *= chart->u.pack.rescale;
+
+		maxside = MAX3(maxside, chart->u.pack.size[0], chart->u.pack.size[1]);
+	}
+
+	return;
+
+	printf("%f\n", maxside);
+
+	/* binary search over pack region size */
+	minside = sqrt(totarea);
+	maxside = (((int)sqrt(phandle->ncharts-1))+1)*maxside;
+
+	for (i = 0; i < PACK_SEARCH_DEPTH; i++) {
+		printf("%f %f\n", minside, maxside);
+		if (p_pack_try(phandle, (minside+maxside)*0.5f))
+			minside = (minside+maxside)*0.5f;
+		else
+			maxside = (minside+maxside)*0.5f;
+	}
+
+	side = maxside + 1e-5; /* prevent floating point errors */
+	if (!p_pack_try(phandle, side))
+		printf("impossible\n");
+
+	for (i = 0; i < phandle->ncharts; i++) {
+		chart = phandle->charts[i];
+
+		p_chart_uv_scale(chart, chart->u.pack.rescale/side);
+		trans[0] = chart->u.pack.trans[0]/side;
+		trans[1] = chart->u.pack.trans[1]/side;
+		p_chart_uv_translate(chart, trans);
+	}
+}
+
diff --git a/source/blender/src/parametrizer.h b/source/blender/src/parametrizer.h
new file mode 100644
index 00000000000..a3904af83b9
--- /dev/null
+++ b/source/blender/src/parametrizer.h
@@ -0,0 +1,73 @@
+
+#ifndef __PARAMETRIZER_H__
+#define __PARAMETRIZER_H__
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef void ParamHandle;	/* handle to a set of charts */
+typedef long ParamKey;		/* (hash) key for identifying verts and faces */
+typedef enum ParamBool {
+	PARAM_TRUE = 1,
+	PARAM_FALSE = 0
+} ParamBool;
+
+/* Chart construction:
+   -------------------
+   - faces and seams may only be added between construct_{begin|end}
+   - the pointers to co and uv are stored, rather than being copied
+   - vertices are implicitly created
+   - in construct_end the mesh will be split up according to the seams
+   - the resulting charts must be:
+      - manifold, connected, open (at least one boundary loop)
+   - output will be written to the uv pointers
+*/
+
+ParamHandle *param_construct_begin();
+
+void param_face_add(ParamHandle *handle,
+                    ParamKey key,
+                    int nverts,	
+                    ParamKey *vkeys,
+                    float **co,
+                    float **uv,
+					ParamBool *pin,
+					ParamBool *select);
+
+void param_edge_set_seam(ParamHandle *handle,
+                         ParamKey *vkeys);
+
+void param_construct_end(ParamHandle *handle, ParamBool fill, ParamBool impl);
+void param_delete(ParamHandle *chart);
+
+/* Least Squares Conformal Maps:
+   -----------------------------
+   - charts with less than two pinned vertices are assigned 2 pins
+   - lscm is divided in three steps:
+      - begin: compute matrix and it's factorization (expensive)
+      - solve using pinned coordinates (cheap)
+	  - end: clean up 
+	- uv coordinates are allowed to change within begin/end, for
+	  quick re-solving
+*/
+
+void param_lscm_begin(ParamHandle *handle);
+void param_lscm_solve(ParamHandle *handle);
+void param_lscm_end(ParamHandle *handle);
+
+/* Stretch */
+
+void param_stretch_begin(ParamHandle *handle);
+void param_stretch_iter(ParamHandle *handle);
+void param_stretch_end(ParamHandle *handle, ParamBool restore);
+
+/* Packing */
+void param_pack(ParamHandle *handle);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /*__PARAMETRIZER_H__*/
+
diff --git a/source/blender/src/parametrizer_intern.h b/source/blender/src/parametrizer_intern.h
new file mode 100644
index 00000000000..4086372e2d5
--- /dev/null
+++ b/source/blender/src/parametrizer_intern.h
@@ -0,0 +1,186 @@
+
+#ifndef __PARAMETRIZER_INTERN_H__
+#define __PARAMETRIZER_INTERN_H__
+
+/* Hash:
+   -----
+   - insert only
+   - elements are all stored in a flat linked list
+*/
+
+typedef long PHashKey;
+
+typedef struct PHashLink {
+	struct PHashLink *next;
+	PHashKey key;
+} PHashLink;
+
+typedef struct PHash {
+	PHashLink *first;
+	PHashLink **buckets;
+	int size, cursize, cursize_id;
+} PHash;
+
+PHash *phash_new(int sizehint);
+void phash_delete_with_links(PHash *ph);
+void phash_delete(PHash *ph);
+
+int phash_size(PHash *ph);
+
+void phash_insert(PHash *ph, PHashLink *link);
+PHashLink *phash_lookup(PHash *ph, PHashKey key);
+PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link);
+
+#if 0
+	#define param_assert(condition)
+	#define param_warning(message);
+#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);
+#endif
+
+typedef enum PBool {
+	P_TRUE = 1,
+	P_FALSE = 0
+} PBool;
+
+struct PVert;
+struct PEdge;
+struct PFace;
+struct PChart;
+struct PHandle;
+
+/* Heap */
+
+typedef struct PHeapLink {
+	void *ptr;
+	float value;
+	int index;
+} PHeapLink;
+
+typedef struct PHeap {
+	unsigned int size;
+	unsigned int bufsize;
+	PHeapLink *links;
+	PHeapLink **tree;
+} PHeap;
+
+/* Simplices */
+
+typedef struct PVert {
+	struct PVertLink {
+		struct PVert *next;
+		PHashKey key;
+	} link;
+
+	struct PEdge *edge;
+	float *co;
+	float uv[2];
+	int flag;
+
+	union PVertUnion {
+		int index; /* lscm matrix index */
+		float distortion; /* area smoothing */
+	} u;
+} PVert; 
+
+typedef struct PEdge {
+	struct PEdgeLink {
+		struct PEdge *next;
+		PHashKey key;
+	} link;
+
+	struct PVert *vert;
+	struct PEdge *pair;
+	struct PEdge *next;
+	struct PFace *face;
+	float *orig_uv, old_uv[2];
+	int flag;
+
+	union PEdgeUnion {
+		PHeapLink *heaplink;
+	} u;
+} PEdge;
+
+typedef struct PFace {
+	struct PFaceLink {
+		struct PFace *next;
+		PHashKey key;
+	} link;
+
+	struct PEdge *edge;
+	int flag;
+
+	union PFaceUnion {
+		int chart; /* chart construction */
+		float area3d; /* stretch */
+	} u;
+} PFace;
+
+enum PVertFlag {
+	PVERT_PIN = 1,
+	PVERT_SELECT = 2
+};
+
+enum PEdgeFlag {
+	PEDGE_SEAM = 1,
+	PEDGE_VERTEX_SPLIT = 2,
+	PEDGE_PIN = 4,
+	PEDGE_SELECT = 8,
+	PEDGE_DONE = 16,
+	PEDGE_FILLED = 32
+};
+
+/* for flipping faces */
+#define PEDGE_VERTEX_FLAGS (PEDGE_PIN)
+
+enum PFaceFlag {
+	PFACE_CONNECTED = 1,
+	PFACE_FILLED = 2
+};
+
+/* Chart */
+
+typedef struct PChart {
+	PHash *verts;
+	PHash *edges;
+	PHash *faces;
+
+	union PChartUnion {
+		struct PChartLscm {
+			NLContext context;
+			float *abf_alpha;
+			PVert *singlepin;
+			PVert *pin1, *pin2;
+		} lscm;
+		struct PChartPack {
+			float rescale;
+			float size[2], trans[2];
+		} pack;
+	} u;
+
+	struct PHandle *handle;
+} PChart;
+
+enum PHandleState {
+	PHANDLE_STATE_ALLOCATED,
+	PHANDLE_STATE_CONSTRUCTED,
+	PHANDLE_STATE_LSCM,
+	PHANDLE_STATE_STRETCH,
+};
+
+typedef struct PHandle {
+	PChart *construction_chart;
+	PChart **charts;
+	int ncharts;
+	enum PHandleState state;
+	MemArena *arena;
+	PBool implicit;
+	RNG *rng;
+} PHandle;
+
+#endif /*__PARAMETRIZER_INTERN_H__*/
+
diff --git a/source/blender/src/space.c b/source/blender/src/space.c
index 47c3b8e4004..73552449fba 100644
--- a/source/blender/src/space.c
+++ b/source/blender/src/space.c
@@ -3867,8 +3867,7 @@ static void winqreadimagespace(ScrArea *sa, void *spacedata, BWinEvent *evt)
 					sample_vpaint();
 				break;
 			case AKEY:
-				if((G.qual==0))
-					select_swap_tface_uv();
+				select_swap_tface_uv();
 				break;
 			case BKEY:
 				if(G.qual==LR_SHIFTKEY)
@@ -3956,6 +3955,8 @@ static void winqreadimagespace(ScrArea *sa, void *spacedata, BWinEvent *evt)
 					stitch_uv_tface(0);
 				else if(G.qual==LR_SHIFTKEY)
 					stitch_uv_tface(1);
+				else if(G.qual==LR_CTRLKEY)
+					minimize_stretch_tface_uv();
 				break;
 			case WKEY:
 				weld_align_menu_tface_uv();
diff --git a/source/blender/src/unwrapper.c b/source/blender/src/unwrapper.c
index 045eea31625..7b61e8f51f0 100644
--- a/source/blender/src/unwrapper.c
+++ b/source/blender/src/unwrapper.c
@@ -43,6 +43,8 @@
 #include "DNA_mesh_types.h"
 #include "DNA_meshdata_types.h"
 #include "DNA_scene_types.h"
+#include "DNA_screen_types.h"
+#include "DNA_space_types.h"
 
 #include "BKE_global.h"
 #include "BKE_mesh.h"
@@ -53,6 +55,7 @@
 
 #include "BIF_editsima.h"
 #include "BIF_space.h"
+#include "BIF_screen.h"
 
 #include "blendef.h"
 #include "mydevice.h"
@@ -60,6 +63,10 @@
 #include "ONL_opennl.h"
 #include "BDR_unwrapper.h"
 
+#include "PIL_time.h"
+
+#include "parametrizer.h"
+
 /* Implementation Least Squares Conformal Maps parameterization, based on
  * chapter 2 of:
  * Bruno Levy, Sylvain Petitjean, Nicolas Ray, Jerome Maillot. Least Squares
@@ -903,9 +910,9 @@ static int unwrap_lscm_face_group(Mesh *me, int *groups, int gid)
 		lscm_build_matrix(me, lscmvert, groups, gid, center, radius);
 		
 		nlEnd(NL_SYSTEM);
-		
+
 		/* LSCM solver magic! */
-		nlSolve();
+		nlSolve(NULL, NL_FALSE);
 		
 		/* load new uv's: will be projected uv's if solving failed  */
 		lscm_load_solution(me, lscmvert, groups, gid);
@@ -1336,3 +1343,162 @@ void select_linked_tfaces_with_seams(int mode, Mesh *me, unsigned int index)
 	object_tface_flags_changed(OBACT, 0);
 }
 
+/* Parametrizer */
+
+ParamHandle *construct_param_handle(Mesh *me, short implicit, short fill)
+{
+	int a;
+	TFace *tf;
+	MFace *mf;
+	MVert *mv;
+	MEdge *medge;
+	ParamHandle *handle;
+	
+	handle = param_construct_begin();
+	
+	mv= me->mvert;
+	mf= me->mface;
+	tf= me->tface;
+	for (a=0; a<me->totface; a++, mf++, tf++) {
+		ParamKey key, vkeys[4];
+		ParamBool pin[4], select[4];
+		float *co[4];
+		float *uv[4];
+		int nverts;
+
+		if ((tf->flag & TF_HIDE) || !(tf->flag & TF_SELECT))
+			continue;
+
+		if (implicit && !(tf->flag & (TF_SEL1|TF_SEL2|TF_SEL3|TF_SEL4)))
+			continue;
+
+		key = (ParamKey)mf;
+		vkeys[0] = (ParamKey)mf->v1;
+		vkeys[1] = (ParamKey)mf->v2;
+		vkeys[2] = (ParamKey)mf->v3;
+
+		co[0] = (mv+mf->v1)->co;
+		co[1] = (mv+mf->v2)->co;
+		co[2] = (mv+mf->v3)->co;
+
+		uv[0] = tf->uv[0];
+		uv[1] = tf->uv[1];
+		uv[2] = tf->uv[2];
+
+		pin[0] = ((tf->flag & TF_PIN1) != 0);
+		pin[1] = ((tf->flag & TF_PIN2) != 0);
+		pin[2] = ((tf->flag & TF_PIN3) != 0);
+
+		select[0] = ((tf->flag & TF_SEL1) != 0);
+		select[1] = ((tf->flag & TF_SEL2) != 0);
+		select[2] = ((tf->flag & TF_SEL3) != 0);
+
+		if (mf->v4) {
+			vkeys[3] = (ParamKey)mf->v4;
+			co[3] = (mv+mf->v4)->co;
+			uv[3] = tf->uv[3];
+			pin[3] = ((tf->flag & TF_PIN4) != 0);
+			select[3] = ((tf->flag & TF_SEL4) != 0);
+			nverts = 4;
+		}
+		else
+			nverts = 3;
+
+		param_face_add(handle, key, nverts, vkeys, co, uv, pin, select);
+	}
+
+	if (!implicit) {
+		for(medge=me->medge, a=me->totedge; a>0; a--, medge++) {
+			if(medge->flag & ME_SEAM) {
+				ParamKey vkeys[2];
+
+				vkeys[0] = (ParamKey)medge->v1;
+				vkeys[1] = (ParamKey)medge->v2;
+				param_edge_set_seam(handle, vkeys);
+			}
+		}
+	}
+
+	param_construct_end(handle, fill, implicit);
+
+	return handle;
+}
+
+#if 0
+void unwrap_lscm(void)
+{
+	Mesh *me;
+	ParamHandle *handle;
+	
+	me= get_mesh(OBACT);
+	if(me==0 || me->tface==0) return;
+
+	handle = construct_param_handle(me, 0, 1);
+
+	param_lscm_begin(handle);
+	param_lscm_solve(handle);
+	param_lscm_end(handle);
+
+	param_pack(handle);
+
+	param_delete(handle);
+
+	BIF_undo_push("UV lscm unwrap");
+
+	object_uvs_changed(OBACT);
+
+	allqueue(REDRAWVIEW3D, 0);
+	allqueue(REDRAWIMAGE, 0);
+}
+#endif
+
+void minimize_stretch_tface_uv(void)
+{
+	Mesh *me;
+	ParamHandle *handle;
+	double lasttime;
+	short doit = 1, val;
+	unsigned short event = 0;
+	
+	me = get_mesh(OBACT);
+	if(me==0 || me->tface==0) return;
+
+	handle = construct_param_handle(me, 1, 0);
+
+	lasttime = PIL_check_seconds_timer();
+
+	param_stretch_begin(handle);
+
+	while (doit) {
+		param_stretch_iter(handle);
+
+		while (qtest()) {
+			event= extern_qread(&val);
+			if (val && (event==ESCKEY || event==RETKEY || event==PADENTER))
+				doit = 0;
+		}
+		
+		if (!doit)
+			break;
+
+		if (PIL_check_seconds_timer() - lasttime > 0.5) {
+			headerprint("Enter to finish. Escape to cancel.");
+
+			lasttime = PIL_check_seconds_timer();
+			if(G.sima->lock) force_draw_plus(SPACE_VIEW3D, 0);
+			else force_draw(0);
+		}
+	}
+
+	param_stretch_end(handle, event==ESCKEY);
+
+	param_delete(handle);
+
+	BIF_undo_push("UV stretch minimize");
+
+	object_uvs_changed(OBACT);
+
+	allqueue(REDRAWVIEW3D, 0);
+	allqueue(REDRAWIMAGE, 0);
+}
+