IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 799


Ignore:
Timestamp:
May 27, 2004, 4:53:27 PM (22 years ago)
Author:
harman
Message:

Added to psMatrix functionality

Location:
trunk/psLib/src
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/dataManip/psMatrix.c

    r770 r799  
    2020 *  @author Ross Harman, MHPCC
    2121 *   
    22  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    23  *  @date $Date: 2004-05-24 23:30:52 $
     22 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     23 *  @date $Date: 2004-05-28 02:53:18 $
    2424 *
    2525 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    6868/*****************************************************************************/
    6969
    70 // None
     70/** Preprocessor macro to generate error a NULL image */
     71#define PS_CHECK_NULL_VECTOR(NAME, RETURN)                                                                    \
     72if (NAME == NULL || NAME->vec.v == NULL) {                                                                \
     73    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                                \
     74    return RETURN;                                                                                        \
     75}
     76
     77/** Preprocessor macro to create vector based on another */
     78#define PS_CHECK_ALLOC_VECTOR(NAME, SIZE, PS_TYPE)                                                            \
     79if(NAME == NULL) {                                                                                        \
     80    NAME = psVectorAlloc(PS_TYPE, SIZE);                                                                  \
     81}
     82
     83/** Preprocessor macro to generate error for zero length vector */
     84#define PS_CHECK_SIZE_VECTOR(NAME, RETURN)                                                                    \
     85if (NAME->n < 1) {                                                                                        \
     86    psError(__func__,"Invalid operation: %s has zero n value.", #NAME);                                   \
     87    return RETURN;                                                                                        \
     88}
     89
     90/** Preprocessor macro to generate error a NULL image */
     91#define PS_CHECK_NULL_IMAGE(NAME, RETURN)                                                                     \
     92if (NAME == NULL || NAME->data.v == NULL) {                                                               \
     93    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                                \
     94    return RETURN;                                                                                        \
     95}
     96
     97/** Preprocessor macro to create image based on another */
     98#define PS_CHECK_ALLOC_IMAGE(NAME, NCOLS, NROWS, PS_TYPE)                                                    \
     99if(NAME == NULL) {                                                                                       \
     100    NAME = psImageAlloc(NCOLS, NROWS, PS_TYPE);                                                          \
     101}
     102
     103/** Preprocessor macro to generate error for zero length rows or columns */
     104#define PS_CHECK_SIZE_IMAGE(NAME, RETURN)                                                                     \
     105if (NAME->numCols < 1 || NAME->numRows < 1) {                                                             \
     106    psError(__func__,"Invalid operation: %s has zero rows or columns (%dx%d).", #NAME,                    \
     107            NAME->numCols, NAME->numRows);                                                                        \
     108    return RETURN;                                                                                        \
     109}
     110
     111/** Preprocessor macro to generate error for image dimensionality not set to PS_DIMEN_IMAGE */
     112#define PS_CHECK_DIMEN_AND_TYPE(NAME, PS_DIMEN, RETURN)                                                       \
     113if (NAME->type.dimen != PS_DIMEN) {                                                                       \
     114    psError(__func__,"Invalid operation: %s incorrect dimensionality %d.", #NAME, PS_DIMEN);              \
     115    return RETURN;                                                                                        \
     116} else if(NAME->type.type != PS_TYPE_F64) {                                                               \
     117    psError(__func__, "Invalid operation: %s not PS_TYPE_F64.", #NAME);                                   \
     118    return RETURN;                                                                                        \
     119}
     120
     121/** Preprocessor macro to check that input is not equal to output */
     122#define PS_CHECK_POINTERS(NAME1, NAME2, RETURN)                                                               \
     123if (NAME1 == NAME2) {                                                                                     \
     124    psError(__func__,"Invalid operation: Pointer to %s is same as %s.", #NAME1, #NAME2);                  \
     125    return RETURN;                                                                                        \
     126}
     127
     128/** Preprocessor macro to check that an image is square */
     129#define PS_CHECK_SQUARE(NAME, RETURN)                                                                         \
     130if (NAME->numCols != NAME->numRows) {                                                                     \
     131    psError(__func__,"Invalid operation: %s not square array.", #NAME);                                   \
     132    return RETURN;                                                                                        \
     133}
     134
     135/** Preprocessor macro to initalize a GSL matrix. */
     136#define PS_GSL_MATRIX_INITIALIZE(LHS_NAME, RHS_NAME)                                                          \
     137LHS_NAME.size1 = numRows;                                                                                 \
     138LHS_NAME.size2 = numCols;                                                                                 \
     139LHS_NAME.tda   = numCols;                                                                                 \
     140LHS_NAME.data  = RHS_NAME;
    71141
    72142/*****************************************************************************/
    73143/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    74144/*****************************************************************************/
    75 
    76 #define PS_MATRIX_TRANSPOSE(PS_TYPE)                                                                          \
    77 {                                                                                                             \
    78     int arraySize = 0;                                                                                        \
    79     int numRows = 0;                                                                                          \
    80     int numCols = 0;                                                                                          \
    81     gsl_matrix##PS_TYPE trans;                                                                                \
    82     \
    83     /* Initialize data */                                                                                     \
    84     numRows = inImage->numRows;                                                                               \
    85     numCols = inImage->numCols;                                                                               \
    86     arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols;                                        \
    87     \
    88     /* Copy psImage input data into psImage output data to keep input data pristine */                        \
    89     memcpy(outImage->data.v[0], inImage->data.v[0], arraySize);                                               \
    90     \
    91     /* Manually fill inverted output matrix so it will be aligned with output image */                        \
    92     trans.size1 = numRows;                                                                                    \
    93     trans.size2 = numCols;                                                                                    \
    94     trans.tda = numCols;                                                                                      \
    95     trans.data = outImage->data.v[0];                                                                         \
    96     \
    97     /* Transpose data */                                                                                      \
    98     gsl_matrix##PS_TYPE##_transpose(&trans);                                                                  \
    99 }
    100 
    101 psImage* psMatrixTranspose(psImage *outImage, const psImage *inImage)
    102 {
    103     psElemType elemType = 0;
    104 
    105     elemType = inImage->type.type;
    106     switch(elemType) {
    107     case PS_TYPE_FLOAT:
    108         PS_MATRIX_TRANSPOSE(_float);
    109         break;
    110     case PS_TYPE_DOUBLE:
    111         PS_MATRIX_TRANSPOSE();
    112         break;
    113     default:
    114         psError(__func__, " : Line %d - Invalid psElemType:  %d\n", __LINE__, elemType);
    115     }
    116 
    117     return outImage;
    118 }
    119 
    120 psImage* psMatrixOp(psImage *outImage, psImage *inImage1, const char op, psImage *inImage2)
    121 {
    122     int arraySize = 0;
    123     int numRows = 0;
    124     int numCols = 0;
    125     gsl_matrix m2;
    126     gsl_matrix m3;
    127 
    128     // Initialize data
    129     numRows = inImage1->numRows;
    130     numCols = inImage1->numCols;
    131     arraySize = PSELEMTYPE_SIZEOF(outImage->type.type)*numRows*numCols;
    132 
    133     m2.size1 = numRows;
    134     m2.size2 = numCols;
    135     m2.tda   = numCols;
    136     m2.data  = inImage2->data.v[0];
    137 
    138     m3.size1 = numRows;
    139     m3.size2 = numCols;
    140     m3.tda   = numCols;
    141     m3.data  = outImage->data.v[0];
    142 
    143     // Copy psImage input data into GSL matrix data to keep input data pristine
    144     memcpy(m3.data, inImage1->data.v[0], arraySize);
    145 
    146 
    147     switch(op) {
    148     case '+':
    149         gsl_matrix_add(&m3, &m2);
    150         break;
    151     case '-':
    152         gsl_matrix_sub(&m3, &m2);
    153         break;
    154     case '*':
    155         gsl_linalg_matmult(&m3, &m2, &m3);
    156         break;
    157     default:
    158         psError(__func__, " : Line %d - Invalid psMatrixOp operation: %c\n", __LINE__, op);
    159         \
    160     }
    161 
    162     return NULL;
    163 }
    164145
    165146psImage *psMatrixLUD(psImage *outImage, psVector *outPerm, psImage *inImage)
     
    171152    gsl_matrix lu;
    172153    gsl_permutation perm;
     154
     155    // Error checks
     156    PS_CHECK_POINTERS(inImage, outImage, outImage);
     157    PS_CHECK_NULL_IMAGE(inImage, outImage);
     158    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     159    PS_CHECK_SIZE_IMAGE(inImage, outImage);
     160    PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     161    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     162    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     163    PS_CHECK_ALLOC_VECTOR(outPerm, inImage->numRows, inImage->type.type);
     164    PS_CHECK_NULL_VECTOR(outPerm, outImage);
     165    PS_CHECK_DIMEN_AND_TYPE(outPerm, PS_DIMEN_VECTOR, outImage);
    173166
    174167    // Initialize data
     
    177170    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols;
    178171
    179     // Manually fill GSL structs so they will be aligned with output data
     172    // Initialize GSL data
    180173    perm.size = numCols;
    181174    outPerm->n = numCols;
    182175    perm.data = outPerm->vec.v;
    183     lu.size1 = numRows;
    184     lu.size2 = numCols;
    185     lu.tda = numCols;
    186     lu.data = outImage->data.v[0];
     176    PS_GSL_MATRIX_INITIALIZE(lu, outImage->data.v[0]);
     177
     178    // Non-square matrices not allowed
     179    PS_CHECK_SQUARE(inImage, outImage);
     180    PS_CHECK_SQUARE(outImage, outImage);
    187181
    188182    // Copy psImage input data into GSL matrix data to keep input data pristine
     
    195189}
    196190
    197 psVector *psMatrixLUSolve(psVector *outVector, const psImage *luImage, const psVector *inVector, const psVector *inPerm)
     191psVector *psMatrixLUSolve(psVector *outVector, const psImage *inImage, const psVector *inVector, const psVector *inPerm)
    198192{
    199193    int arraySize = 0;
     
    205199    gsl_vector x;
    206200
    207     // Initialize data
    208     numRows = luImage->numRows;
    209     numCols = luImage->numCols;
    210     arraySize = PSELEMTYPE_SIZEOF(luImage->type.type)*numRows*numCols;
    211 
    212     // Manually fill GSL structs so they will be aligned with output data
    213     lu.size1 = luImage->numRows;
    214     lu.size2 = luImage->numCols;
    215     lu.tda = luImage->numCols;
    216     lu.data = luImage->data.v[0];
     201    // Error checks
     202    PS_CHECK_POINTERS(outVector, inVector, outVector);
     203    PS_CHECK_POINTERS(inVector, inPerm, outVector);
     204    PS_CHECK_POINTERS(outVector, inPerm, outVector);
     205    PS_CHECK_NULL_IMAGE(inImage, outVector);
     206    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outVector);
     207    PS_CHECK_SIZE_IMAGE(inImage, outVector);
     208    PS_CHECK_NULL_VECTOR(outVector, outVector);
     209    PS_CHECK_DIMEN_AND_TYPE(outVector, PS_DIMEN_VECTOR, outVector);
     210    PS_CHECK_NULL_VECTOR(inVector, outVector);
     211    PS_CHECK_DIMEN_AND_TYPE(inVector, PS_DIMEN_VECTOR, outVector);
     212    PS_CHECK_NULL_VECTOR(inPerm, outVector);
     213    PS_CHECK_DIMEN_AND_TYPE(inPerm, PS_DIMEN_VECTOR, outVector);
     214
     215    // Initialize data
     216    numRows = inImage->numRows;
     217    numCols = inImage->numCols;
     218    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols;
     219
     220    // Initialize GSL data
     221    PS_GSL_MATRIX_INITIALIZE(lu, inImage->data.v[0]);
    217222
    218223    outVector->n = numCols;
     
    245250    gsl_permutation *perm = NULL;
    246251
     252    // Error checks
     253    PS_CHECK_POINTERS(inImage, outImage, outImage);
     254    PS_CHECK_NULL_IMAGE(inImage, outImage);
     255    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     256    PS_CHECK_SIZE_IMAGE(inImage, outImage);
     257    PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     258    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     259    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     260
    247261    // Initialize data
    248262    numRows = inImage->numRows;
     
    254268    lu = gsl_matrix_alloc(numRows, numCols);
    255269
     270    // Initialize GSL data
     271    PS_GSL_MATRIX_INITIALIZE(inv, outImage->data.v[0]);
     272
     273    // Non-square matrices not allowed
     274    PS_CHECK_SQUARE(inImage, outImage);
     275    PS_CHECK_SQUARE(outImage, outImage);
     276
    256277    // Copy psImage input data into GSL matrix data to keep input data pristine
    257278    memcpy(lu->data, inImage->data.v[0], arraySize);
    258 
    259     // Manually fill inverted output matrix so it will be aligned with output image
    260     inv.size1 = outImage->numRows;
    261     inv.size2 = outImage->numCols;
    262     inv.tda = outImage->numCols;
    263     inv.data = outImage->data.v[0];
    264279
    265280    // Invert data and calculate determinant
     
    285300    gsl_permutation *perm = NULL;
    286301
     302    // Error checks
     303    PS_CHECK_NULL_IMAGE(inImage, 0);
     304    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, 0);
     305    PS_CHECK_SIZE_IMAGE(inImage, 0);
     306
    287307    // Initialize data
    288308    numRows = inImage->numRows;
     
    294314    lu = gsl_matrix_alloc(numRows, numCols);
    295315
     316    // Non-square matrices not allowed
     317    PS_CHECK_SQUARE(inImage, 0);
     318
    296319    // Copy psImage input data into GSL matrix data to keep input data pristine
    297320    memcpy(lu->data, inImage->data.v[0], arraySize);
     
    306329
    307330    return det;
     331}
     332
     333psImage* psMatrixMultiply(psImage *outImage, psImage *inImage1, psImage *inImage2)
     334{
     335    int arraySize = 0;
     336    int numRows = 0;
     337    int numCols = 0;
     338    gsl_matrix m1;
     339    gsl_matrix m2;
     340    gsl_matrix m3;
     341
     342    // Error checks
     343    PS_CHECK_POINTERS(inImage1, outImage, outImage);
     344    PS_CHECK_POINTERS(inImage1, inImage2, outImage);
     345    PS_CHECK_NULL_IMAGE(inImage1, outImage);
     346    PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, outImage);
     347    PS_CHECK_SIZE_IMAGE(inImage1, outImage);
     348    PS_CHECK_NULL_IMAGE(inImage2, outImage);
     349    PS_CHECK_DIMEN_AND_TYPE(inImage2, PS_DIMEN_IMAGE, outImage);
     350    PS_CHECK_SIZE_IMAGE(inImage2, outImage);
     351    PS_CHECK_ALLOC_IMAGE(outImage, inImage2->numCols, inImage2->numRows, inImage2->type.type);
     352    PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, outImage);
     353    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     354
     355    // Initialize data
     356    numRows = inImage1->numRows;
     357    numCols = inImage1->numCols;
     358    arraySize = PSELEMTYPE_SIZEOF(outImage->type.type)*numRows*numCols;
     359
     360    // Initialize GSL data
     361    PS_GSL_MATRIX_INITIALIZE(m1, inImage1->data.v[0]);
     362    PS_GSL_MATRIX_INITIALIZE(m2, inImage2->data.v[0]);
     363    PS_GSL_MATRIX_INITIALIZE(m3, outImage->data.v[0]);
     364
     365    // Non-square matrices not allowed
     366    PS_CHECK_SQUARE(inImage1, outImage);
     367    PS_CHECK_SQUARE(inImage2, outImage);
     368    PS_CHECK_SQUARE(outImage, outImage);
     369
     370    gsl_linalg_matmult(&m1, &m2, &m3);
     371
     372    return outImage;
     373}
     374
     375psImage* psMatrixTranspose(psImage *outImage, const psImage *inImage)
     376{
     377    int arraySize = 0;
     378    int numRows = 0;
     379    int numCols = 0;
     380    gsl_matrix trans;
     381
     382    // Error checks
     383    PS_CHECK_POINTERS(inImage, outImage, outImage);
     384    PS_CHECK_NULL_IMAGE(inImage, outImage);
     385    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     386    PS_CHECK_SIZE_IMAGE(inImage, outImage);
     387    PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     388    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     389    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     390
     391    // Initialize data
     392    numRows = inImage->numRows;
     393    numCols = inImage->numCols;
     394    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols;
     395
     396    // Initialize GSL data
     397    PS_GSL_MATRIX_INITIALIZE(trans, outImage->data.v[0]);
     398
     399    // Non-square matrices not allowed
     400    PS_CHECK_SQUARE(inImage, outImage);
     401    PS_CHECK_SQUARE(outImage, outImage);
     402
     403    // Copy psImage input data into psImage output data to keep input data pristine
     404    memcpy(outImage->data.v[0], inImage->data.v[0], arraySize);
     405
     406    // Transpose data
     407    gsl_matrix_transpose(&trans);
     408
     409    return outImage;
    308410}
    309411
     
    317419    gsl_matrix in;
    318420
     421    // Error checks
     422    PS_CHECK_POINTERS(inImage, outImage, outImage);
     423    PS_CHECK_NULL_IMAGE(inImage, outImage);
     424    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     425    PS_CHECK_SIZE_IMAGE(inImage, outImage);
     426    PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     427    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     428    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     429
    319430    // Initialize data
    320431    numRows = inImage->numRows;
    321432    numCols = inImage->numCols;
    322     out.data = outImage->data.v[0];
    323 
    324     // Manually fill GSL structs so they will be aligned with output data
    325     in.size1 = numRows;
    326     in.size2 = numCols;
    327     in.tda = numCols;
    328     in.data = inImage->data.v[0];
    329 
    330     out.size1 = numRows;
    331     out.size2 = numCols;
    332     out.tda = numCols;
    333     out.data = outImage->data.v[0];
     433
     434    // Initialize GSL data
     435    PS_GSL_MATRIX_INITIALIZE(in, inImage->data.v[0]);
     436    PS_GSL_MATRIX_INITIALIZE(out, outImage->data.v[0]);
    334437
    335438    // Allocate GSL structs
     
    337440    w = gsl_eigen_symmv_alloc(numRows);
    338441
    339     // Calculate Eigenvalues and Eigenvectors. Eigenvalues not currently used
     442    // Non-square matrices not allowed
     443    PS_CHECK_SQUARE(inImage, outImage);
     444    PS_CHECK_SQUARE(outImage, outImage);
     445
     446    // Calculate Eigenvalues and Eigenvectors...Eigenvalues not currently used
    340447    gsl_eigen_symmv(&in, eVals, &out, w);
    341448
     
    349456psVector *psMatrixToVector(psVector *outVector, psImage *inImage)
    350457{
    351     // if inimage isn't 1d
    352     // if inimage nrows != vector len
    353 
    354458    int colSize = 0;
     459
     460    // Error checks
     461    PS_CHECK_NULL_IMAGE(inImage, outVector);
     462    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outVector);
     463    PS_CHECK_SIZE_IMAGE(inImage, outVector);
     464    PS_CHECK_ALLOC_VECTOR(outVector, inImage->numRows, inImage->type.type);
     465    PS_CHECK_NULL_VECTOR(outVector, outVector);
     466    PS_CHECK_DIMEN_AND_TYPE(outVector, PS_DIMEN_VECTOR, outVector);
     467
     468    // Set n if allocated, but empty
     469    if(outVector->n == 0) {
     470        outVector->n = inImage->numRows;
     471    }
     472
     473    // More checks
     474    if(inImage->numCols > 1) {
     475        psError(__func__, "Image has more than 1 column: numCols = %d.", inImage->numCols);
     476        return outVector;
     477    } else if(outVector->n != inImage->numRows) {
     478        psError(__func__, "Image and vector sizes differ: (%d vs %d).", inImage->numRows, outVector->n);
     479        return outVector;
     480    }
    355481
    356482    colSize = PSELEMTYPE_SIZEOF(inImage->type.type)*inImage->numRows;
     
    362488psImage *psVectorToMatrix(psImage *outImage, psVector *inVector)
    363489{
    364     // if inimage isn't 1d
    365     // if inimage nrows != vector len
    366 
    367490    int colSize = 0;
     491
     492    // Error checks
     493    PS_CHECK_NULL_VECTOR(inVector, outImage);
     494    PS_CHECK_DIMEN_AND_TYPE(inVector, PS_DIMEN_VECTOR, outImage);
     495    PS_CHECK_SIZE_VECTOR(inVector, outImage);
     496    PS_CHECK_ALLOC_IMAGE(outImage, 1, inVector->n, PS_TYPE_F32)
     497    PS_CHECK_NULL_IMAGE(outImage, outImage);
     498    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     499
     500    // More checks
     501    if(outImage->numCols > 1) {
     502        psError(__func__, "Image has more than 1 column: numCols = %d.", outImage->numCols);
     503        return outImage;
     504    } else if(outImage->numRows != inVector->n) {
     505        psError(__func__, "Image and vector sizes differ: (%d vs %d).", outImage->numRows, inVector->n);
     506        return outImage;
     507    }
    368508
    369509    colSize = PSELEMTYPE_SIZEOF(outImage->type.type)*outImage->numRows;
  • trunk/psLib/src/dataManip/psMatrix.h

    r760 r799  
    2121 *  @author Ross Harman, MHPCC
    2222 *   
    23  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    24  *  @date $Date: 2004-05-24 21:10:03 $
     23 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     24 *  @date $Date: 2004-05-28 02:53:27 $
    2525 *
    2626 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    4949 *
    5050 *  Performs a LU decomposition on a psImage matrix and returns the LU matrix. If the user specifies NULL as the
    51  *  outImage argument, then the decomposition is done in-place and inImage will hold the inverted image. The input
    52  *  image must be square. This function operates only with psF32 and psF64 data types.
     51 *  outImage argument, then the output argument will be created based on input values. The input image must
     52 *  be square. This function operates only with the psF32 data type.
    5353 *
    5454 *  @return  psImage*: Pointer to LU decomposed psImage.
     
    102102/** Performs basic psImage matrix operations.
    103103 *
    104  *  Performs psImage matrix operations for: addition, subtraction, multiplication. For addition, each element
    105  *  of the input image is added to the corresponding element of the output matrix. Subtraction works in a
    106  *  similar manner. For multiplication, the function performs a classical matrix multiplication involving row
    107  *  and column operations. For matrix multiplication, the number of columns must match the number of rows for
    108  *  inImage1 and inImage2, respectively. Matrix division is not defined for this function. If the user
    109  *  specifies NULL as the outImage argument, then a new psImage will be created and returned. This function
    110  *  operates only with psF32 and psF64 data types. 
     104 *  Performs a classical matrix multiplication involving row and column operations. This function assumes both
     105 *  matrices are the same size. If the user specifies NULL as the outImage argument, then a new psImage will be
     106 *  created and returned. This function operates only with the psF64 data type. GSL indexes the top row as the
     107 *  zero row, not the bottom. 
    111108 *
    112109 *  @return  psImage*: Pointer to resulting psImage.
    113110 */
    114 psImage *psMatrixOp(
     111psImage *psMatrixMultiply(
    115112    psImage *outImage,  ///< Matrix to return, or NULL.
    116113    psImage *inImage1,  ///< First input image.
    117     const char op,      ///< Operation to perform: "+", "-", "*".
    118114    psImage *inImage2   ///< Second input image.
    119115);
     
    134130/** Calculate matrix eigenvectors.
    135131 *
    136  *  Calculates the eigenvectors for a matrix. The input image must be square. If the user specifies NULL as
     132 *  Calculates the eigenvectors for a matrix. The input image must be symmetric. If the user specifies NULL as
    137133 *  the outImage argument, then a new psImage will be created and returned. This function operates only with
    138  *  psF32 and psF64 data types.
     134 *  the psF64 data type.
    139135 *
    140136 *  @return  psImage*: Pointer to matrix of Eigenvectors.
    141137 */
    142138psImage *psMatrixEigenvectors(
    143     psImage *outImage, ///< Eigenvectors to return, or NULL.
     139    psImage *outImage,  ///< Eigenvectors to return, or NULL.
    144140    psImage *inImage    ///< Input image.
    145141);
     
    147143/** Convert matrix to vector.
    148144 *
    149  *  Converts a 1-d matrix into a vector. If the user specifies NULL as the outVector argument, then a new psImage will be created
    150  *  and returned. This function operates only with psF32 and psF64 data types. 
     145 *  Converts a 1-d psImage matrix into a vector. If the user specifies NULL as the outVector argument, then a
     146 *  new psVector will be created. The input matrix must be a 1-d column matrix. This function operates only with
     147 *  the psF64 data type
    151148 *
    152149 *  @return  psVector*: Pointer to psVector.
     
    159156/** Convert vector to matrix.
    160157 *
    161  *  Converts a vector to a 1-d psImage matrix into a vector. If the user specifies NULL as the outImage
    162  *  argument, then a new psImage will be created and returned. This function operates only with psF32 and
    163  *  psF64 data types. 
     158 *  Converts a vector into a 1-d column psImage matrix. If the user specifies NULL as the outImage argument,
     159 *  then a new psImage will be created. This function operates only with the psF64 data type. 
    164160 *
    165161 *  @return  psVector*: Pointer to psIamge.
  • trunk/psLib/src/math/psMatrix.c

    r770 r799  
    2020 *  @author Ross Harman, MHPCC
    2121 *   
    22  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    23  *  @date $Date: 2004-05-24 23:30:52 $
     22 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     23 *  @date $Date: 2004-05-28 02:53:18 $
    2424 *
    2525 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    6868/*****************************************************************************/
    6969
    70 // None
     70/** Preprocessor macro to generate error a NULL image */
     71#define PS_CHECK_NULL_VECTOR(NAME, RETURN)                                                                    \
     72if (NAME == NULL || NAME->vec.v == NULL) {                                                                \
     73    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                                \
     74    return RETURN;                                                                                        \
     75}
     76
     77/** Preprocessor macro to create vector based on another */
     78#define PS_CHECK_ALLOC_VECTOR(NAME, SIZE, PS_TYPE)                                                            \
     79if(NAME == NULL) {                                                                                        \
     80    NAME = psVectorAlloc(PS_TYPE, SIZE);                                                                  \
     81}
     82
     83/** Preprocessor macro to generate error for zero length vector */
     84#define PS_CHECK_SIZE_VECTOR(NAME, RETURN)                                                                    \
     85if (NAME->n < 1) {                                                                                        \
     86    psError(__func__,"Invalid operation: %s has zero n value.", #NAME);                                   \
     87    return RETURN;                                                                                        \
     88}
     89
     90/** Preprocessor macro to generate error a NULL image */
     91#define PS_CHECK_NULL_IMAGE(NAME, RETURN)                                                                     \
     92if (NAME == NULL || NAME->data.v == NULL) {                                                               \
     93    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                                \
     94    return RETURN;                                                                                        \
     95}
     96
     97/** Preprocessor macro to create image based on another */
     98#define PS_CHECK_ALLOC_IMAGE(NAME, NCOLS, NROWS, PS_TYPE)                                                    \
     99if(NAME == NULL) {                                                                                       \
     100    NAME = psImageAlloc(NCOLS, NROWS, PS_TYPE);                                                          \
     101}
     102
     103/** Preprocessor macro to generate error for zero length rows or columns */
     104#define PS_CHECK_SIZE_IMAGE(NAME, RETURN)                                                                     \
     105if (NAME->numCols < 1 || NAME->numRows < 1) {                                                             \
     106    psError(__func__,"Invalid operation: %s has zero rows or columns (%dx%d).", #NAME,                    \
     107            NAME->numCols, NAME->numRows);                                                                        \
     108    return RETURN;                                                                                        \
     109}
     110
     111/** Preprocessor macro to generate error for image dimensionality not set to PS_DIMEN_IMAGE */
     112#define PS_CHECK_DIMEN_AND_TYPE(NAME, PS_DIMEN, RETURN)                                                       \
     113if (NAME->type.dimen != PS_DIMEN) {                                                                       \
     114    psError(__func__,"Invalid operation: %s incorrect dimensionality %d.", #NAME, PS_DIMEN);              \
     115    return RETURN;                                                                                        \
     116} else if(NAME->type.type != PS_TYPE_F64) {                                                               \
     117    psError(__func__, "Invalid operation: %s not PS_TYPE_F64.", #NAME);                                   \
     118    return RETURN;                                                                                        \
     119}
     120
     121/** Preprocessor macro to check that input is not equal to output */
     122#define PS_CHECK_POINTERS(NAME1, NAME2, RETURN)                                                               \
     123if (NAME1 == NAME2) {                                                                                     \
     124    psError(__func__,"Invalid operation: Pointer to %s is same as %s.", #NAME1, #NAME2);                  \
     125    return RETURN;                                                                                        \
     126}
     127
     128/** Preprocessor macro to check that an image is square */
     129#define PS_CHECK_SQUARE(NAME, RETURN)                                                                         \
     130if (NAME->numCols != NAME->numRows) {                                                                     \
     131    psError(__func__,"Invalid operation: %s not square array.", #NAME);                                   \
     132    return RETURN;                                                                                        \
     133}
     134
     135/** Preprocessor macro to initalize a GSL matrix. */
     136#define PS_GSL_MATRIX_INITIALIZE(LHS_NAME, RHS_NAME)                                                          \
     137LHS_NAME.size1 = numRows;                                                                                 \
     138LHS_NAME.size2 = numCols;                                                                                 \
     139LHS_NAME.tda   = numCols;                                                                                 \
     140LHS_NAME.data  = RHS_NAME;
    71141
    72142/*****************************************************************************/
    73143/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    74144/*****************************************************************************/
    75 
    76 #define PS_MATRIX_TRANSPOSE(PS_TYPE)                                                                          \
    77 {                                                                                                             \
    78     int arraySize = 0;                                                                                        \
    79     int numRows = 0;                                                                                          \
    80     int numCols = 0;                                                                                          \
    81     gsl_matrix##PS_TYPE trans;                                                                                \
    82     \
    83     /* Initialize data */                                                                                     \
    84     numRows = inImage->numRows;                                                                               \
    85     numCols = inImage->numCols;                                                                               \
    86     arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols;                                        \
    87     \
    88     /* Copy psImage input data into psImage output data to keep input data pristine */                        \
    89     memcpy(outImage->data.v[0], inImage->data.v[0], arraySize);                                               \
    90     \
    91     /* Manually fill inverted output matrix so it will be aligned with output image */                        \
    92     trans.size1 = numRows;                                                                                    \
    93     trans.size2 = numCols;                                                                                    \
    94     trans.tda = numCols;                                                                                      \
    95     trans.data = outImage->data.v[0];                                                                         \
    96     \
    97     /* Transpose data */                                                                                      \
    98     gsl_matrix##PS_TYPE##_transpose(&trans);                                                                  \
    99 }
    100 
    101 psImage* psMatrixTranspose(psImage *outImage, const psImage *inImage)
    102 {
    103     psElemType elemType = 0;
    104 
    105     elemType = inImage->type.type;
    106     switch(elemType) {
    107     case PS_TYPE_FLOAT:
    108         PS_MATRIX_TRANSPOSE(_float);
    109         break;
    110     case PS_TYPE_DOUBLE:
    111         PS_MATRIX_TRANSPOSE();
    112         break;
    113     default:
    114         psError(__func__, " : Line %d - Invalid psElemType:  %d\n", __LINE__, elemType);
    115     }
    116 
    117     return outImage;
    118 }
    119 
    120 psImage* psMatrixOp(psImage *outImage, psImage *inImage1, const char op, psImage *inImage2)
    121 {
    122     int arraySize = 0;
    123     int numRows = 0;
    124     int numCols = 0;
    125     gsl_matrix m2;
    126     gsl_matrix m3;
    127 
    128     // Initialize data
    129     numRows = inImage1->numRows;
    130     numCols = inImage1->numCols;
    131     arraySize = PSELEMTYPE_SIZEOF(outImage->type.type)*numRows*numCols;
    132 
    133     m2.size1 = numRows;
    134     m2.size2 = numCols;
    135     m2.tda   = numCols;
    136     m2.data  = inImage2->data.v[0];
    137 
    138     m3.size1 = numRows;
    139     m3.size2 = numCols;
    140     m3.tda   = numCols;
    141     m3.data  = outImage->data.v[0];
    142 
    143     // Copy psImage input data into GSL matrix data to keep input data pristine
    144     memcpy(m3.data, inImage1->data.v[0], arraySize);
    145 
    146 
    147     switch(op) {
    148     case '+':
    149         gsl_matrix_add(&m3, &m2);
    150         break;
    151     case '-':
    152         gsl_matrix_sub(&m3, &m2);
    153         break;
    154     case '*':
    155         gsl_linalg_matmult(&m3, &m2, &m3);
    156         break;
    157     default:
    158         psError(__func__, " : Line %d - Invalid psMatrixOp operation: %c\n", __LINE__, op);
    159         \
    160     }
    161 
    162     return NULL;
    163 }
    164145
    165146psImage *psMatrixLUD(psImage *outImage, psVector *outPerm, psImage *inImage)
     
    171152    gsl_matrix lu;
    172153    gsl_permutation perm;
     154
     155    // Error checks
     156    PS_CHECK_POINTERS(inImage, outImage, outImage);
     157    PS_CHECK_NULL_IMAGE(inImage, outImage);
     158    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     159    PS_CHECK_SIZE_IMAGE(inImage, outImage);
     160    PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     161    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     162    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     163    PS_CHECK_ALLOC_VECTOR(outPerm, inImage->numRows, inImage->type.type);
     164    PS_CHECK_NULL_VECTOR(outPerm, outImage);
     165    PS_CHECK_DIMEN_AND_TYPE(outPerm, PS_DIMEN_VECTOR, outImage);
    173166
    174167    // Initialize data
     
    177170    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols;
    178171
    179     // Manually fill GSL structs so they will be aligned with output data
     172    // Initialize GSL data
    180173    perm.size = numCols;
    181174    outPerm->n = numCols;
    182175    perm.data = outPerm->vec.v;
    183     lu.size1 = numRows;
    184     lu.size2 = numCols;
    185     lu.tda = numCols;
    186     lu.data = outImage->data.v[0];
     176    PS_GSL_MATRIX_INITIALIZE(lu, outImage->data.v[0]);
     177
     178    // Non-square matrices not allowed
     179    PS_CHECK_SQUARE(inImage, outImage);
     180    PS_CHECK_SQUARE(outImage, outImage);
    187181
    188182    // Copy psImage input data into GSL matrix data to keep input data pristine
     
    195189}
    196190
    197 psVector *psMatrixLUSolve(psVector *outVector, const psImage *luImage, const psVector *inVector, const psVector *inPerm)
     191psVector *psMatrixLUSolve(psVector *outVector, const psImage *inImage, const psVector *inVector, const psVector *inPerm)
    198192{
    199193    int arraySize = 0;
     
    205199    gsl_vector x;
    206200
    207     // Initialize data
    208     numRows = luImage->numRows;
    209     numCols = luImage->numCols;
    210     arraySize = PSELEMTYPE_SIZEOF(luImage->type.type)*numRows*numCols;
    211 
    212     // Manually fill GSL structs so they will be aligned with output data
    213     lu.size1 = luImage->numRows;
    214     lu.size2 = luImage->numCols;
    215     lu.tda = luImage->numCols;
    216     lu.data = luImage->data.v[0];
     201    // Error checks
     202    PS_CHECK_POINTERS(outVector, inVector, outVector);
     203    PS_CHECK_POINTERS(inVector, inPerm, outVector);
     204    PS_CHECK_POINTERS(outVector, inPerm, outVector);
     205    PS_CHECK_NULL_IMAGE(inImage, outVector);
     206    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outVector);
     207    PS_CHECK_SIZE_IMAGE(inImage, outVector);
     208    PS_CHECK_NULL_VECTOR(outVector, outVector);
     209    PS_CHECK_DIMEN_AND_TYPE(outVector, PS_DIMEN_VECTOR, outVector);
     210    PS_CHECK_NULL_VECTOR(inVector, outVector);
     211    PS_CHECK_DIMEN_AND_TYPE(inVector, PS_DIMEN_VECTOR, outVector);
     212    PS_CHECK_NULL_VECTOR(inPerm, outVector);
     213    PS_CHECK_DIMEN_AND_TYPE(inPerm, PS_DIMEN_VECTOR, outVector);
     214
     215    // Initialize data
     216    numRows = inImage->numRows;
     217    numCols = inImage->numCols;
     218    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols;
     219
     220    // Initialize GSL data
     221    PS_GSL_MATRIX_INITIALIZE(lu, inImage->data.v[0]);
    217222
    218223    outVector->n = numCols;
     
    245250    gsl_permutation *perm = NULL;
    246251
     252    // Error checks
     253    PS_CHECK_POINTERS(inImage, outImage, outImage);
     254    PS_CHECK_NULL_IMAGE(inImage, outImage);
     255    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     256    PS_CHECK_SIZE_IMAGE(inImage, outImage);
     257    PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     258    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     259    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     260
    247261    // Initialize data
    248262    numRows = inImage->numRows;
     
    254268    lu = gsl_matrix_alloc(numRows, numCols);
    255269
     270    // Initialize GSL data
     271    PS_GSL_MATRIX_INITIALIZE(inv, outImage->data.v[0]);
     272
     273    // Non-square matrices not allowed
     274    PS_CHECK_SQUARE(inImage, outImage);
     275    PS_CHECK_SQUARE(outImage, outImage);
     276
    256277    // Copy psImage input data into GSL matrix data to keep input data pristine
    257278    memcpy(lu->data, inImage->data.v[0], arraySize);
    258 
    259     // Manually fill inverted output matrix so it will be aligned with output image
    260     inv.size1 = outImage->numRows;
    261     inv.size2 = outImage->numCols;
    262     inv.tda = outImage->numCols;
    263     inv.data = outImage->data.v[0];
    264279
    265280    // Invert data and calculate determinant
     
    285300    gsl_permutation *perm = NULL;
    286301
     302    // Error checks
     303    PS_CHECK_NULL_IMAGE(inImage, 0);
     304    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, 0);
     305    PS_CHECK_SIZE_IMAGE(inImage, 0);
     306
    287307    // Initialize data
    288308    numRows = inImage->numRows;
     
    294314    lu = gsl_matrix_alloc(numRows, numCols);
    295315
     316    // Non-square matrices not allowed
     317    PS_CHECK_SQUARE(inImage, 0);
     318
    296319    // Copy psImage input data into GSL matrix data to keep input data pristine
    297320    memcpy(lu->data, inImage->data.v[0], arraySize);
     
    306329
    307330    return det;
     331}
     332
     333psImage* psMatrixMultiply(psImage *outImage, psImage *inImage1, psImage *inImage2)
     334{
     335    int arraySize = 0;
     336    int numRows = 0;
     337    int numCols = 0;
     338    gsl_matrix m1;
     339    gsl_matrix m2;
     340    gsl_matrix m3;
     341
     342    // Error checks
     343    PS_CHECK_POINTERS(inImage1, outImage, outImage);
     344    PS_CHECK_POINTERS(inImage1, inImage2, outImage);
     345    PS_CHECK_NULL_IMAGE(inImage1, outImage);
     346    PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, outImage);
     347    PS_CHECK_SIZE_IMAGE(inImage1, outImage);
     348    PS_CHECK_NULL_IMAGE(inImage2, outImage);
     349    PS_CHECK_DIMEN_AND_TYPE(inImage2, PS_DIMEN_IMAGE, outImage);
     350    PS_CHECK_SIZE_IMAGE(inImage2, outImage);
     351    PS_CHECK_ALLOC_IMAGE(outImage, inImage2->numCols, inImage2->numRows, inImage2->type.type);
     352    PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, outImage);
     353    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     354
     355    // Initialize data
     356    numRows = inImage1->numRows;
     357    numCols = inImage1->numCols;
     358    arraySize = PSELEMTYPE_SIZEOF(outImage->type.type)*numRows*numCols;
     359
     360    // Initialize GSL data
     361    PS_GSL_MATRIX_INITIALIZE(m1, inImage1->data.v[0]);
     362    PS_GSL_MATRIX_INITIALIZE(m2, inImage2->data.v[0]);
     363    PS_GSL_MATRIX_INITIALIZE(m3, outImage->data.v[0]);
     364
     365    // Non-square matrices not allowed
     366    PS_CHECK_SQUARE(inImage1, outImage);
     367    PS_CHECK_SQUARE(inImage2, outImage);
     368    PS_CHECK_SQUARE(outImage, outImage);
     369
     370    gsl_linalg_matmult(&m1, &m2, &m3);
     371
     372    return outImage;
     373}
     374
     375psImage* psMatrixTranspose(psImage *outImage, const psImage *inImage)
     376{
     377    int arraySize = 0;
     378    int numRows = 0;
     379    int numCols = 0;
     380    gsl_matrix trans;
     381
     382    // Error checks
     383    PS_CHECK_POINTERS(inImage, outImage, outImage);
     384    PS_CHECK_NULL_IMAGE(inImage, outImage);
     385    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     386    PS_CHECK_SIZE_IMAGE(inImage, outImage);
     387    PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     388    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     389    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     390
     391    // Initialize data
     392    numRows = inImage->numRows;
     393    numCols = inImage->numCols;
     394    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols;
     395
     396    // Initialize GSL data
     397    PS_GSL_MATRIX_INITIALIZE(trans, outImage->data.v[0]);
     398
     399    // Non-square matrices not allowed
     400    PS_CHECK_SQUARE(inImage, outImage);
     401    PS_CHECK_SQUARE(outImage, outImage);
     402
     403    // Copy psImage input data into psImage output data to keep input data pristine
     404    memcpy(outImage->data.v[0], inImage->data.v[0], arraySize);
     405
     406    // Transpose data
     407    gsl_matrix_transpose(&trans);
     408
     409    return outImage;
    308410}
    309411
     
    317419    gsl_matrix in;
    318420
     421    // Error checks
     422    PS_CHECK_POINTERS(inImage, outImage, outImage);
     423    PS_CHECK_NULL_IMAGE(inImage, outImage);
     424    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     425    PS_CHECK_SIZE_IMAGE(inImage, outImage);
     426    PS_CHECK_ALLOC_IMAGE(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     427    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     428    PS_CHECK_SIZE_IMAGE(outImage, outImage);
     429
    319430    // Initialize data
    320431    numRows = inImage->numRows;
    321432    numCols = inImage->numCols;
    322     out.data = outImage->data.v[0];
    323 
    324     // Manually fill GSL structs so they will be aligned with output data
    325     in.size1 = numRows;
    326     in.size2 = numCols;
    327     in.tda = numCols;
    328     in.data = inImage->data.v[0];
    329 
    330     out.size1 = numRows;
    331     out.size2 = numCols;
    332     out.tda = numCols;
    333     out.data = outImage->data.v[0];
     433
     434    // Initialize GSL data
     435    PS_GSL_MATRIX_INITIALIZE(in, inImage->data.v[0]);
     436    PS_GSL_MATRIX_INITIALIZE(out, outImage->data.v[0]);
    334437
    335438    // Allocate GSL structs
     
    337440    w = gsl_eigen_symmv_alloc(numRows);
    338441
    339     // Calculate Eigenvalues and Eigenvectors. Eigenvalues not currently used
     442    // Non-square matrices not allowed
     443    PS_CHECK_SQUARE(inImage, outImage);
     444    PS_CHECK_SQUARE(outImage, outImage);
     445
     446    // Calculate Eigenvalues and Eigenvectors...Eigenvalues not currently used
    340447    gsl_eigen_symmv(&in, eVals, &out, w);
    341448
     
    349456psVector *psMatrixToVector(psVector *outVector, psImage *inImage)
    350457{
    351     // if inimage isn't 1d
    352     // if inimage nrows != vector len
    353 
    354458    int colSize = 0;
     459
     460    // Error checks
     461    PS_CHECK_NULL_IMAGE(inImage, outVector);
     462    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outVector);
     463    PS_CHECK_SIZE_IMAGE(inImage, outVector);
     464    PS_CHECK_ALLOC_VECTOR(outVector, inImage->numRows, inImage->type.type);
     465    PS_CHECK_NULL_VECTOR(outVector, outVector);
     466    PS_CHECK_DIMEN_AND_TYPE(outVector, PS_DIMEN_VECTOR, outVector);
     467
     468    // Set n if allocated, but empty
     469    if(outVector->n == 0) {
     470        outVector->n = inImage->numRows;
     471    }
     472
     473    // More checks
     474    if(inImage->numCols > 1) {
     475        psError(__func__, "Image has more than 1 column: numCols = %d.", inImage->numCols);
     476        return outVector;
     477    } else if(outVector->n != inImage->numRows) {
     478        psError(__func__, "Image and vector sizes differ: (%d vs %d).", inImage->numRows, outVector->n);
     479        return outVector;
     480    }
    355481
    356482    colSize = PSELEMTYPE_SIZEOF(inImage->type.type)*inImage->numRows;
     
    362488psImage *psVectorToMatrix(psImage *outImage, psVector *inVector)
    363489{
    364     // if inimage isn't 1d
    365     // if inimage nrows != vector len
    366 
    367490    int colSize = 0;
     491
     492    // Error checks
     493    PS_CHECK_NULL_VECTOR(inVector, outImage);
     494    PS_CHECK_DIMEN_AND_TYPE(inVector, PS_DIMEN_VECTOR, outImage);
     495    PS_CHECK_SIZE_VECTOR(inVector, outImage);
     496    PS_CHECK_ALLOC_IMAGE(outImage, 1, inVector->n, PS_TYPE_F32)
     497    PS_CHECK_NULL_IMAGE(outImage, outImage);
     498    PS_CHECK_DIMEN_AND_TYPE(outImage, PS_DIMEN_IMAGE, outImage);
     499
     500    // More checks
     501    if(outImage->numCols > 1) {
     502        psError(__func__, "Image has more than 1 column: numCols = %d.", outImage->numCols);
     503        return outImage;
     504    } else if(outImage->numRows != inVector->n) {
     505        psError(__func__, "Image and vector sizes differ: (%d vs %d).", outImage->numRows, inVector->n);
     506        return outImage;
     507    }
    368508
    369509    colSize = PSELEMTYPE_SIZEOF(outImage->type.type)*outImage->numRows;
  • trunk/psLib/src/math/psMatrix.h

    r760 r799  
    2121 *  @author Ross Harman, MHPCC
    2222 *   
    23  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    24  *  @date $Date: 2004-05-24 21:10:03 $
     23 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     24 *  @date $Date: 2004-05-28 02:53:27 $
    2525 *
    2626 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    4949 *
    5050 *  Performs a LU decomposition on a psImage matrix and returns the LU matrix. If the user specifies NULL as the
    51  *  outImage argument, then the decomposition is done in-place and inImage will hold the inverted image. The input
    52  *  image must be square. This function operates only with psF32 and psF64 data types.
     51 *  outImage argument, then the output argument will be created based on input values. The input image must
     52 *  be square. This function operates only with the psF32 data type.
    5353 *
    5454 *  @return  psImage*: Pointer to LU decomposed psImage.
     
    102102/** Performs basic psImage matrix operations.
    103103 *
    104  *  Performs psImage matrix operations for: addition, subtraction, multiplication. For addition, each element
    105  *  of the input image is added to the corresponding element of the output matrix. Subtraction works in a
    106  *  similar manner. For multiplication, the function performs a classical matrix multiplication involving row
    107  *  and column operations. For matrix multiplication, the number of columns must match the number of rows for
    108  *  inImage1 and inImage2, respectively. Matrix division is not defined for this function. If the user
    109  *  specifies NULL as the outImage argument, then a new psImage will be created and returned. This function
    110  *  operates only with psF32 and psF64 data types. 
     104 *  Performs a classical matrix multiplication involving row and column operations. This function assumes both
     105 *  matrices are the same size. If the user specifies NULL as the outImage argument, then a new psImage will be
     106 *  created and returned. This function operates only with the psF64 data type. GSL indexes the top row as the
     107 *  zero row, not the bottom. 
    111108 *
    112109 *  @return  psImage*: Pointer to resulting psImage.
    113110 */
    114 psImage *psMatrixOp(
     111psImage *psMatrixMultiply(
    115112    psImage *outImage,  ///< Matrix to return, or NULL.
    116113    psImage *inImage1,  ///< First input image.
    117     const char op,      ///< Operation to perform: "+", "-", "*".
    118114    psImage *inImage2   ///< Second input image.
    119115);
     
    134130/** Calculate matrix eigenvectors.
    135131 *
    136  *  Calculates the eigenvectors for a matrix. The input image must be square. If the user specifies NULL as
     132 *  Calculates the eigenvectors for a matrix. The input image must be symmetric. If the user specifies NULL as
    137133 *  the outImage argument, then a new psImage will be created and returned. This function operates only with
    138  *  psF32 and psF64 data types.
     134 *  the psF64 data type.
    139135 *
    140136 *  @return  psImage*: Pointer to matrix of Eigenvectors.
    141137 */
    142138psImage *psMatrixEigenvectors(
    143     psImage *outImage, ///< Eigenvectors to return, or NULL.
     139    psImage *outImage,  ///< Eigenvectors to return, or NULL.
    144140    psImage *inImage    ///< Input image.
    145141);
     
    147143/** Convert matrix to vector.
    148144 *
    149  *  Converts a 1-d matrix into a vector. If the user specifies NULL as the outVector argument, then a new psImage will be created
    150  *  and returned. This function operates only with psF32 and psF64 data types. 
     145 *  Converts a 1-d psImage matrix into a vector. If the user specifies NULL as the outVector argument, then a
     146 *  new psVector will be created. The input matrix must be a 1-d column matrix. This function operates only with
     147 *  the psF64 data type
    151148 *
    152149 *  @return  psVector*: Pointer to psVector.
     
    159156/** Convert vector to matrix.
    160157 *
    161  *  Converts a vector to a 1-d psImage matrix into a vector. If the user specifies NULL as the outImage
    162  *  argument, then a new psImage will be created and returned. This function operates only with psF32 and
    163  *  psF64 data types. 
     158 *  Converts a vector into a 1-d column psImage matrix. If the user specifies NULL as the outImage argument,
     159 *  then a new psImage will be created. This function operates only with the psF64 data type. 
    164160 *
    165161 *  @return  psVector*: Pointer to psIamge.
Note: See TracChangeset for help on using the changeset viewer.