diff options
| -rw-r--r-- | intern/opennl/extern/ONL_opennl.h | 8 | ||||
| -rw-r--r-- | intern/opennl/intern/opennl.c | 1631 | ||||
| -rw-r--r-- | source/blender/include/BDR_unwrapper.h | 1 | ||||
| -rw-r--r-- | source/blender/src/header_image.c | 4 | ||||
| -rw-r--r-- | source/blender/src/parametrizer.c | 1785 | ||||
| -rw-r--r-- | source/blender/src/parametrizer.h | 73 | ||||
| -rw-r--r-- | source/blender/src/parametrizer_intern.h | 186 | ||||
| -rw-r--r-- | source/blender/src/space.c | 5 | ||||
| -rw-r--r-- | source/blender/src/unwrapper.c | 170 |
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); +} + |