Index: trunk/psLib/src/math/psMatrix.c =================================================================== --- trunk/psLib/src/math/psMatrix.c (revision 770) +++ trunk/psLib/src/math/psMatrix.c (revision 799) @@ -20,6 +20,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-05-24 23:30:52 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-05-28 02:53:18 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -68,98 +68,79 @@ /*****************************************************************************/ -// None +/** Preprocessor macro to generate error a NULL image */ +#define PS_CHECK_NULL_VECTOR(NAME, RETURN) \ +if (NAME == NULL || NAME->vec.v == NULL) { \ + psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \ + return RETURN; \ +} + +/** Preprocessor macro to create vector based on another */ +#define PS_CHECK_ALLOC_VECTOR(NAME, SIZE, PS_TYPE) \ +if(NAME == NULL) { \ + NAME = psVectorAlloc(PS_TYPE, SIZE); \ +} + +/** Preprocessor macro to generate error for zero length vector */ +#define PS_CHECK_SIZE_VECTOR(NAME, RETURN) \ +if (NAME->n < 1) { \ + psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \ + return RETURN; \ +} + +/** Preprocessor macro to generate error a NULL image */ +#define PS_CHECK_NULL_IMAGE(NAME, RETURN) \ +if (NAME == NULL || NAME->data.v == NULL) { \ + psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \ + return RETURN; \ +} + +/** Preprocessor macro to create image based on another */ +#define PS_CHECK_ALLOC_IMAGE(NAME, NCOLS, NROWS, PS_TYPE) \ +if(NAME == NULL) { \ + NAME = psImageAlloc(NCOLS, NROWS, PS_TYPE); \ +} + +/** Preprocessor macro to generate error for zero length rows or columns */ +#define PS_CHECK_SIZE_IMAGE(NAME, RETURN) \ +if (NAME->numCols < 1 || NAME->numRows < 1) { \ + psError(__func__,"Invalid operation: %s has zero rows or columns (%dx%d).", #NAME, \ + NAME->numCols, NAME->numRows); \ + return RETURN; \ +} + +/** Preprocessor macro to generate error for image dimensionality not set to PS_DIMEN_IMAGE */ +#define PS_CHECK_DIMEN_AND_TYPE(NAME, PS_DIMEN, RETURN) \ +if (NAME->type.dimen != PS_DIMEN) { \ + psError(__func__,"Invalid operation: %s incorrect dimensionality %d.", #NAME, PS_DIMEN); \ + return RETURN; \ +} else if(NAME->type.type != PS_TYPE_F64) { \ + psError(__func__, "Invalid operation: %s not PS_TYPE_F64.", #NAME); \ + return RETURN; \ +} + +/** Preprocessor macro to check that input is not equal to output */ +#define PS_CHECK_POINTERS(NAME1, NAME2, RETURN) \ +if (NAME1 == NAME2) { \ + psError(__func__,"Invalid operation: Pointer to %s is same as %s.", #NAME1, #NAME2); \ + return RETURN; \ +} + +/** Preprocessor macro to check that an image is square */ +#define PS_CHECK_SQUARE(NAME, RETURN) \ +if (NAME->numCols != NAME->numRows) { \ + psError(__func__,"Invalid operation: %s not square array.", #NAME); \ + return RETURN; \ +} + +/** Preprocessor macro to initalize a GSL matrix. */ +#define PS_GSL_MATRIX_INITIALIZE(LHS_NAME, RHS_NAME) \ +LHS_NAME.size1 = numRows; \ +LHS_NAME.size2 = numCols; \ +LHS_NAME.tda = numCols; \ +LHS_NAME.data = RHS_NAME; /*****************************************************************************/ /* FUNCTION IMPLEMENTATION - PUBLIC */ /*****************************************************************************/ - -#define PS_MATRIX_TRANSPOSE(PS_TYPE) \ -{ \ - int arraySize = 0; \ - int numRows = 0; \ - int numCols = 0; \ - gsl_matrix##PS_TYPE trans; \ - \ - /* Initialize data */ \ - numRows = inImage->numRows; \ - numCols = inImage->numCols; \ - arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols; \ - \ - /* Copy psImage input data into psImage output data to keep input data pristine */ \ - memcpy(outImage->data.v[0], inImage->data.v[0], arraySize); \ - \ - /* Manually fill inverted output matrix so it will be aligned with output image */ \ - trans.size1 = numRows; \ - trans.size2 = numCols; \ - trans.tda = numCols; \ - trans.data = outImage->data.v[0]; \ - \ - /* Transpose data */ \ - gsl_matrix##PS_TYPE##_transpose(&trans); \ -} - -psImage* psMatrixTranspose(psImage *outImage, const psImage *inImage) -{ - psElemType elemType = 0; - - elemType = inImage->type.type; - switch(elemType) { - case PS_TYPE_FLOAT: - PS_MATRIX_TRANSPOSE(_float); - break; - case PS_TYPE_DOUBLE: - PS_MATRIX_TRANSPOSE(); - break; - default: - psError(__func__, " : Line %d - Invalid psElemType: %d\n", __LINE__, elemType); - } - - return outImage; -} - -psImage* psMatrixOp(psImage *outImage, psImage *inImage1, const char op, psImage *inImage2) -{ - int arraySize = 0; - int numRows = 0; - int numCols = 0; - gsl_matrix m2; - gsl_matrix m3; - - // Initialize data - numRows = inImage1->numRows; - numCols = inImage1->numCols; - arraySize = PSELEMTYPE_SIZEOF(outImage->type.type)*numRows*numCols; - - m2.size1 = numRows; - m2.size2 = numCols; - m2.tda = numCols; - m2.data = inImage2->data.v[0]; - - m3.size1 = numRows; - m3.size2 = numCols; - m3.tda = numCols; - m3.data = outImage->data.v[0]; - - // Copy psImage input data into GSL matrix data to keep input data pristine - memcpy(m3.data, inImage1->data.v[0], arraySize); - - - switch(op) { - case '+': - gsl_matrix_add(&m3, &m2); - break; - case '-': - gsl_matrix_sub(&m3, &m2); - break; - case '*': - gsl_linalg_matmult(&m3, &m2, &m3); - break; - default: - psError(__func__, " : Line %d - Invalid psMatrixOp operation: %c\n", __LINE__, op); - \ - } - - return NULL; -} psImage *psMatrixLUD(psImage *outImage, psVector *outPerm, psImage *inImage) @@ -171,4 +152,16 @@ gsl_matrix lu; gsl_permutation perm; + + // Error checks + PS_CHECK_POINTERS(inImage, outImage, outImage); + PS_CHECK_NULL_IMAGE(inImage, outImage); + PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(inImage, outImage); + PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type); + PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(outImage, outImage); + PS_CHECK_ALLOC_VECTOR(outPerm, inImage->numRows, inImage->type.type); + PS_CHECK_NULL_VECTOR(outPerm, outImage); + PS_CHECK_DIMEN_AND_TYPE(outPerm, PS_DIMEN_VECTOR, outImage); // Initialize data @@ -177,12 +170,13 @@ arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols; - // Manually fill GSL structs so they will be aligned with output data + // Initialize GSL data perm.size = numCols; outPerm->n = numCols; perm.data = outPerm->vec.v; - lu.size1 = numRows; - lu.size2 = numCols; - lu.tda = numCols; - lu.data = outImage->data.v[0]; + PS_GSL_MATRIX_INITIALIZE(lu, outImage->data.v[0]); + + // Non-square matrices not allowed + PS_CHECK_SQUARE(inImage, outImage); + PS_CHECK_SQUARE(outImage, outImage); // Copy psImage input data into GSL matrix data to keep input data pristine @@ -195,5 +189,5 @@ } -psVector *psMatrixLUSolve(psVector *outVector, const psImage *luImage, const psVector *inVector, const psVector *inPerm) +psVector *psMatrixLUSolve(psVector *outVector, const psImage *inImage, const psVector *inVector, const psVector *inPerm) { int arraySize = 0; @@ -205,14 +199,25 @@ gsl_vector x; - // Initialize data - numRows = luImage->numRows; - numCols = luImage->numCols; - arraySize = PSELEMTYPE_SIZEOF(luImage->type.type)*numRows*numCols; - - // Manually fill GSL structs so they will be aligned with output data - lu.size1 = luImage->numRows; - lu.size2 = luImage->numCols; - lu.tda = luImage->numCols; - lu.data = luImage->data.v[0]; + // Error checks + PS_CHECK_POINTERS(outVector, inVector, outVector); + PS_CHECK_POINTERS(inVector, inPerm, outVector); + PS_CHECK_POINTERS(outVector, inPerm, outVector); + PS_CHECK_NULL_IMAGE(inImage, outVector); + PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outVector); + PS_CHECK_SIZE_IMAGE(inImage, outVector); + PS_CHECK_NULL_VECTOR(outVector, outVector); + PS_CHECK_DIMEN_AND_TYPE(outVector, PS_DIMEN_VECTOR, outVector); + PS_CHECK_NULL_VECTOR(inVector, outVector); + PS_CHECK_DIMEN_AND_TYPE(inVector, PS_DIMEN_VECTOR, outVector); + PS_CHECK_NULL_VECTOR(inPerm, outVector); + PS_CHECK_DIMEN_AND_TYPE(inPerm, PS_DIMEN_VECTOR, outVector); + + // Initialize data + numRows = inImage->numRows; + numCols = inImage->numCols; + arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols; + + // Initialize GSL data + PS_GSL_MATRIX_INITIALIZE(lu, inImage->data.v[0]); outVector->n = numCols; @@ -245,4 +250,13 @@ gsl_permutation *perm = NULL; + // Error checks + PS_CHECK_POINTERS(inImage, outImage, outImage); + PS_CHECK_NULL_IMAGE(inImage, outImage); + PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(inImage, outImage); + PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type); + PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(outImage, outImage); + // Initialize data numRows = inImage->numRows; @@ -254,12 +268,13 @@ lu = gsl_matrix_alloc(numRows, numCols); + // Initialize GSL data + PS_GSL_MATRIX_INITIALIZE(inv, outImage->data.v[0]); + + // Non-square matrices not allowed + PS_CHECK_SQUARE(inImage, outImage); + PS_CHECK_SQUARE(outImage, outImage); + // Copy psImage input data into GSL matrix data to keep input data pristine memcpy(lu->data, inImage->data.v[0], arraySize); - - // Manually fill inverted output matrix so it will be aligned with output image - inv.size1 = outImage->numRows; - inv.size2 = outImage->numCols; - inv.tda = outImage->numCols; - inv.data = outImage->data.v[0]; // Invert data and calculate determinant @@ -285,4 +300,9 @@ gsl_permutation *perm = NULL; + // Error checks + PS_CHECK_NULL_IMAGE(inImage, 0); + PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, 0); + PS_CHECK_SIZE_IMAGE(inImage, 0); + // Initialize data numRows = inImage->numRows; @@ -294,4 +314,7 @@ lu = gsl_matrix_alloc(numRows, numCols); + // Non-square matrices not allowed + PS_CHECK_SQUARE(inImage, 0); + // Copy psImage input data into GSL matrix data to keep input data pristine memcpy(lu->data, inImage->data.v[0], arraySize); @@ -306,4 +329,83 @@ return det; +} + +psImage* psMatrixMultiply(psImage *outImage, psImage *inImage1, psImage *inImage2) +{ + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix m1; + gsl_matrix m2; + gsl_matrix m3; + + // Error checks + PS_CHECK_POINTERS(inImage1, outImage, outImage); + PS_CHECK_POINTERS(inImage1, inImage2, outImage); + PS_CHECK_NULL_IMAGE(inImage1, outImage); + PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(inImage1, outImage); + PS_CHECK_NULL_IMAGE(inImage2, outImage); + PS_CHECK_DIMEN_AND_TYPE(inImage2, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(inImage2, outImage); + PS_CHECK_ALLOC_IMAGE(outImage, inImage2->numCols, inImage2->numRows, inImage2->type.type); + PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(outImage, outImage); + + // Initialize data + numRows = inImage1->numRows; + numCols = inImage1->numCols; + arraySize = PSELEMTYPE_SIZEOF(outImage->type.type)*numRows*numCols; + + // Initialize GSL data + PS_GSL_MATRIX_INITIALIZE(m1, inImage1->data.v[0]); + PS_GSL_MATRIX_INITIALIZE(m2, inImage2->data.v[0]); + PS_GSL_MATRIX_INITIALIZE(m3, outImage->data.v[0]); + + // Non-square matrices not allowed + PS_CHECK_SQUARE(inImage1, outImage); + PS_CHECK_SQUARE(inImage2, outImage); + PS_CHECK_SQUARE(outImage, outImage); + + gsl_linalg_matmult(&m1, &m2, &m3); + + return outImage; +} + +psImage* psMatrixTranspose(psImage *outImage, const psImage *inImage) +{ + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix trans; + + // Error checks + PS_CHECK_POINTERS(inImage, outImage, outImage); + PS_CHECK_NULL_IMAGE(inImage, outImage); + PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(inImage, outImage); + PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type); + PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(outImage, outImage); + + // Initialize data + numRows = inImage->numRows; + numCols = inImage->numCols; + arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols; + + // Initialize GSL data + PS_GSL_MATRIX_INITIALIZE(trans, outImage->data.v[0]); + + // Non-square matrices not allowed + PS_CHECK_SQUARE(inImage, outImage); + PS_CHECK_SQUARE(outImage, outImage); + + // Copy psImage input data into psImage output data to keep input data pristine + memcpy(outImage->data.v[0], inImage->data.v[0], arraySize); + + // Transpose data + gsl_matrix_transpose(&trans); + + return outImage; } @@ -317,19 +419,20 @@ gsl_matrix in; + // Error checks + PS_CHECK_POINTERS(inImage, outImage, outImage); + PS_CHECK_NULL_IMAGE(inImage, outImage); + PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(inImage, outImage); + PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type); + PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage); + PS_CHECK_SIZE_IMAGE(outImage, outImage); + // Initialize data numRows = inImage->numRows; numCols = inImage->numCols; - out.data = outImage->data.v[0]; - - // Manually fill GSL structs so they will be aligned with output data - in.size1 = numRows; - in.size2 = numCols; - in.tda = numCols; - in.data = inImage->data.v[0]; - - out.size1 = numRows; - out.size2 = numCols; - out.tda = numCols; - out.data = outImage->data.v[0]; + + // Initialize GSL data + PS_GSL_MATRIX_INITIALIZE(in, inImage->data.v[0]); + PS_GSL_MATRIX_INITIALIZE(out, outImage->data.v[0]); // Allocate GSL structs @@ -337,5 +440,9 @@ w = gsl_eigen_symmv_alloc(numRows); - // Calculate Eigenvalues and Eigenvectors. Eigenvalues not currently used + // Non-square matrices not allowed + PS_CHECK_SQUARE(inImage, outImage); + PS_CHECK_SQUARE(outImage, outImage); + + // Calculate Eigenvalues and Eigenvectors...Eigenvalues not currently used gsl_eigen_symmv(&in, eVals, &out, w); @@ -349,8 +456,27 @@ psVector *psMatrixToVector(psVector *outVector, psImage *inImage) { - // if inimage isn't 1d - // if inimage nrows != vector len - int colSize = 0; + + // Error checks + PS_CHECK_NULL_IMAGE(inImage, outVector); + PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outVector); + PS_CHECK_SIZE_IMAGE(inImage, outVector); + PS_CHECK_ALLOC_VECTOR(outVector, inImage->numRows, inImage->type.type); + PS_CHECK_NULL_VECTOR(outVector, outVector); + PS_CHECK_DIMEN_AND_TYPE(outVector, PS_DIMEN_VECTOR, outVector); + + // Set n if allocated, but empty + if(outVector->n == 0) { + outVector->n = inImage->numRows; + } + + // More checks + if(inImage->numCols > 1) { + psError(__func__, "Image has more than 1 column: numCols = %d.", inImage->numCols); + return outVector; + } else if(outVector->n != inImage->numRows) { + psError(__func__, "Image and vector sizes differ: (%d vs %d).", inImage->numRows, outVector->n); + return outVector; + } colSize = PSELEMTYPE_SIZEOF(inImage->type.type)*inImage->numRows; @@ -362,8 +488,22 @@ psImage *psVectorToMatrix(psImage *outImage, psVector *inVector) { - // if inimage isn't 1d - // if inimage nrows != vector len - int colSize = 0; + + // Error checks + PS_CHECK_NULL_VECTOR(inVector, outImage); + PS_CHECK_DIMEN_AND_TYPE(inVector, PS_DIMEN_VECTOR, outImage); + PS_CHECK_SIZE_VECTOR(inVector, outImage); + PS_CHECK_ALLOC_IMAGE(outImage, 1, inVector->n, PS_TYPE_F32) + PS_CHECK_NULL_IMAGE(outImage, outImage); + PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage); + + // More checks + if(outImage->numCols > 1) { + psError(__func__, "Image has more than 1 column: numCols = %d.", outImage->numCols); + return outImage; + } else if(outImage->numRows != inVector->n) { + psError(__func__, "Image and vector sizes differ: (%d vs %d).", outImage->numRows, inVector->n); + return outImage; + } colSize = PSELEMTYPE_SIZEOF(outImage->type.type)*outImage->numRows;