IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 9, 2004, 10:51:35 AM (22 years ago)
Author:
harman
Message:

Added support for psF32 types

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psMatrix.c

    r2273 r2671  
    11/** @file  psMatrix.c
    22 *
    3  *  @brief Provides functions for linear algebra operations on psImages and psVectors. 
     3 *  @brief Provides functions for linear algebra operations on psImages and psVectors.
    44 *
    55 *  Functions are provided to:
    … …  
    1919 *
    2020 *  @author Ross Harman, MHPCC
    21  *   
    22  *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
    23  *  @date $Date: 2004-11-04 01:04:59 $
     21 *
     22 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
     23 *  @date $Date: 2004-12-09 20:51:16 $
    2424 *
    2525 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4141#include "psMatrix.h"
    4242#include "psConstants.h"
    43 
    44 
     43#include "psDataManipErrors.h"
     44
     45
     46/*****************************************************************************/
     47/* DEFINE STATEMENTS                                                         */
     48/*****************************************************************************/
    4549
    4650/** Preprocessor macro to generate error for image dimensionality not set to PS_DIMEN_IMAGE */
    … …  
    5054            "Invalid operation. %s has incorrect dimensionality %d.", #NAME, PS_DIMEN);             \
    5155    return RETURN;                                                                                  \
    52 } else if(NAME->type.type != PS_TYPE_F64) {                                                         \
     56} else if(NAME->type.type!=PS_TYPE_F64 && NAME->type.type!=PS_TYPE_F32) {                           \
    5357    psError(PS_ERR_BAD_PARAMETER_TYPE, true,                                                        \
    5458            "Invalid operation. %s not PS_TYPE_F64.", #NAME);                                       \
    … …  
    6771#define PS_CHECK_SQUARE(NAME, RETURN)                                                               \
    6872if (NAME->numCols != NAME->numRows) {                                                               \
    69     psError(PS_ERR_BAD_PARAMETER_SIZE, true,                                                       \
    70             "Invalid operation: %s not square array.", #NAME);                             \
     73    psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid operation: %s not square array.", #NAME);     \
    7174    return RETURN;                                                                                  \
    7275}
    … …  
    7982LHS_NAME.data  = RHS_NAME;
    8083
     84
     85/*****************************************************************************/
     86/* FILE STATIC FUNCTIONS                                                     */
     87/*****************************************************************************/
     88
     89static void  psVectorToGslVector(gsl_vector *outGslVector, const psVector *inVector);
     90static void gslVectorToPsVector(psVector *outVector, gsl_vector *inGslVector);
     91static void  psImageToGslMatrix(gsl_matrix *outGslMatrix, const psImage *inImage);
     92static void gslMatrixToPsImage(psImage *outImage, gsl_matrix *inGslMatrix);
     93
     94/** Static function to copy psF32 or psF64 vector data to a GSL vector */
     95static void  psVectorToGslVector(gsl_vector *outGslVector, const psVector *inVector)
     96{
     97    psU32 i = 0;
     98    psU32 n = 0;
     99
     100
     101    n = inVector->n;
     102    for(i=0; i<n; i++) {
     103        if(inVector->type.type == PS_TYPE_F32) {
     104            outGslVector->data[i] = (psS32)inVector->data.F32[i];
     105        } else {
     106            outGslVector->data[i] = inVector->data.F64[i];
     107        }
     108    }
     109}
     110
     111/** Static function to copy GSL vector data to a psF32 or psF64 vector */
     112static void gslVectorToPsVector(psVector *outVector, gsl_vector *inGslVector)
     113{
     114    psU32 i = 0;
     115    psU32 n = 0;
     116
     117
     118    n = outVector->n;
     119    for(i=0; i<n; i++) {
     120        if(outVector->type.type == PS_TYPE_F32) {
     121            outVector->data.F32[i] = (psF32)inGslVector->data[i];
     122        } else {
     123            outVector->data.F64[i] = inGslVector->data[i];
     124        }
     125    }
     126}
     127
     128/** Static function to copy psF32 or psF64 image data to a GSL matrix */
     129static void  psImageToGslMatrix(gsl_matrix *outGslMatrix, const psImage *inImage)
     130{
     131    psU32 i = 0;
     132    psU32 j = 0;
     133    psU32 numRows = 0;
     134    psU32 numCols = 0;
     135
     136
     137    numRows = inImage->numRows;
     138    numCols = inImage->numCols;
     139    for(i=0; i<numRows; i++) {
     140        for(j=0; j<numCols; j++) {
     141            if(inImage->type.type == PS_TYPE_F32) {
     142                outGslMatrix->data[i*numCols+j] = (psS32)inImage->data.F32[i][j];
     143            } else {
     144                outGslMatrix->data[i*numCols+j] = inImage->data.F64[i][j];
     145            }
     146        }
     147    }
     148}
     149
     150/** Static function to copy GSL matrix data to a psF32 or psF64 image */
     151static void gslMatrixToPsImage(psImage *outImage, gsl_matrix *inGslMatrix)
     152{
     153    psU32 i = 0;
     154    psU32 j = 0;
     155    psU32 numRows = 0;
     156    psU32 numCols = 0;
     157
     158
     159    numRows = outImage->numRows;
     160    numCols = outImage->numCols;
     161    for(i=0; i<numRows; i++) {
     162        for(j=0; j<numCols; j++) {
     163            if(outImage->type.type == PS_TYPE_F32) {
     164                outImage->data.F32[i][j] = (psF32)inGslMatrix->data[i*numCols+j];
     165            } else {
     166                outImage->data.F64[i][j] = inGslMatrix->data[i*numCols+j];
     167            }
     168        }
     169    }
     170}
     171
     172
    81173/*****************************************************************************/
    82174/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    … …  
    86178{
    87179    psS32 signum = 0;
    88     psS32 arraySize = 0;
    89180    psS32 numRows = 0;
    90181    psS32 numCols = 0;
    91     gsl_matrix lu;
     182    gsl_matrix *lu = NULL;
    92183    gsl_permutation perm;
    93184
    94     #define psMatrixLUD_EXIT { \
    95                                psFree(outImage); \
    96                                return NULL; \
    97                              }
    98 
    99     // Error checks
     185
     186    #define psMatrixLUD_EXIT {psFree(outImage); return NULL;}
     187
     188    // Error checks
     189    PS_IMAGE_CHECK_NULL_GENERAL(inImage, psMatrixLUD_EXIT);
    100190    PS_CHECK_POINTERS(inImage, outImage, outImage);
    101191    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
    102 
    103     PS_IMAGE_CHECK_NULL_GENERAL(inImage, psMatrixLUD_EXIT);
    104192    PS_VECTOR_CHECK_NULL_GENERAL(outPerm, psMatrixLUD_EXIT);
    105 
    106193    PS_CHECK_DIMEN_AND_TYPE(outPerm, PS_DIMEN_VECTOR, outImage);
    107 
    108194    psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
    109195    psVectorRecycle(outPerm, inImage->numRows, inImage->type.type);
     196    PS_CHECK_SQUARE(inImage, outImage);
     197    PS_CHECK_SQUARE(outImage, outImage);
    110198
    111199    // Initialize data
    112200    numRows = inImage->numRows;
    113201    numCols = inImage->numCols;
    114     arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
    115202
    116203    // Initialize GSL data
    … …  
    118205    outPerm->n = numCols;
    119206    perm.data = outPerm->data.V;
    120     PS_GSL_MATRIX_INITIALIZE(lu, outImage->data.V[0]);
    121 
    122     // Non-square matrices not allowed
    123     PS_CHECK_SQUARE(inImage, outImage);
    124     PS_CHECK_SQUARE(outImage, outImage);
    125 
    126     // Copy psImage input data into GSL matrix data to keep input data pristine
    127     memcpy(lu.data, inImage->data.V[0], arraySize);
     207    lu = gsl_matrix_alloc(numRows, numCols);
     208
     209    // Copy psImage data into GSL matrix data
     210    psImageToGslMatrix(lu, inImage);
    128211
    129212    // Calculate LU decomposition
    130     gsl_linalg_LU_decomp(&lu, &perm, &signum);
     213    gsl_linalg_LU_decomp(lu, &perm, &signum);
     214
     215    // Copy GSL matrix data to psImage data
     216    gslMatrixToPsImage(outImage, lu);
     217
     218    // Free GSL data
     219    gsl_matrix_free(lu);
    131220
    132221    return outImage;
    133222}
    134223
    135 psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage,
    136                           const psVector* inVector, const psVector* inPerm)
    137 {
    138     psS32 arraySize = 0;
     224psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage, const psVector* inVector,
     225                          const psVector* inPerm)
     226{
    139227    psS32 numRows = 0;
    140228    psS32 numCols = 0;
    141     gsl_matrix lu;
     229    gsl_matrix *lu;
    142230    gsl_permutation perm;
    143     gsl_vector b;
    144     gsl_vector x;
     231    gsl_vector *b = NULL;
     232    gsl_vector *x = NULL;
    145233
    146234    // Error checks
    … …  
    157245    PS_VECTOR_CHECK_NULL(inPerm, outVector);
    158246    PS_CHECK_DIMEN_AND_TYPE(inPerm, PS_DIMEN_VECTOR, outVector);
    159 
    160247    psVectorRecycle(outVector, inImage->numRows, inImage->type.type);
    161 
    162248
    163249    // Initialize data
    164250    numRows = inImage->numRows;
    165251    numCols = inImage->numCols;
    166     arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
    167252
    168253    // Initialize GSL data
    169     PS_GSL_MATRIX_INITIALIZE(lu, inImage->data.V[0]);
     254    lu = gsl_matrix_alloc(numRows, numCols);
     255    psImageToGslMatrix(lu, inImage);
     256    b = gsl_vector_alloc(inVector->n);
     257    psVectorToGslVector(b, inVector);
     258    x = gsl_vector_alloc(inVector->n);
    170259
    171260    outVector->n = numCols;
    172 
    173261    perm.size = inPerm->n;
    174262    perm.data = inPerm->data.V;
    175263
    176     b.size = inVector->n;
    177     b.stride = 1;
    178     b.data = inVector->data.V;
    179 
    180     x.size = numCols;
    181     x.stride = 1;
    182     x.data = outVector->data.V;
    183 
    184264    // Solve for {x} in equation: {b} = [A]{x}
    185     gsl_linalg_LU_solve(&lu, &perm, &b, &x);
     265    gsl_linalg_LU_solve(lu, &perm, b, x);
     266
     267    // Copy GSL vector data to psVector data
     268    gslVectorToPsVector(outVector, x);
     269
     270    // Free GSL data
     271    gsl_vector_free(b);
     272    gsl_vector_free(x);
     273    gsl_matrix_free(lu);
    186274
    187275    return outVector;
    188276}
    189277
    190 psImage* psMatrixInvert(psImage* outImage, const psImage* inImage, float *restrict det)
     278psImage* psMatrixInvert(psImage* outImage, const psImage* inImage, psF32 *det)
    191279{
    192280    psS32 signum = 0;
    193     psS32 arraySize = 0;
    194281    psS32 numRows = 0;
    195282    psS32 numCols = 0;
    196     gsl_matrix inv;
    197     gsl_matrix *lu;
     283    gsl_matrix *inv = NULL;
     284    gsl_matrix *lu = NULL;
    198285    gsl_permutation *perm = NULL;
    199286
    … …  
    204291    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
    205292    PS_IMAGE_CHECK_EMPTY(inImage, outImage);
    206 
    207293    psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     294    PS_CHECK_SQUARE(inImage, outImage);
     295    PS_CHECK_SQUARE(outImage, outImage);
    208296
    209297    // Initialize data
    210298    numRows = inImage->numRows;
    211299    numCols = inImage->numCols;
    212     arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
    213 
    214     // Allocate GSL structs
    215     perm = gsl_permutation_alloc(inImage->numRows);
     300
     301    // Initialize GSL data
     302    perm = gsl_permutation_alloc(numRows);
    216303    lu = gsl_matrix_alloc(numRows, numCols);
    217 
    218     // Initialize GSL data
    219     PS_GSL_MATRIX_INITIALIZE(inv, outImage->data.V[0]);
    220 
    221     // Non-square matrices not allowed
    222     PS_CHECK_SQUARE(inImage, outImage);
    223     PS_CHECK_SQUARE(outImage, outImage);
    224 
    225     // Copy psImage input data into GSL matrix data to keep input data pristine
    226     memcpy(lu->data, inImage->data.V[0], arraySize);
     304    inv = gsl_matrix_alloc(numRows, numCols);
     305    psImageToGslMatrix(lu, inImage);
    227306
    228307    // Invert data and calculate determinant
    229308    gsl_linalg_LU_decomp(lu, perm, &signum);
    230     gsl_linalg_LU_invert(lu, perm, &inv);
     309    gsl_linalg_LU_invert(lu, perm, inv);
    231310    *det = (float)gsl_linalg_LU_det(lu, signum);
     311
     312    // Copy GSL matrix data to psImage data
     313    gslMatrixToPsImage(outImage, inv);
    232314
    233315    // Free GSL structs
    234316    gsl_permutation_free(perm);
    235317    gsl_matrix_free(lu);
     318    gsl_matrix_free(inv);
    236319
    237320    return outImage;
    238321}
    239322
    240 float *psMatrixDeterminant(const psImage* restrict inImage)
     323psF32 *psMatrixDeterminant(const psImage* inImage)
    241324{
    242325    psS32 signum = 0;
    243     psS32 arraySize = 0;
    244326    psS32 numRows = 0;
    245327    psS32 numCols = 0;
    246     float *det = NULL;
     328    psF32 *det = NULL;
    247329    gsl_matrix *lu = NULL;
    248330    gsl_permutation *perm = NULL;
    … …  
    252334    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, NULL);
    253335    PS_IMAGE_CHECK_EMPTY(inImage, NULL);
     336    PS_CHECK_SQUARE(inImage, 0);
    254337
    255338    // Initialize data
    256339    numRows = inImage->numRows;
    257340    numCols = inImage->numCols;
    258     arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
    259341
    260342    // Allocate GSL structs
    261343    perm = gsl_permutation_alloc(numRows);
    262344    lu = gsl_matrix_alloc(numRows, numCols);
    263 
    264     // Non-square matrices not allowed
    265     PS_CHECK_SQUARE(inImage, 0);
    266 
    267     // Copy psImage input data into GSL matrix data to keep input data pristine
    268     memcpy(lu->data, inImage->data.V[0], arraySize);
     345    psImageToGslMatrix(lu, inImage);
    269346
    270347    // Calculate determinant
    271     det = (float *)psAlloc(sizeof(float));
     348    det = (psF32*)psAlloc(sizeof(psF32));
    272349    gsl_linalg_LU_decomp(lu, perm, &signum);
    273     *det = (float)gsl_linalg_LU_det(lu, signum);
     350    *det = (psF32)gsl_linalg_LU_det(lu, signum);
    274351
    275352    // Free GSL structs
    … …  
    282359psImage* psMatrixMultiply(psImage* outImage, psImage* inImage1, psImage* inImage2)
    283360{
    284     psS32 arraySize = 0;
    285361    psS32 numRows = 0;
    286362    psS32 numCols = 0;
    287     gsl_matrix m1;
    288     gsl_matrix m2;
    289     gsl_matrix m3;
     363    gsl_matrix *m1 = NULL;
     364    gsl_matrix *m2 = NULL;
     365    gsl_matrix *m3 = NULL;
    290366
    291367    // Error checks
    … …  
    299375    PS_IMAGE_CHECK_EMPTY(inImage2, outImage);
    300376    PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, outImage);
    301 
    302377    psImageRecycle(outImage, inImage2->numCols, inImage2->numRows, inImage2->type.type);
     378    PS_CHECK_SQUARE(inImage1, outImage);
     379    PS_CHECK_SQUARE(inImage2, outImage);
     380    PS_CHECK_SQUARE(outImage, outImage);
    303381
    304382    // Initialize data
    305383    numRows = inImage1->numRows;
    306384    numCols = inImage1->numCols;
    307     arraySize = PSELEMTYPE_SIZEOF(outImage->type.type) * numRows * numCols;
    308385
    309386    // Initialize GSL data
    310     PS_GSL_MATRIX_INITIALIZE(m1, inImage1->data.V[0]);
    311     PS_GSL_MATRIX_INITIALIZE(m2, inImage2->data.V[0]);
    312     PS_GSL_MATRIX_INITIALIZE(m3, outImage->data.V[0]);
    313 
    314     // Non-square matrices not allowed
    315     PS_CHECK_SQUARE(inImage1, outImage);
    316     PS_CHECK_SQUARE(inImage2, outImage);
    317     PS_CHECK_SQUARE(outImage, outImage);
    318 
    319     gsl_linalg_matmult(&m1, &m2, &m3);
     387    m1 = gsl_matrix_alloc(numRows, numCols);
     388    psImageToGslMatrix(m1, inImage1);
     389    m2 = gsl_matrix_alloc(numRows, numCols);
     390    psImageToGslMatrix(m2, inImage2);
     391    m3 = gsl_matrix_alloc(numRows, numCols);
     392    psImageToGslMatrix(m3, outImage);
     393
     394    // Perform multiplication
     395    gsl_linalg_matmult(m1, m2, m3);
     396
     397    // Copy GSL matrix data to psImage data
     398    gslMatrixToPsImage(outImage, m3);
     399
     400    // Free GSL structs
     401    gsl_matrix_free(m1);
     402    gsl_matrix_free(m2);
     403    gsl_matrix_free(m3);
    320404
    321405    return outImage;
    … …  
    324408psImage* psMatrixTranspose(psImage* outImage, const psImage* inImage)
    325409{
    326     psS32 arraySize = 0;
    327     psS32 numRows = 0;
    328     psS32 numCols = 0;
    329     gsl_matrix trans;
     410    psU32 i = 0;
     411    psU32 j = 0;
     412    psS32 numRowsIn = 0;
     413    psS32 numColsIn = 0;
     414    psS32 numRowsOut = 0;
     415    psS32 numColsOut = 0;
     416
    330417
    331418    // Error checks
    … …  
    334421    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
    335422    PS_IMAGE_CHECK_EMPTY(inImage, outImage);
    336 
    337     psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     423    outImage = psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
    338424
    339425    // Initialize data
    340     numRows = inImage->numRows;
    341     numCols = inImage->numCols;
    342     arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
    343 
    344     // Initialize GSL data
    345     PS_GSL_MATRIX_INITIALIZE(trans, outImage->data.V[0]);
    346 
    347     // Non-square matrices not allowed
    348     PS_CHECK_SQUARE(inImage, outImage);
    349     PS_CHECK_SQUARE(outImage, outImage);
    350 
    351     // Copy psImage input data into psImage output data to keep input data pristine
    352     memcpy(outImage->data.V[0], inImage->data.V[0], arraySize);
    353 
    354     // Transpose data
    355     gsl_matrix_transpose(&trans);
     426    numRowsIn = inImage->numRows;
     427    numColsIn = inImage->numCols;
     428    numRowsOut = outImage->numRows;
     429    numColsOut = outImage->numCols;
     430
     431    if(numRowsIn!=numColsOut && numRowsOut!=numColsIn) {
     432        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psMatrix_TRANSPOSE_MISMATCH);
     433    }
     434
     435    if(outImage->type.type == PS_TYPE_F32) {
     436        for(i=0; i<numRowsOut; i++) {
     437            for(j=0; j<numColsOut; j++) {
     438                outImage->data.F32[i][j] = inImage->data.F32[j][i];
     439            }
     440        }
     441    } else {
     442        for(i=0; i<numRowsOut; i++) {
     443            for(j=0; j<numColsOut; j++) {
     444                outImage->data.F64[i][j] = inImage->data.F64[j][i];
     445            }
     446        }
     447    }
    356448
    357449    return outImage;
    … …  
    364456    gsl_vector *eVals = NULL;
    365457    gsl_eigen_symmv_workspace *w = NULL;
    366     gsl_matrix out;
    367     gsl_matrix in;
     458    gsl_matrix *out = NULL;
     459    gsl_matrix *in = NULL;
     460
    368461
    369462    // Error checks
    … …  
    372465    PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
    373466    PS_IMAGE_CHECK_EMPTY(inImage, outImage);
    374 
    375467    psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
    376468
    … …  
    379471    numCols = inImage->numCols;
    380472
    381     // Initialize GSL data
    382     PS_GSL_MATRIX_INITIALIZE(in, inImage->data.V[0]);
    383     PS_GSL_MATRIX_INITIALIZE(out, outImage->data.V[0]);
     473    in = gsl_matrix_alloc(numRows, numCols);
     474    psImageToGslMatrix(in, inImage);
     475    out = gsl_matrix_alloc(numRows, numCols);
    384476
    385477    // Allocate GSL structs
    … …  
    392484
    393485    // Calculate Eigenvalues and Eigenvectors...Eigenvalues not currently used
    394     gsl_eigen_symmv(&in, eVals, &out, w);
     486    gsl_eigen_symmv(in, eVals, out, w);
     487
     488    // Copy GSL matrix data to psImage data
     489    gslMatrixToPsImage(outImage, out);
    395490
    396491    // Free GSL structs
     492    gsl_matrix_free(in);
     493    gsl_matrix_free(out);
    397494    gsl_eigen_symmv_free(w);
    398495    gsl_vector_free(eVals);
    … …  
    405502    psS32 size = 0;
    406503
    407     #define psMatrixToVector_EXIT { \
    408                                     psFree(outVector); \
    409                                     return NULL; \
    410                                   }
     504    #define psMatrixToVector_EXIT {psFree(outVector); return NULL;}
    411505
    412506    // Error checks
Note: See TracChangeset for help on using the changeset viewer.