IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 29, 2004, 6:03:51 PM (22 years ago)
Author:
evanalst
Message:

Add check for input type psF32 if forward FFT.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psFFT.c

    r1264 r1349  
    11/** @file  psFFT.c
    2  *
    3  *  @brief Contains FFT transforms functions
    4  *
    5  *  @author Robert DeSonia, MHPCC
    6  *
    7  *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2004-07-22 20:42:54 $
    9  *
    10  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    11  */
     2*
     3*  @brief Contains FFT transforms functions
     4*
     5*  @author Robert DeSonia, MHPCC
     6*
     7*  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2004-07-30 04:03:51 $
     9*
     10*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     11*/
    1212
    1313#include <unistd.h>
    … …  
    2727static bool p_fftwWisdomImported = false;
    2828
    29 psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction)
     29psImage* psImageFFT( psImage* out, const psImage* in, psFftDirection direction )
    3030{
    3131    unsigned int numCols;
    … …  
    3333    psElemType type;
    3434    fftwf_plan plan;
    35 
     35   
    3636    /* got good image data? */
    37     if (in==NULL) {
    38         psFree(out);
    39         return NULL;
    40     }
    41 
    42     type = in->type.type;
    43 
    44     if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) {
    45         psError(__func__,"Input image must be a 32-bit float or complex image (type=%d)",
    46                 type);
    47         psFree(out);
    48         return NULL;
    49     }
    50 
     37    if ( in == NULL ) {
     38            psFree( out );
     39            return NULL;
     40        }
     41       
     42    type = in->type.type;
     43   
     44    if ( ( type != PS_TYPE_F32 ) && ( type != PS_TYPE_C32 ) ) {
     45            psError( __func__, "Input image must be a 32-bit float or complex image (type=%d)",
     46                     type );
     47            psFree( out );
     48            return NULL;
     49        }
     50       
    5151    if ( type != PS_TYPE_C32 && direction == PS_FFT_REVERSE ) {
    52         psError(__func__,"Input image must be complex image for reverse FFT (type=%d).",
    53                 type);
    54         psFree(out);
    55         return NULL;
    56 
    57     }
    58 
     52            psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
     53                     type );
     54            psFree( out );
     55            return NULL;
     56           
     57        }
     58       
    5959    /* make sure the system-level wisdom information is imported. */
    60     if (! p_fftwWisdomImported) {
    61         fftwf_import_system_wisdom();
    62         p_fftwWisdomImported = true;
    63     }
    64 
     60    if ( ! p_fftwWisdomImported ) {
     61            fftwf_import_system_wisdom();
     62            p_fftwWisdomImported = true;
     63        }
     64       
    6565    numRows = in->numRows;
    6666    numCols = in->numCols;
    67 
    68     out = psImageCopy(out,in, PS_TYPE_C32);
    69 
    70     plan = fftwf_plan_dft_2d(numCols, numRows,
    71                              (fftwf_complex*)out->data.C32[0],
    72                              (fftwf_complex*)out->data.C32[0],
    73                              direction,
    74                              P_FFTW_PLAN_RIGOR);
    75 
     67   
     68    out = psImageCopy( out, in, PS_TYPE_C32 );
     69   
     70    plan = fftwf_plan_dft_2d( numCols, numRows,
     71                              ( fftwf_complex* ) out->data.C32[ 0 ],
     72                              ( fftwf_complex* ) out->data.C32[ 0 ],
     73                              direction,
     74                              P_FFTW_PLAN_RIGOR );
     75                             
    7676    /* check if a plan exists now*/
    77     if (plan == NULL) {
    78         psError(__func__,"Failed to create FFTW plan.");
    79         psFree(out);
    80         return NULL;
    81     }
    82 
     77    if ( plan == NULL ) {
     78            psError( __func__, "Failed to create FFTW plan." );
     79            psFree( out );
     80            return NULL;
     81        }
     82       
    8383    /* finally, call FFTW with the plan made above */
    84     fftwf_execute(plan);
    85 
    86     fftwf_destroy_plan(plan);
    87 
    88     return out;
    89 
    90 }
    91 
    92 
    93 psImage *psImageReal(psImage *out, const psImage* in)
     84    fftwf_execute( plan );
     85   
     86    fftwf_destroy_plan( plan );
     87   
     88    return out;
     89   
     90}
     91
     92
     93psImage *psImageReal( psImage *out, const psImage* in )
    9494{
    9595    psElemType type;
    9696    unsigned int numCols;
    9797    unsigned int numRows;
    98 
    99 
    100     if (in == NULL) {
    101         psFree(out);
    102         return NULL;
    103     }
    104 
     98   
     99   
     100    if ( in == NULL ) {
     101            psFree( out );
     102            return NULL;
     103        }
     104       
    105105    type = in->type.type;
    106106    numCols = in->numCols;
    107107    numRows = in->numRows;
    108 
     108   
    109109    /* if not a complex number, this is logically just a copy */
    110     if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    111         // Warn user, as this is probably not expected
    112         psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. "
    113                  "Just an image copy was performed.");
    114         return psImageCopy(out,in,type);
    115     }
    116 
    117     if (type == PS_TYPE_C32) {
    118         psF32* outRow;
    119         psC32* inRow;
    120 
    121         out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
    122         for (unsigned int row=0;row<numRows;row++) {
    123             outRow = out->data.F32[row];
    124             inRow = in->data.C32[row];
    125 
    126             for (unsigned int col=0;col<numCols;col++) {
    127                 outRow[col] = crealf(inRow[col]);
    128             }
    129         }
    130     } else if (type == PS_TYPE_C64) {
    131         psF64* outRow;
    132         psC64* inRow;
    133 
    134         out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
    135         for (unsigned int row=0;row<numRows;row++) {
    136             outRow = out->data.F64[row];
    137             inRow = in->data.C64[row];
    138 
    139             for (unsigned int col=0;col<numCols;col++) {
    140                 outRow[col] = creal(inRow[col]);
    141             }
    142         }
    143     } else {
    144         psError(__func__,"Can not extract real component from given image type (%d).",
    145                 type);
    146         psFree(out);
    147         return NULL;
    148     }
    149 
    150     return out;
    151 }
    152 
    153 psImage *psImageImaginary(psImage *out, const psImage* in)
     110    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     111            // Warn user, as this is probably not expected
     112            psLogMsg( __func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
     113                      "Just an image copy was performed." );
     114            return psImageCopy( out, in, type );
     115        }
     116       
     117    if ( type == PS_TYPE_C32 ) {
     118            psF32 * outRow;
     119            psC32* inRow;
     120           
     121            out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
     122            for ( unsigned int row = 0;row < numRows;row++ ) {
     123                    outRow = out->data.F32[ row ];
     124                    inRow = in->data.C32[ row ];
     125                   
     126                    for ( unsigned int col = 0;col < numCols;col++ ) {
     127                            outRow[ col ] = crealf( inRow[ col ] );
     128                        }
     129                }
     130        } else if ( type == PS_TYPE_C64 ) {
     131            psF64 * outRow;
     132            psC64* inRow;
     133           
     134            out = psImageRecycle( out, numCols, numRows, PS_TYPE_F64 );
     135            for ( unsigned int row = 0;row < numRows;row++ ) {
     136                    outRow = out->data.F64[ row ];
     137                    inRow = in->data.C64[ row ];
     138                   
     139                    for ( unsigned int col = 0;col < numCols;col++ ) {
     140                            outRow[ col ] = creal( inRow[ col ] );
     141                        }
     142                }
     143        } else {
     144            psError( __func__, "Can not extract real component from given image type (%d).",
     145                     type );
     146            psFree( out );
     147            return NULL;
     148        }
     149       
     150    return out;
     151}
     152
     153psImage *psImageImaginary( psImage *out, const psImage* in )
    154154{
    155155    psElemType type;
    156156    unsigned int numCols;
    157157    unsigned int numRows;
    158 
    159 
    160     if (in == NULL) {
    161         psFree(out);
    162         return NULL;
    163     }
    164 
     158   
     159   
     160    if ( in == NULL ) {
     161            psFree( out );
     162            return NULL;
     163        }
     164       
    165165    type = in->type.type;
    166166    numCols = in->numCols;
    167167    numRows = in->numRows;
    168 
     168   
    169169    /* if not a complex number, this is logically just zeroed image of same size */
    170     if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    171         // Warn user, as this is probably not expected
    172         psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. "
    173                  "A zero image was returned.");
    174         out = psImageRecycle(out,numCols,numRows,type);
    175         memset(out->data.V[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows);
    176         return out;
    177     }
    178 
    179     if (type == PS_TYPE_C32) {
    180         psF32* outRow;
    181         psC32* inRow;
    182 
    183         out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
    184         for (unsigned int row=0;row<numRows;row++) {
    185             outRow = out->data.F32[row];
    186             inRow = in->data.C32[row];
    187 
    188             for (unsigned int col=0;col<numCols;col++) {
    189                 outRow[col] = cimagf(inRow[col]);
    190             }
    191         }
    192     } else if (type == PS_TYPE_C64) {
    193         psF64* outRow;
    194         psC64* inRow;
    195 
    196         out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
    197         for (unsigned int row=0;row<numRows;row++) {
    198             outRow = out->data.F64[row];
    199             inRow = in->data.C64[row];
    200 
    201             for (unsigned int col=0;col<numCols;col++) {
    202                 outRow[col] = cimag(inRow[col]);
    203             }
    204         }
    205     } else {
    206         psError(__func__,"Can not extract imaginary component from given image type (%d).",
    207                 type);
    208         psFree(out);
    209         return NULL;
    210     }
    211 
    212     return out;
    213 }
    214 
    215 psImage *psImageComplex(psImage* out, psImage *real, const psImage *imag)
     170    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     171            // Warn user, as this is probably not expected
     172            psLogMsg( __func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
     173                      "A zero image was returned." );
     174            out = psImageRecycle( out, numCols, numRows, type );
     175            memset( out->data.V[ 0 ], 0, PSELEMTYPE_SIZEOF( type ) * numCols * numRows );
     176            return out;
     177        }
     178       
     179    if ( type == PS_TYPE_C32 ) {
     180            psF32 * outRow;
     181            psC32* inRow;
     182           
     183            out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
     184            for ( unsigned int row = 0;row < numRows;row++ ) {
     185                    outRow = out->data.F32[ row ];
     186                    inRow = in->data.C32[ row ];
     187                   
     188                    for ( unsigned int col = 0;col < numCols;col++ ) {
     189                            outRow[ col ] = cimagf( inRow[ col ] );
     190                        }
     191                }
     192        } else if ( type == PS_TYPE_C64 ) {
     193            psF64 * outRow;
     194            psC64* inRow;
     195           
     196            out = psImageRecycle( out, numCols, numRows, PS_TYPE_F64 );
     197            for ( unsigned int row = 0;row < numRows;row++ ) {
     198                    outRow = out->data.F64[ row ];
     199                    inRow = in->data.C64[ row ];
     200                   
     201                    for ( unsigned int col = 0;col < numCols;col++ ) {
     202                            outRow[ col ] = cimag( inRow[ col ] );
     203                        }
     204                }
     205        } else {
     206            psError( __func__, "Can not extract imaginary component from given image type (%d).",
     207                     type );
     208            psFree( out );
     209            return NULL;
     210        }
     211       
     212    return out;
     213}
     214
     215psImage *psImageComplex( psImage* out, psImage *real, const psImage *imag )
    216216{
    217217    psElemType type;
    218218    unsigned int numCols;
    219219    unsigned int numRows;
    220 
    221 
    222     if (real == NULL || imag == NULL) {
    223         psFree(out);
    224         return NULL;
    225     }
    226 
     220   
     221   
     222    if ( real == NULL || imag == NULL ) {
     223            psFree( out );
     224            return NULL;
     225        }
     226       
    227227    type = real->type.type;
    228228    numCols = real->numCols;
    229229    numRows = real->numRows;
    230 
    231     if (imag->type.type != type) {
    232         psError(__func__,"The inputs to psImageComplex must be the same type.");
    233         psFree(out);
    234         return NULL;
    235     }
    236 
    237     if (imag->numCols != numCols ||
    238             imag->numRows != numRows) {
    239         psError(__func__,"The inputs to psImageComplex must be the same dimensions.");
    240         psFree(out);
    241         return NULL;
    242     }
    243 
    244     if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
    245         psError(__func__,"The inputs to psImageComplex can not be complex.");
    246         psFree(out);
    247         return NULL;
    248     }
    249 
    250     if (type != PS_TYPE_F32 && type != PS_TYPE_F64) {
    251         psError(__func__,"The input type to psImageComplex must be a floating point.");
    252         psFree(out);
    253         return NULL;
    254     }
    255 
    256     if (type == PS_TYPE_F32) {
    257         psC32* outRow;
    258         psF32* realRow;
    259         psF32* imagRow;
    260 
    261         out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
    262 
    263         for (unsigned int row=0;row<numRows;row++) {
    264             outRow = out->data.C32[row];
    265             realRow = real->data.F32[row];
    266             imagRow = imag->data.F32[row];
    267 
    268             for (unsigned int col=0;col<numCols;col++) {
    269                 outRow[col] = realRow[col] + I*imagRow[col];
    270             }
    271         }
    272     } else if (type == PS_TYPE_F64) {
    273         psC64* outRow;
    274         psF64* realRow;
    275         psF64* imagRow;
    276 
    277         out = psImageRecycle(out, numCols,numRows,PS_TYPE_C64);
    278         for (unsigned int row=0;row<numRows;row++) {
    279             outRow = out->data.C64[row];
    280             realRow = real->data.F64[row];
    281             imagRow = imag->data.F64[row];
    282 
    283             for (unsigned int col=0;col<numCols;col++) {
    284                 outRow[col] = realRow[col] + I*imagRow[col];
    285             }
    286         }
    287     } else {
    288         psError(__func__,"Can not merge real and imaginary portions for given image type (%d).",
    289                 type);
    290         psFree(out);
    291         return NULL;
    292     }
    293 
    294     return out;
    295 }
    296 
    297 psImage *psImageConjugate(psImage *out, const psImage *in)
     230   
     231    if ( imag->type.type != type ) {
     232            psError( __func__, "The inputs to psImageComplex must be the same type." );
     233            psFree( out );
     234            return NULL;
     235        }
     236       
     237    if ( imag->numCols != numCols ||
     238            imag->numRows != numRows ) {
     239            psError( __func__, "The inputs to psImageComplex must be the same dimensions." );
     240            psFree( out );
     241            return NULL;
     242        }
     243       
     244    if ( PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     245            psError( __func__, "The inputs to psImageComplex can not be complex." );
     246            psFree( out );
     247            return NULL;
     248        }
     249       
     250    if ( type != PS_TYPE_F32 && type != PS_TYPE_F64 ) {
     251            psError( __func__, "The input type to psImageComplex must be a floating point." );
     252            psFree( out );
     253            return NULL;
     254        }
     255       
     256    if ( type == PS_TYPE_F32 ) {
     257            psC32 * outRow;
     258            psF32* realRow;
     259            psF32* imagRow;
     260           
     261            out = psImageRecycle( out, numCols, numRows, PS_TYPE_C32 );
     262           
     263            for ( unsigned int row = 0;row < numRows;row++ ) {
     264                    outRow = out->data.C32[ row ];
     265                    realRow = real->data.F32[ row ];
     266                    imagRow = imag->data.F32[ row ];
     267                   
     268                    for ( unsigned int col = 0;col < numCols;col++ ) {
     269                            outRow[ col ] = realRow[ col ] + I * imagRow[ col ];
     270                        }
     271                }
     272        } else if ( type == PS_TYPE_F64 ) {
     273            psC64 * outRow;
     274            psF64* realRow;
     275            psF64* imagRow;
     276           
     277            out = psImageRecycle( out, numCols, numRows, PS_TYPE_C64 );
     278            for ( unsigned int row = 0;row < numRows;row++ ) {
     279                    outRow = out->data.C64[ row ];
     280                    realRow = real->data.F64[ row ];
     281                    imagRow = imag->data.F64[ row ];
     282                   
     283                    for ( unsigned int col = 0;col < numCols;col++ ) {
     284                            outRow[ col ] = realRow[ col ] + I * imagRow[ col ];
     285                        }
     286                }
     287        } else {
     288            psError( __func__, "Can not merge real and imaginary portions for given image type (%d).",
     289                     type );
     290            psFree( out );
     291            return NULL;
     292        }
     293       
     294    return out;
     295}
     296
     297psImage *psImageConjugate( psImage *out, const psImage *in )
    298298{
    299299    psElemType type;
    300300    unsigned int numCols;
    301301    unsigned int numRows;
    302 
    303 
    304     if (in == NULL) {
    305         psFree(out);
    306         return NULL;
    307     }
    308 
     302   
     303   
     304    if ( in == NULL ) {
     305            psFree( out );
     306            return NULL;
     307        }
     308       
    309309    type = in->type.type;
    310310    numCols = in->numCols;
    311311    numRows = in->numRows;
    312 
     312   
    313313    /* if not a complex number, this is logically just a image copy */
    314     if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    315         // Warn user, as this is probably not expected
    316         psLogMsg(__func__,PS_LOG_WARN,"Complex Conjugate of a non-Complex type called for. "
    317                  "Image copy was performed instead.");
    318         return psImageCopy(out,in,type);
    319     }
    320 
    321     if (type == PS_TYPE_C32) {
    322         psC32* outRow;
    323         psC32* inRow;
    324 
    325         out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
    326         for (unsigned int row=0;row<numRows;row++) {
    327             outRow = out->data.C32[row];
    328             inRow = in->data.C32[row];
    329 
    330             for (unsigned int col=0;col<numCols;col++) {
    331                 outRow[col] = crealf(inRow[col]) - I*cimagf(inRow[col]);
    332             }
    333         }
    334     } else if (type == PS_TYPE_C64) {
    335         psC64* outRow;
    336         psC64* inRow;
    337 
    338         out = psImageRecycle(out,numCols,numRows,PS_TYPE_C64);
    339         for (unsigned int row=0;row<numRows;row++) {
    340             outRow = out->data.C64[row];
    341             inRow = in->data.C64[row];
    342 
    343             for (unsigned int col=0;col<numCols;col++) {
    344                 outRow[col] = creal(inRow[col]) - I*cimag(inRow[col]);
    345             }
    346         }
    347     } else {
    348         psError(__func__,"Can not compute complex conjugate for given image type (%d).",
    349                 type);
    350         psFree(out);
    351         return NULL;
    352     }
    353 
    354     return out;
    355 }
    356 
    357 psImage *psImagePowerSpectrum(psImage* out, const psImage* in)
     314    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     315            // Warn user, as this is probably not expected
     316            psLogMsg( __func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
     317                      "Image copy was performed instead." );
     318            return psImageCopy( out, in, type );
     319        }
     320       
     321    if ( type == PS_TYPE_C32 ) {
     322            psC32 * outRow;
     323            psC32* inRow;
     324           
     325            out = psImageRecycle( out, numCols, numRows, PS_TYPE_C32 );
     326            for ( unsigned int row = 0;row < numRows;row++ ) {
     327                    outRow = out->data.C32[ row ];
     328                    inRow = in->data.C32[ row ];
     329                   
     330                    for ( unsigned int col = 0;col < numCols;col++ ) {
     331                            outRow[ col ] = crealf( inRow[ col ] ) - I * cimagf( inRow[ col ] );
     332                        }
     333                }
     334        } else if ( type == PS_TYPE_C64 ) {
     335            psC64 * outRow;
     336            psC64* inRow;
     337           
     338            out = psImageRecycle( out, numCols, numRows, PS_TYPE_C64 );
     339            for ( unsigned int row = 0;row < numRows;row++ ) {
     340                    outRow = out->data.C64[ row ];
     341                    inRow = in->data.C64[ row ];
     342                   
     343                    for ( unsigned int col = 0;col < numCols;col++ ) {
     344                            outRow[ col ] = creal( inRow[ col ] ) - I * cimag( inRow[ col ] );
     345                        }
     346                }
     347        } else {
     348            psError( __func__, "Can not compute complex conjugate for given image type (%d).",
     349                     type );
     350            psFree( out );
     351            return NULL;
     352        }
     353       
     354    return out;
     355}
     356
     357psImage *psImagePowerSpectrum( psImage* out, const psImage* in )
    358358{
    359359    psElemType type;
    … …  
    361361    unsigned int numRows;
    362362    int numElementsSquared;
    363 
    364     if (in == NULL) {
    365         psFree(out);
    366         return NULL;
    367     }
    368 
     363   
     364    if ( in == NULL ) {
     365            psFree( out );
     366            return NULL;
     367        }
     368       
    369369    type = in->type.type;
    370370    numCols = in->numCols;
    371371    numRows = in->numRows;
    372     numElementsSquared = numCols*numCols*numRows*numRows;
    373 
     372    numElementsSquared = numCols * numCols * numRows * numRows;
     373   
    374374    /* if not a complex number, this is not implemented */
    375     if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    376         psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
    377         psFree(out);
    378         return NULL;
    379     }
    380 
    381     if (type == PS_TYPE_C32) {
    382         psF32* outRow;
    383         psC32* inRow;
    384         psF32 real;
    385         psF32 imag;
    386 
    387 
    388         out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
    389         for (unsigned int row=0;row<numRows;row++) {
    390             outRow = out->data.F32[row];
    391             inRow = in->data.C32[row];
    392 
    393             for (unsigned int col=0;col<numCols;col++) {
    394                 real = crealf(inRow[col]);
    395                 imag = cimagf(inRow[col]);
    396                 outRow[col] = (real*real+imag*imag)/numElementsSquared;
    397             }
    398         }
    399     } else if (type == PS_TYPE_C64) {
    400         psF64* outRow;
    401         psC64* inRow;
    402         psF64 real;
    403         psF64 imag;
    404 
    405 
    406         out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
    407         for (unsigned int row=0;row<numRows;row++) {
    408             outRow = out->data.F64[row];
    409             inRow = in->data.C64[row];
    410 
    411             for (unsigned int col=0;col<numCols;col++) {
    412                 real = crealf(inRow[col]);
    413                 imag = cimagf(inRow[col]);
    414                 outRow[col] = real*real+imag*imag/numElementsSquared;
    415             }
    416         }
    417     } else {
    418         psError(__func__,"Can not power spectrum for given image type (%d).",
    419                 type);
    420         psFree(out);
    421         return NULL;
    422     }
    423 
    424     return out;
    425 
     375    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     376            psError( __func__, "Power Spectrum for non-complex inputs is not implemented." );
     377            psFree( out );
     378            return NULL;
     379        }
     380       
     381    if ( type == PS_TYPE_C32 ) {
     382            psF32 * outRow;
     383            psC32* inRow;
     384            psF32 real;
     385            psF32 imag;
     386           
     387           
     388            out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
     389            for ( unsigned int row = 0;row < numRows;row++ ) {
     390                    outRow = out->data.F32[ row ];
     391                    inRow = in->data.C32[ row ];
     392                   
     393                    for ( unsigned int col = 0;col < numCols;col++ ) {
     394                            real = crealf( inRow[ col ] );
     395                            imag = cimagf( inRow[ col ] );
     396                            outRow[ col ] = ( real * real + imag * imag ) / numElementsSquared;
     397                        }
     398                }
     399        } else if ( type == PS_TYPE_C64 ) {
     400            psF64 * outRow;
     401            psC64* inRow;
     402            psF64 real;
     403            psF64 imag;
     404           
     405           
     406            out = psImageRecycle( out, numCols, numRows, PS_TYPE_F64 );
     407            for ( unsigned int row = 0;row < numRows;row++ ) {
     408                    outRow = out->data.F64[ row ];
     409                    inRow = in->data.C64[ row ];
     410                   
     411                    for ( unsigned int col = 0;col < numCols;col++ ) {
     412                            real = crealf( inRow[ col ] );
     413                            imag = cimagf( inRow[ col ] );
     414                            outRow[ col ] = real * real + imag * imag / numElementsSquared;
     415                        }
     416                }
     417        } else {
     418            psError( __func__, "Can not power spectrum for given image type (%d).",
     419                     type );
     420            psFree( out );
     421            return NULL;
     422        }
     423       
     424    return out;
     425   
    426426}
    427427
    428428/************************************** Vector Functions ***************************************/
    429429
    430 psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction)
     430psVector* psVectorFFT( psVector* out, const psVector* in, psFftDirection direction )
    431431{
    432432    unsigned int numElements;
    433433    psElemType type;
    434434    fftwf_plan plan;
    435 
     435   
    436436    /* got good image data? */
    437     if (in==NULL) {
    438         psFree(out);
    439         return NULL;
    440     }
    441 
    442     type = in->type.type;
    443 
    444     if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) {
    445         psError(__func__,"Input image must be a 32-bit float or complex image (type=%d)",
    446                 type);
    447         psFree(out);
    448         return NULL;
    449     }
    450 
    451     if ( (type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE) ) {
    452         psError(__func__,"Input image must be complex image for reverse FFT (type=%d).",
    453                 type);
    454         psFree(out);
    455         return NULL;
    456 
    457     }
    458 
     437    if ( in == NULL ) {
     438            psFree( out );
     439            return NULL;
     440        }
     441       
     442    type = in->type.type;
     443   
     444    if ( ( type != PS_TYPE_F32 ) && ( type != PS_TYPE_C32 ) ) {
     445            psError( __func__, "Input image must be a 32-bit float or complex image (type=%d)",
     446                     type );
     447            psFree( out );
     448            return NULL;
     449        }
     450       
     451    if ( ( type != PS_TYPE_C32 ) && ( direction == PS_FFT_REVERSE ) ) {
     452            psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
     453                     type );
     454            psFree( out );
     455            return NULL;
     456           
     457        }
     458       
     459    if ( ( type != PS_TYPE_F32 ) && ( direction == PS_FFT_FORWARD ) ) {
     460            psError( __func__, "Input image must be real image for forward FFT (type=%d).",
     461                     type );
     462            psFree( out );
     463            return NULL;
     464        }
     465       
    459466    /* make sure the system-level wisdom information is imported. */
    460     if (! p_fftwWisdomImported) {
    461         fftwf_import_system_wisdom();
    462         p_fftwWisdomImported = true;
    463     }
    464 
     467    if ( ! p_fftwWisdomImported ) {
     468            fftwf_import_system_wisdom();
     469            p_fftwWisdomImported = true;
     470        }
     471       
    465472    numElements = in->n;
    466 
    467     out = psVectorRecycle(out, PS_TYPE_C32, numElements);
     473   
     474    out = psVectorRecycle( out, PS_TYPE_C32, numElements );
    468475    out->n = numElements;
    469 
    470     if (type == PS_TYPE_F32) {
    471         // need to convert to complex
    472         psC32* outVec = out->data.C32;
    473         psF32* inVec = in->data.F32;
    474         for (unsigned int i=0;i<numElements;i++) {
    475             outVec[i] = inVec[i];
    476         }
    477     } else {
    478         psC32* outVec = out->data.C32;
    479         psC32* inVec = in->data.C32;
    480         for (unsigned int i=0;i<numElements;i++) {
    481             outVec[i] = inVec[i];
    482         }
    483     }
    484 
    485     plan = fftwf_plan_dft_1d(numElements,
    486                              (fftwf_complex*)out->data.C32,
    487                              (fftwf_complex*)out->data.C32,
    488                              direction,
    489                              P_FFTW_PLAN_RIGOR);
    490 
     476   
     477    if ( type == PS_TYPE_F32 ) {
     478            // need to convert to complex
     479            psC32 * outVec = out->data.C32;
     480            psF32* inVec = in->data.F32;
     481            for ( unsigned int i = 0;i < numElements;i++ ) {
     482                    outVec[ i ] = inVec[ i ];
     483                }
     484        } else {
     485            psC32* outVec = out->data.C32;
     486            psC32* inVec = in->data.C32;
     487            for ( unsigned int i = 0;i < numElements;i++ ) {
     488                    outVec[ i ] = inVec[ i ];
     489                }
     490        }
     491       
     492    plan = fftwf_plan_dft_1d( numElements,
     493                              ( fftwf_complex* ) out->data.C32,
     494                              ( fftwf_complex* ) out->data.C32,
     495                              direction,
     496                              P_FFTW_PLAN_RIGOR );
     497                             
    491498    /* check if a plan exists now*/
    492     if (plan == NULL) {
    493         psError(__func__,"Failed to create FFTW plan.");
    494         psFree(out);
    495         return NULL;
    496     }
    497 
     499    if ( plan == NULL ) {
     500            psError( __func__, "Failed to create FFTW plan." );
     501            psFree( out );
     502            return NULL;
     503        }
     504       
    498505    /* finally, call FFTW with the plan made above */
    499     fftwf_execute(plan);
    500 
    501     fftwf_destroy_plan(plan);
    502 
    503     return out;
    504 }
    505 
    506 
    507 psVector *psVectorReal(psVector *out, const psVector* in)
     506    fftwf_execute( plan );
     507   
     508    fftwf_destroy_plan( plan );
     509   
     510    return out;
     511}
     512
     513
     514psVector *psVectorReal( psVector *out, const psVector* in )
    508515{
    509516    psElemType type;
    510517    unsigned int numElements;
    511 
    512     if (in == NULL) {
    513         psFree(out);
    514         return NULL;
    515     }
    516 
     518   
     519    if ( in == NULL ) {
     520            psFree( out );
     521            return NULL;
     522        }
     523       
    517524    type = in->type.type;
    518525    numElements = in->n;
    519 
     526   
    520527    /* if not a complex number, this is logically just a copy */
    521     if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    522         // Warn user, as this is probably not expected
    523         psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. "
    524                  "Just a vector copy was performed.");
    525         out = psVectorRecycle(out,type,numElements);
    526         out->n = numElements;
    527         memcpy(out->data.V,in->data.V,numElements*PSELEMTYPE_SIZEOF(type));
    528         return out;
    529     }
    530 
    531     if (type == PS_TYPE_C32) {
    532         psF32* outVec;
    533         psC32* inVec = in->data.C32;
    534 
    535         out = psVectorRecycle(out,PS_TYPE_F32,numElements);
    536         out->n = numElements;
    537         outVec = out->data.F32;
    538 
    539         for (unsigned int i=0;i<numElements;i++) {
    540             outVec[i] = crealf(inVec[i]);
    541         }
    542     } else {
    543         psError(__func__,"Can not extract real component from given vector type (%d).",
    544                 type);
    545         psFree(out);
    546         return NULL;
    547     }
    548 
    549     return out;
    550 }
    551 
    552 psVector *psVectorImaginary(psVector *out, const psVector* in)
     528    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     529            // Warn user, as this is probably not expected
     530            psLogMsg( __func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
     531                      "Just a vector copy was performed." );
     532            out = psVectorRecycle( out, type, numElements );
     533            out->n = numElements;
     534            memcpy( out->data.V, in->data.V, numElements * PSELEMTYPE_SIZEOF( type ) );
     535            return out;
     536        }
     537       
     538    if ( type == PS_TYPE_C32 ) {
     539            psF32 * outVec;
     540            psC32* inVec = in->data.C32;
     541           
     542            out = psVectorRecycle( out, PS_TYPE_F32, numElements );
     543            out->n = numElements;
     544            outVec = out->data.F32;
     545           
     546            for ( unsigned int i = 0;i < numElements;i++ ) {
     547                    outVec[ i ] = crealf( inVec[ i ] );
     548                }
     549        } else {
     550            psError( __func__, "Can not extract real component from given vector type (%d).",
     551                     type );
     552            psFree( out );
     553            return NULL;
     554        }
     555       
     556    return out;
     557}
     558
     559psVector *psVectorImaginary( psVector *out, const psVector* in )
    553560{
    554561    psElemType type;
    555562    unsigned int numElements;
    556 
    557 
    558     if (in == NULL) {
    559         psFree(out);
    560         return NULL;
    561     }
    562 
     563   
     564   
     565    if ( in == NULL ) {
     566            psFree( out );
     567            return NULL;
     568        }
     569       
    563570    type = in->type.type;
    564571    numElements = in->n;
    565 
     572   
    566573    /* if not a complex number, this is logically just zeroed image of same size */
    567     if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    568         // Warn user, as this is probably not expected
    569         psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. "
    570                  "A zeroed vector was returned.");
    571         out = psVectorRecycle(out,type,numElements);
    572         out->n = numElements;
    573         memset(out->data.V,0,PSELEMTYPE_SIZEOF(type)*numElements);
    574         return out;
    575     }
    576 
    577     if (type == PS_TYPE_C32) {
    578         psF32* outVec;
    579         psC32* inVec = in->data.C32;
    580 
    581         out = psVectorRecycle(out,PS_TYPE_F32,numElements);
    582         out->n = numElements;
    583         outVec = out->data.F32;
    584 
    585         for (unsigned int i=0;i<numElements;i++) {
    586             outVec[i] = cimagf(inVec[i]);
    587         }
    588     } else {
    589         psError(__func__,"Can not extract imaginary component from given vector type (%d).",
    590                 type);
    591         psFree(out);
    592         return NULL;
    593     }
    594 
    595     return out;
    596 }
    597 
    598 psVector *psVectorComplex(psVector* out, psVector *real, const psVector *imag)
     574    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     575            // Warn user, as this is probably not expected
     576            psLogMsg( __func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
     577                      "A zeroed vector was returned." );
     578            out = psVectorRecycle( out, type, numElements );
     579            out->n = numElements;
     580            memset( out->data.V, 0, PSELEMTYPE_SIZEOF( type ) * numElements );
     581            return out;
     582        }
     583       
     584    if ( type == PS_TYPE_C32 ) {
     585            psF32 * outVec;
     586            psC32* inVec = in->data.C32;
     587           
     588            out = psVectorRecycle( out, PS_TYPE_F32, numElements );
     589            out->n = numElements;
     590            outVec = out->data.F32;
     591           
     592            for ( unsigned int i = 0;i < numElements;i++ ) {
     593                    outVec[ i ] = cimagf( inVec[ i ] );
     594                }
     595        } else {
     596            psError( __func__, "Can not extract imaginary component from given vector type (%d).",
     597                     type );
     598            psFree( out );
     599            return NULL;
     600        }
     601       
     602    return out;
     603}
     604
     605psVector *psVectorComplex( psVector* out, psVector *real, const psVector *imag )
    599606{
    600607    psElemType type;
    601608    unsigned int numElements;
    602 
    603 
    604     if (real == NULL || imag == NULL) {
    605         psFree(out);
    606         return NULL;
    607     }
    608 
     609   
     610   
     611    if ( real == NULL || imag == NULL ) {
     612            psFree( out );
     613            return NULL;
     614        }
     615       
    609616    type = real->type.type;
    610     if (real->n >= imag->n) {
    611         numElements = real->n;
    612     } else {
    613         numElements = imag->n;
    614     }
    615 
    616     if (imag->type.type != type) {
    617         psError(__func__,"The inputs to psVectorComplex must be the same type.");
    618         psFree(out);
    619         return NULL;
    620     }
    621 
    622     if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
    623         psError(__func__,"The inputs to psVectorComplex can not be complex.");
    624         psFree(out);
    625         return NULL;
    626     }
    627 
    628     if (type == PS_TYPE_F32) {
    629         psC32* outVec;
    630         psF32* realVec = real->data.F32;
    631         psF32* imagVec = imag->data.F32;
    632 
    633         out = psVectorRecycle(out,PS_TYPE_C32,numElements);
    634         out->n = numElements;
    635         outVec = out->data.C32;
    636 
    637         for (unsigned int i=0;i<numElements;i++) {
    638             outVec[i] = realVec[i] + I*imagVec[i];
    639         }
    640     } else {
    641         psError(__func__,"Can not merge real and imaginary portions for given vector type (%d).",
    642                 type);
    643         psFree(out);
    644         return NULL;
    645     }
    646 
    647     return out;
    648 }
    649 
    650 psVector *psVectorConjugate(psVector *out, const psVector *in)
     617    if ( real->n >= imag->n ) {
     618            numElements = real->n;
     619        } else {
     620            numElements = imag->n;
     621        }
     622       
     623    if ( imag->type.type != type ) {
     624            psError( __func__, "The inputs to psVectorComplex must be the same type." );
     625            psFree( out );
     626            return NULL;
     627        }
     628       
     629    if ( PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     630            psError( __func__, "The inputs to psVectorComplex can not be complex." );
     631            psFree( out );
     632            return NULL;
     633        }
     634       
     635    if ( type == PS_TYPE_F32 ) {
     636            psC32 * outVec;
     637            psF32* realVec = real->data.F32;
     638            psF32* imagVec = imag->data.F32;
     639           
     640            out = psVectorRecycle( out, PS_TYPE_C32, numElements );
     641            out->n = numElements;
     642            outVec = out->data.C32;
     643           
     644            for ( unsigned int i = 0;i < numElements;i++ ) {
     645                    outVec[ i ] = realVec[ i ] + I * imagVec[ i ];
     646                }
     647        } else {
     648            psError( __func__, "Can not merge real and imaginary portions for given vector type (%d).",
     649                     type );
     650            psFree( out );
     651            return NULL;
     652        }
     653       
     654    return out;
     655}
     656
     657psVector *psVectorConjugate( psVector *out, const psVector *in )
    651658{
    652659    psElemType type;
    653660    unsigned int numElements;
    654 
    655 
    656     if (in == NULL) {
    657         psFree(out);
    658         return NULL;
    659     }
    660 
     661   
     662   
     663    if ( in == NULL ) {
     664            psFree( out );
     665            return NULL;
     666        }
     667       
    661668    type = in->type.type;
    662669    numElements = in->n;
    663 
     670   
    664671    /* if not a complex number, this is logically just a image copy */
    665     if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    666         // Warn user, as this is probably not expected
    667         psLogMsg(__func__,PS_LOG_WARN,"Complex Conjugate of a non-Complex type called for. "
    668                  "Vector copy was performed instead.");
    669 
    670         out = psVectorRecycle(out,type,numElements);
    671         out->n = numElements;
    672         memcpy(out->data.V,in->data.V,PSELEMTYPE_SIZEOF(type)*numElements);
    673         return out;
    674     }
    675 
    676     if (type == PS_TYPE_C32) {
    677         psC32* outVec;
    678         psC32* inVec = in->data.C32;
    679 
    680         out = psVectorRecycle(out,PS_TYPE_C32, numElements);
    681         out->n = numElements;
    682         outVec = out->data.C32;
    683 
    684         for (unsigned int i=0;i<numElements;i++) {
    685             outVec[i] = crealf(inVec[i]) - I*cimagf(inVec[i]);
    686         }
    687     } else {
    688         psError(__func__,"Can not compute complex conjugate for given vector type (%d).",
    689                 type);
    690         psFree(out);
    691         return NULL;
    692     }
    693 
    694     return out;
    695 }
    696 
    697 psVector *psVectorPowerSpectrum(psVector* out, const psVector* in)
     672    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     673            // Warn user, as this is probably not expected
     674            psLogMsg( __func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
     675                      "Vector copy was performed instead." );
     676                     
     677            out = psVectorRecycle( out, type, numElements );
     678            out->n = numElements;
     679            memcpy( out->data.V, in->data.V, PSELEMTYPE_SIZEOF( type ) * numElements );
     680            return out;
     681        }
     682       
     683    if ( type == PS_TYPE_C32 ) {
     684            psC32 * outVec;
     685            psC32* inVec = in->data.C32;
     686           
     687            out = psVectorRecycle( out, PS_TYPE_C32, numElements );
     688            out->n = numElements;
     689            outVec = out->data.C32;
     690           
     691            for ( unsigned int i = 0;i < numElements;i++ ) {
     692                    outVec[ i ] = crealf( inVec[ i ] ) - I * cimagf( inVec[ i ] );
     693                }
     694        } else {
     695            psError( __func__, "Can not compute complex conjugate for given vector type (%d).",
     696                     type );
     697            psFree( out );
     698            return NULL;
     699        }
     700       
     701    return out;
     702}
     703
     704psVector *psVectorPowerSpectrum( psVector* out, const psVector* in )
    698705{
    699706    psElemType type;
    … …  
    702709    unsigned int inHalfNumElements;
    703710    unsigned int inNumElementsSquared;
    704 
    705     if (in == NULL) {
    706         psFree(out);
    707         return NULL;
    708     }
    709 
     711   
     712    if ( in == NULL ) {
     713            psFree( out );
     714            return NULL;
     715        }
     716       
    710717    type = in->type.type;
    711718    inNumElements = in->n;
    712     inNumElementsSquared = inNumElements*inNumElements;
    713     inHalfNumElements = inNumElements/2;
    714     outNumElements = inHalfNumElements+1;
    715 
     719    inNumElementsSquared = inNumElements * inNumElements;
     720    inHalfNumElements = inNumElements / 2;
     721    outNumElements = inHalfNumElements + 1;
     722   
    716723    /* if not a complex number, this is not implemented */
    717     if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    718         psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
    719         psFree(out);
    720         return NULL;
    721     }
    722 
    723     if (type == PS_TYPE_C32) {
    724         psF32* outVec;
    725         psC32* inVec = in->data.C32;
    726         psF32 inAbs1;
    727         psF32 inAbs2;
    728 
    729         out = psVectorRecycle(out,PS_TYPE_F32,outNumElements);
    730         out->n = outNumElements;
    731         outVec = out->data.F32;
    732 
    733         // from ADD: P_0 = |C_0|^2/N^2
    734         inAbs1 = cabsf(inVec[0]);
    735         outVec[0] = inAbs1*inAbs1/inNumElementsSquared;
    736 
    737         // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
    738         for (unsigned int i=1;i<inHalfNumElements;i++) {
    739             inAbs1 = cabsf(inVec[i]);
    740             inAbs2 = cabsf(inVec[inNumElements-i]);
    741             outVec[i] = (inAbs1*inAbs1+inAbs2*inAbs2)/inNumElementsSquared;
    742         }
    743 
    744         // from ADD: P_N/2 = |C_N/2|^2/N^2
    745         inAbs1 = cabsf(inVec[inHalfNumElements]);
    746         outVec[inHalfNumElements] = inAbs1*inAbs1/inNumElementsSquared;
    747     } else {
    748         psError(__func__,"Can not power spectrum for given vector type (%d).",
    749                 type);
    750         psFree(out);
    751         return NULL;
    752     }
    753 
    754     return out;
    755 
    756 }
     724    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     725            psError( __func__, "Power Spectrum for non-complex inputs is not implemented." );
     726            psFree( out );
     727            return NULL;
     728        }
     729       
     730    if ( type == PS_TYPE_C32 ) {
     731            psF32 * outVec;
     732            psC32* inVec = in->data.C32;
     733            psF32 inAbs1;
     734            psF32 inAbs2;
     735           
     736            out = psVectorRecycle( out, PS_TYPE_F32, outNumElements );
     737            out->n = outNumElements;
     738            outVec = out->data.F32;
     739           
     740            // from ADD: P_0 = |C_0|^2/N^2
     741            inAbs1 = cabsf( inVec[ 0 ] );
     742            outVec[ 0 ] = inAbs1 * inAbs1 / inNumElementsSquared;
     743           
     744            // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
     745            for ( unsigned int i = 1;i < inHalfNumElements;i++ ) {
     746                    inAbs1 = cabsf( inVec[ i ] );
     747                    inAbs2 = cabsf( inVec[ inNumElements - i ] );
     748                    outVec[ i ] = ( inAbs1 * inAbs1 + inAbs2 * inAbs2 ) / inNumElementsSquared;
     749                }
     750               
     751            // from ADD: P_N/2 = |C_N/2|^2/N^2
     752            inAbs1 = cabsf( inVec[ inHalfNumElements ] );
     753            outVec[ inHalfNumElements ] = inAbs1 * inAbs1 / inNumElementsSquared;
     754        } else {
     755            psError( __func__, "Can not power spectrum for given vector type (%d).",
     756                     type );
     757            psFree( out );
     758            return NULL;
     759        }
     760       
     761    return out;
     762   
     763}
Note: See TracChangeset for help on using the changeset viewer.