IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 6, 2004, 2:06:06 PM (22 years ago)
Author:
desonia
Message:

another attempt to get astyle to get it right.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/fft/psVectorFFT.c

    r1406 r1407  
     1
    12/** @file  psFFT.c
    23*
    … …  
    56*  @author Robert DeSonia, MHPCC
    67*
    7 *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2004-08-06 22:34:05 $
     8*  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
     9*  @date $Date: 2004-08-07 00:06:06 $
    910*
    1011*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2728static bool p_fftwWisdomImported = false;
    2829
    29 psImage* psImageFFT( psImage* out, const psImage* in, psFftDirection direction )
     30psImage *psImageFFT(psImage * out, const psImage * in, psFftDirection direction)
    3031{
    3132    unsigned int numCols;
    … …  
    3536
    3637    /* 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 
    51     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 
    59     if ( type != PS_TYPE_F32 && direction == PS_FFT_FORWARD ) {
    60         psError( __func__, "Input image must be real image for forward FFT (type=%d).",
    61                  type );
    62         psFree( out );
     38    if (in == NULL) {
     39        psFree(out);
     40        return NULL;
     41    }
     42
     43    type = in->type.type;
     44
     45    if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) {
     46        psError(__func__, "Input image must be a 32-bit float or complex image (type=%d)", type);
     47        psFree(out);
     48        return NULL;
     49    }
     50
     51    if (type != PS_TYPE_C32 && direction == PS_FFT_REVERSE) {
     52        psError(__func__, "Input image must be complex image for reverse FFT (type=%d).", type);
     53        psFree(out);
     54        return NULL;
     55
     56    }
     57
     58    if (type != PS_TYPE_F32 && direction == PS_FFT_FORWARD) {
     59        psError(__func__, "Input image must be real image for forward FFT (type=%d).", type);
     60        psFree(out);
    6361        return NULL;
    6462    }
    6563
    6664    /* make sure the system-level wisdom information is imported. */
    67     if ( ! p_fftwWisdomImported ) {
     65    if (!p_fftwWisdomImported) {
    6866        fftwf_import_system_wisdom();
    6967        p_fftwWisdomImported = true;
    … …  
    7371    numCols = in->numCols;
    7472
    75     out = psImageCopy( out, in, PS_TYPE_C32 );
    76 
    77     plan = fftwf_plan_dft_2d( numCols, numRows,
    78                               ( fftwf_complex* ) out->data.C32[ 0 ],
    79                               ( fftwf_complex* ) out->data.C32[ 0 ],
    80                               direction,
    81                               P_FFTW_PLAN_RIGOR );
    82 
    83     /* check if a plan exists now*/
    84     if ( plan == NULL ) {
    85         psError( __func__, "Failed to create FFTW plan." );
    86         psFree( out );
     73    out = psImageCopy(out, in, PS_TYPE_C32);
     74
     75    plan = fftwf_plan_dft_2d(numCols, numRows,
     76                             (fftwf_complex *) out->data.C32[0],
     77                             (fftwf_complex *) out->data.C32[0], direction, P_FFTW_PLAN_RIGOR);
     78
     79    /* check if a plan exists now */
     80    if (plan == NULL) {
     81        psError(__func__, "Failed to create FFTW plan.");
     82        psFree(out);
    8783        return NULL;
    8884    }
    8985
    9086    /* finally, call FFTW with the plan made above */
    91     fftwf_execute( plan );
    92 
    93     fftwf_destroy_plan( plan );
    94 
    95     return out;
    96 
    97 }
    98 
    99 
    100 psImage *psImageReal( psImage *out, const psImage* in )
     87    fftwf_execute(plan);
     88
     89    fftwf_destroy_plan(plan);
     90
     91    return out;
     92
     93}
     94
     95psImage *psImageReal(psImage * out, const psImage * in)
    10196{
    10297    psElemType type;
    … …  
    10499    unsigned int numRows;
    105100
    106 
    107     if ( in == NULL ) {
    108         psFree( out );
     101    if (in == NULL) {
     102        psFree(out);
    109103        return NULL;
    110104    }
    … …  
    115109
    116110    /* if not a complex number, this is logically just a copy */
    117     if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     111    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    118112        // Warn user, as this is probably not expected
    119         psLogMsg( __func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
    120                   "Just an image copy was performed." );
    121         return psImageCopy( out, in, type );
    122     }
    123 
    124     if ( type == PS_TYPE_C32 ) {
    125         psF32 * outRow;
    126         psC32* inRow;
    127 
    128         out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
    129         for ( unsigned int row = 0;row < numRows;row++ ) {
    130             outRow = out->data.F32[ row ];
    131             inRow = in->data.C32[ row ];
    132 
    133             for ( unsigned int col = 0;col < numCols;col++ ) {
    134                 outRow[ col ] = crealf( inRow[ col ] );
    135             }
    136         }
    137     } else
    138         if ( type == PS_TYPE_C64 ) {
    139             psF64 * outRow;
    140             psC64* inRow;
    141 
    142             out = psImageRecycle( out, numCols, numRows, PS_TYPE_F64 );
    143             for ( unsigned int row = 0;row < numRows;row++ ) {
    144                 outRow = out->data.F64[ row ];
    145                 inRow = in->data.C64[ row ];
    146 
    147                 for ( unsigned int col = 0;col < numCols;col++ ) {
    148                     outRow[ col ] = creal( inRow[ col ] );
    149                 }
    150             }
    151         } else {
    152             psError( __func__, "Can not extract real component from given image type (%d).",
    153                      type );
    154             psFree( out );
    155             return NULL;
    156         }
    157 
    158     return out;
    159 }
    160 
    161 psImage *psImageImaginary( psImage *out, const psImage* in )
     113        psLogMsg(__func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
     114                 "Just an image copy was performed.");
     115        return psImageCopy(out, in, type);
     116    }
     117
     118    if (type == PS_TYPE_C32) {
     119        psF32 *outRow;
     120        psC32 *inRow;
     121
     122        out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
     123        for (unsigned int row = 0; row < numRows; row++) {
     124            outRow = out->data.F32[row];
     125            inRow = in->data.C32[row];
     126
     127            for (unsigned int col = 0; col < numCols; col++) {
     128                outRow[col] = crealf(inRow[col]);
     129            }
     130        }
     131    } else if (type == PS_TYPE_C64) {
     132        psF64 *outRow;
     133        psC64 *inRow;
     134
     135        out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
     136        for (unsigned int row = 0; row < numRows; row++) {
     137            outRow = out->data.F64[row];
     138            inRow = in->data.C64[row];
     139
     140            for (unsigned int col = 0; col < numCols; col++) {
     141                outRow[col] = creal(inRow[col]);
     142            }
     143        }
     144    } else {
     145        psError(__func__, "Can not extract real component from given image type (%d).", type);
     146        psFree(out);
     147        return NULL;
     148    }
     149
     150    return out;
     151}
     152
     153psImage *psImageImaginary(psImage * out, const psImage * in)
    162154{
    163155    psElemType type;
    … …  
    165157    unsigned int numRows;
    166158
    167 
    168     if ( in == NULL ) {
    169         psFree( out );
     159    if (in == NULL) {
     160        psFree(out);
    170161        return NULL;
    171162    }
    … …  
    176167
    177168    /* if not a complex number, this is logically just zeroed image of same size */
    178     if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     169    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    179170        // Warn user, as this is probably not expected
    180         psLogMsg( __func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
    181                   "A zero image was returned." );
    182         out = psImageRecycle( out, numCols, numRows, type );
    183         memset( out->data.V[ 0 ], 0, PSELEMTYPE_SIZEOF( type ) * numCols * numRows );
     171        psLogMsg(__func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
     172                 "A zero image was returned.");
     173        out = psImageRecycle(out, numCols, numRows, type);
     174        memset(out->data.V[0], 0, PSELEMTYPE_SIZEOF(type) * numCols * numRows);
    184175        return out;
    185176    }
    186177
    187     if ( type == PS_TYPE_C32 ) {
    188         psF32 * outRow;
    189         psC32* inRow;
    190 
    191         out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
    192         for ( unsigned int row = 0;row < numRows;row++ ) {
    193             outRow = out->data.F32[ row ];
    194             inRow = in->data.C32[ row ];
    195 
    196             for ( unsigned int col = 0;col < numCols;col++ ) {
    197                 outRow[ col ] = cimagf( inRow[ col ] );
    198             }
    199         }
    200     } else
    201         if ( type == PS_TYPE_C64 ) {
    202             psF64 * outRow;
    203             psC64* inRow;
    204 
    205             out = psImageRecycle( out, numCols, numRows, PS_TYPE_F64 );
    206             for ( unsigned int row = 0;row < numRows;row++ ) {
    207                 outRow = out->data.F64[ row ];
    208                 inRow = in->data.C64[ row ];
    209 
    210                 for ( unsigned int col = 0;col < numCols;col++ ) {
    211                     outRow[ col ] = cimag( inRow[ col ] );
    212                 }
    213             }
    214         } else {
    215             psError( __func__, "Can not extract imaginary component from given image type (%d).",
    216                      type );
    217             psFree( out );
    218             return NULL;
    219         }
    220 
    221     return out;
    222 }
    223 
    224 psImage *psImageComplex( psImage* out, psImage *real, const psImage *imag )
     178    if (type == PS_TYPE_C32) {
     179        psF32 *outRow;
     180        psC32 *inRow;
     181
     182        out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
     183        for (unsigned int row = 0; row < numRows; row++) {
     184            outRow = out->data.F32[row];
     185            inRow = in->data.C32[row];
     186
     187            for (unsigned int col = 0; col < numCols; col++) {
     188                outRow[col] = cimagf(inRow[col]);
     189            }
     190        }
     191    } else if (type == PS_TYPE_C64) {
     192        psF64 *outRow;
     193        psC64 *inRow;
     194
     195        out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
     196        for (unsigned int row = 0; row < numRows; row++) {
     197            outRow = out->data.F64[row];
     198            inRow = in->data.C64[row];
     199
     200            for (unsigned int col = 0; col < numCols; col++) {
     201                outRow[col] = cimag(inRow[col]);
     202            }
     203        }
     204    } else {
     205        psError(__func__, "Can not extract imaginary component from given image type (%d).", type);
     206        psFree(out);
     207        return NULL;
     208    }
     209
     210    return out;
     211}
     212
     213psImage *psImageComplex(psImage * out, psImage * real, const psImage * imag)
    225214{
    226215    psElemType type;
    … …  
    228217    unsigned int numRows;
    229218
    230 
    231     if ( real == NULL || imag == NULL ) {
    232         psFree( out );
     219    if (real == NULL || imag == NULL) {
     220        psFree(out);
    233221        return NULL;
    234222    }
    … …  
    238226    numRows = real->numRows;
    239227
    240     if ( imag->type.type != type ) {
    241         psError( __func__, "The inputs to psImageComplex must be the same type." );
    242         psFree( out );
    243         return NULL;
    244     }
    245 
    246     if ( imag->numCols != numCols ||
    247             imag->numRows != numRows ) {
    248         psError( __func__, "The inputs to psImageComplex must be the same dimensions." );
    249         psFree( out );
    250         return NULL;
    251     }
    252 
    253     if ( PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    254         psError( __func__, "The inputs to psImageComplex can not be complex." );
    255         psFree( out );
    256         return NULL;
    257     }
    258 
    259     if ( type != PS_TYPE_F32 && type != PS_TYPE_F64 ) {
    260         psError( __func__, "The input type to psImageComplex must be a floating point." );
    261         psFree( out );
    262         return NULL;
    263     }
    264 
    265     if ( type == PS_TYPE_F32 ) {
    266         psC32 * outRow;
    267         psF32* realRow;
    268         psF32* imagRow;
    269 
    270         out = psImageRecycle( out, numCols, numRows, PS_TYPE_C32 );
    271 
    272         for ( unsigned int row = 0;row < numRows;row++ ) {
    273             outRow = out->data.C32[ row ];
    274             realRow = real->data.F32[ row ];
    275             imagRow = imag->data.F32[ row ];
    276 
    277             for ( unsigned int col = 0;col < numCols;col++ ) {
    278                 outRow[ col ] = realRow[ col ] + I * imagRow[ col ];
    279             }
    280         }
    281     } else
    282         if ( type == PS_TYPE_F64 ) {
    283             psC64 * outRow;
    284             psF64* realRow;
    285             psF64* imagRow;
    286 
    287             out = psImageRecycle( out, numCols, numRows, PS_TYPE_C64 );
    288             for ( unsigned int row = 0;row < numRows;row++ ) {
    289                 outRow = out->data.C64[ row ];
    290                 realRow = real->data.F64[ row ];
    291                 imagRow = imag->data.F64[ row ];
    292 
    293                 for ( unsigned int col = 0;col < numCols;col++ ) {
    294                     outRow[ col ] = realRow[ col ] + I * imagRow[ col ];
    295                 }
    296             }
    297         } else {
    298             psError( __func__, "Can not merge real and imaginary portions for given image type (%d).",
    299                      type );
    300             psFree( out );
    301             return NULL;
    302         }
    303 
    304     return out;
    305 }
    306 
    307 psImage *psImageConjugate( psImage *out, const psImage *in )
     228    if (imag->type.type != type) {
     229        psError(__func__, "The inputs to psImageComplex must be the same type.");
     230        psFree(out);
     231        return NULL;
     232    }
     233
     234    if (imag->numCols != numCols || imag->numRows != numRows) {
     235        psError(__func__, "The inputs to psImageComplex must be the same dimensions.");
     236        psFree(out);
     237        return NULL;
     238    }
     239
     240    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
     241        psError(__func__, "The inputs to psImageComplex can not be complex.");
     242        psFree(out);
     243        return NULL;
     244    }
     245
     246    if (type != PS_TYPE_F32 && type != PS_TYPE_F64) {
     247        psError(__func__, "The input type to psImageComplex must be a floating point.");
     248        psFree(out);
     249        return NULL;
     250    }
     251
     252    if (type == PS_TYPE_F32) {
     253        psC32 *outRow;
     254        psF32 *realRow;
     255        psF32 *imagRow;
     256
     257        out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32);
     258
     259        for (unsigned int row = 0; row < numRows; row++) {
     260            outRow = out->data.C32[row];
     261            realRow = real->data.F32[row];
     262            imagRow = imag->data.F32[row];
     263
     264            for (unsigned int col = 0; col < numCols; col++) {
     265                outRow[col] = realRow[col] + I * imagRow[col];
     266            }
     267        }
     268    } else if (type == PS_TYPE_F64) {
     269        psC64 *outRow;
     270        psF64 *realRow;
     271        psF64 *imagRow;
     272
     273        out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64);
     274        for (unsigned int row = 0; row < numRows; row++) {
     275            outRow = out->data.C64[row];
     276            realRow = real->data.F64[row];
     277            imagRow = imag->data.F64[row];
     278
     279            for (unsigned int col = 0; col < numCols; col++) {
     280                outRow[col] = realRow[col] + I * imagRow[col];
     281            }
     282        }
     283    } else {
     284        psError(__func__, "Can not merge real and imaginary portions for given image type (%d).", type);
     285        psFree(out);
     286        return NULL;
     287    }
     288
     289    return out;
     290}
     291
     292psImage *psImageConjugate(psImage * out, const psImage * in)
    308293{
    309294    psElemType type;
    … …  
    311296    unsigned int numRows;
    312297
    313 
    314     if ( in == NULL ) {
    315         psFree( out );
     298    if (in == NULL) {
     299        psFree(out);
    316300        return NULL;
    317301    }
    … …  
    322306
    323307    /* if not a complex number, this is logically just a image copy */
    324     if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     308    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    325309        // Warn user, as this is probably not expected
    326         psLogMsg( __func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
    327                   "Image copy was performed instead." );
    328         return psImageCopy( out, in, type );
    329     }
    330 
    331     if ( type == PS_TYPE_C32 ) {
    332         psC32 * outRow;
    333         psC32* inRow;
    334 
    335         out = psImageRecycle( out, numCols, numRows, PS_TYPE_C32 );
    336         for ( unsigned int row = 0;row < numRows;row++ ) {
    337             outRow = out->data.C32[ row ];
    338             inRow = in->data.C32[ row ];
    339 
    340             for ( unsigned int col = 0;col < numCols;col++ ) {
    341                 outRow[ col ] = crealf( inRow[ col ] ) - I * cimagf( inRow[ col ] );
    342             }
    343         }
    344     } else
    345         if ( type == PS_TYPE_C64 ) {
    346             psC64 * outRow;
    347             psC64* inRow;
    348 
    349             out = psImageRecycle( out, numCols, numRows, PS_TYPE_C64 );
    350             for ( unsigned int row = 0;row < numRows;row++ ) {
    351                 outRow = out->data.C64[ row ];
    352                 inRow = in->data.C64[ row ];
    353 
    354                 for ( unsigned int col = 0;col < numCols;col++ ) {
    355                     outRow[ col ] = creal( inRow[ col ] ) - I * cimag( inRow[ col ] );
    356                 }
    357             }
    358         } else {
    359             psError( __func__, "Can not compute complex conjugate for given image type (%d).",
    360                      type );
    361             psFree( out );
    362             return NULL;
    363         }
    364 
    365     return out;
    366 }
    367 
    368 psImage *psImagePowerSpectrum( psImage* out, const psImage* in )
     310        psLogMsg(__func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
     311                 "Image copy was performed instead.");
     312        return psImageCopy(out, in, type);
     313    }
     314
     315    if (type == PS_TYPE_C32) {
     316        psC32 *outRow;
     317        psC32 *inRow;
     318
     319        out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32);
     320        for (unsigned int row = 0; row < numRows; row++) {
     321            outRow = out->data.C32[row];
     322            inRow = in->data.C32[row];
     323
     324            for (unsigned int col = 0; col < numCols; col++) {
     325                outRow[col] = crealf(inRow[col]) - I * cimagf(inRow[col]);
     326            }
     327        }
     328    } else if (type == PS_TYPE_C64) {
     329        psC64 *outRow;
     330        psC64 *inRow;
     331
     332        out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64);
     333        for (unsigned int row = 0; row < numRows; row++) {
     334            outRow = out->data.C64[row];
     335            inRow = in->data.C64[row];
     336
     337            for (unsigned int col = 0; col < numCols; col++) {
     338                outRow[col] = creal(inRow[col]) - I * cimag(inRow[col]);
     339            }
     340        }
     341    } else {
     342        psError(__func__, "Can not compute complex conjugate for given image type (%d).", type);
     343        psFree(out);
     344        return NULL;
     345    }
     346
     347    return out;
     348}
     349
     350psImage *psImagePowerSpectrum(psImage * out, const psImage * in)
    369351{
    370352    psElemType type;
    … …  
    373355    int numElementsSquared;
    374356
    375     if ( in == NULL ) {
    376         psFree( out );
     357    if (in == NULL) {
     358        psFree(out);
    377359        return NULL;
    378360    }
    … …  
    384366
    385367    /* if not a complex number, this is not implemented */
    386     if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    387         psError( __func__, "Power Spectrum for non-complex inputs is not implemented." );
    388         psFree( out );
    389         return NULL;
    390     }
    391 
    392     if ( type == PS_TYPE_C32 ) {
    393         psF32 * outRow;
    394         psC32* inRow;
     368    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
     369        psError(__func__, "Power Spectrum for non-complex inputs is not implemented.");
     370        psFree(out);
     371        return NULL;
     372    }
     373
     374    if (type == PS_TYPE_C32) {
     375        psF32 *outRow;
     376        psC32 *inRow;
    395377        psF32 real;
    396378        psF32 imag;
    397379
    398 
    399         out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
    400         for ( unsigned int row = 0;row < numRows;row++ ) {
    401             outRow = out->data.F32[ row ];
    402             inRow = in->data.C32[ row ];
    403 
    404             for ( unsigned int col = 0;col < numCols;col++ ) {
    405                 real = crealf( inRow[ col ] );
    406                 imag = cimagf( inRow[ col ] );
    407                 outRow[ col ] = ( real * real + imag * imag ) / numElementsSquared;
    408             }
    409         }
    410     } else
    411         if ( type == PS_TYPE_C64 ) {
    412             psF64 * outRow;
    413             psC64* inRow;
    414             psF64 real;
    415             psF64 imag;
    416 
    417 
    418             out = psImageRecycle( out, numCols, numRows, PS_TYPE_F64 );
    419             for ( unsigned int row = 0;row < numRows;row++ ) {
    420                 outRow = out->data.F64[ row ];
    421                 inRow = in->data.C64[ row ];
    422 
    423                 for ( unsigned int col = 0;col < numCols;col++ ) {
    424                     real = crealf( inRow[ col ] );
    425                     imag = cimagf( inRow[ col ] );
    426                     outRow[ col ] = real * real + imag * imag / numElementsSquared;
    427                 }
    428             }
    429         } else {
    430             psError( __func__, "Can not power spectrum for given image type (%d).",
    431                      type );
    432             psFree( out );
    433             return NULL;
    434         }
     380        out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
     381        for (unsigned int row = 0; row < numRows; row++) {
     382            outRow = out->data.F32[row];
     383            inRow = in->data.C32[row];
     384
     385            for (unsigned int col = 0; col < numCols; col++) {
     386                real = crealf(inRow[col]);
     387                imag = cimagf(inRow[col]);
     388                outRow[col] = (real * real + imag * imag) / numElementsSquared;
     389            }
     390        }
     391    } else if (type == PS_TYPE_C64) {
     392        psF64 *outRow;
     393        psC64 *inRow;
     394        psF64 real;
     395        psF64 imag;
     396
     397        out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
     398        for (unsigned int row = 0; row < numRows; row++) {
     399            outRow = out->data.F64[row];
     400            inRow = in->data.C64[row];
     401
     402            for (unsigned int col = 0; col < numCols; col++) {
     403                real = crealf(inRow[col]);
     404                imag = cimagf(inRow[col]);
     405                outRow[col] = real * real + imag * imag / numElementsSquared;
     406            }
     407        }
     408    } else {
     409        psError(__func__, "Can not power spectrum for given image type (%d).", type);
     410        psFree(out);
     411        return NULL;
     412    }
    435413
    436414    return out;
    … …  
    440418/************************************** Vector Functions ***************************************/
    441419
    442 psVector* psVectorFFT( psVector* out, const psVector* in, psFftDirection direction )
     420psVector *psVectorFFT(psVector * out, const psVector * in, psFftDirection direction)
    443421{
    444422    unsigned int numElements;
    … …  
    447425
    448426    /* got good image data? */
    449     if ( in == NULL ) {
    450         psFree( out );
    451         return NULL;
    452     }
    453 
    454     type = in->type.type;
    455 
    456     if ( ( type != PS_TYPE_F32 ) && ( type != PS_TYPE_C32 ) ) {
    457         psError( __func__, "Input image must be a 32-bit float or complex image (type=%d)",
    458                  type );
    459         psFree( out );
    460         return NULL;
    461     }
    462 
    463     if ( ( type != PS_TYPE_C32 ) && ( direction == PS_FFT_REVERSE ) ) {
    464         psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
    465                  type );
    466         psFree( out );
    467         return NULL;
    468 
    469     }
    470 
    471     if ( ( type != PS_TYPE_F32 ) && ( direction == PS_FFT_FORWARD ) ) {
    472         psError( __func__, "Input image must be real image for forward FFT (type=%d).",
    473                  type );
    474         psFree( out );
     427    if (in == NULL) {
     428        psFree(out);
     429        return NULL;
     430    }
     431
     432    type = in->type.type;
     433
     434    if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) {
     435        psError(__func__, "Input image must be a 32-bit float or complex image (type=%d)", type);
     436        psFree(out);
     437        return NULL;
     438    }
     439
     440    if ((type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE)) {
     441        psError(__func__, "Input image must be complex image for reverse FFT (type=%d).", type);
     442        psFree(out);
     443        return NULL;
     444
     445    }
     446
     447    if ((type != PS_TYPE_F32) && (direction == PS_FFT_FORWARD)) {
     448        psError(__func__, "Input image must be real image for forward FFT (type=%d).", type);
     449        psFree(out);
    475450        return NULL;
    476451    }
    477452
    478453    /* make sure the system-level wisdom information is imported. */
    479     if ( ! p_fftwWisdomImported ) {
     454    if (!p_fftwWisdomImported) {
    480455        fftwf_import_system_wisdom();
    481456        p_fftwWisdomImported = true;
    … …  
    484459    numElements = in->n;
    485460
    486     out = psVectorRecycle( out, numElements, PS_TYPE_C32 );
     461    out = psVectorRecycle(out, numElements, PS_TYPE_C32);
    487462    out->n = numElements;
    488463
    489     if ( type == PS_TYPE_F32 ) {
     464    if (type == PS_TYPE_F32) {
    490465        // need to convert to complex
    491         psC32 * outVec = out->data.C32;
    492         psF32* inVec = in->data.F32;
    493         for ( unsigned int i = 0;i < numElements;i++ ) {
    494             outVec[ i ] = inVec[ i ];
    495         }
    496     } else {
    497         psC32* outVec = out->data.C32;
    498         psC32* inVec = in->data.C32;
    499         for ( unsigned int i = 0;i < numElements;i++ ) {
    500             outVec[ i ] = inVec[ i ];
    501         }
    502     }
    503 
    504     plan = fftwf_plan_dft_1d( numElements,
    505                               ( fftwf_complex* ) out->data.C32,
    506                               ( fftwf_complex* ) out->data.C32,
    507                               direction,
    508                               P_FFTW_PLAN_RIGOR );
    509 
    510     /* check if a plan exists now*/
    511     if ( plan == NULL ) {
    512         psError( __func__, "Failed to create FFTW plan." );
    513         psFree( out );
     466        psC32 *outVec = out->data.C32;
     467        psF32 *inVec = in->data.F32;
     468
     469        for (unsigned int i = 0; i < numElements; i++) {
     470            outVec[i] = inVec[i];
     471        }
     472    } else {
     473        psC32 *outVec = out->data.C32;
     474        psC32 *inVec = in->data.C32;
     475
     476        for (unsigned int i = 0; i < numElements; i++) {
     477            outVec[i] = inVec[i];
     478        }
     479    }
     480
     481    plan = fftwf_plan_dft_1d(numElements,
     482                             (fftwf_complex *) out->data.C32,
     483                             (fftwf_complex *) out->data.C32, direction, P_FFTW_PLAN_RIGOR);
     484
     485    /* check if a plan exists now */
     486    if (plan == NULL) {
     487        psError(__func__, "Failed to create FFTW plan.");
     488        psFree(out);
    514489        return NULL;
    515490    }
    516491
    517492    /* finally, call FFTW with the plan made above */
    518     fftwf_execute( plan );
    519 
    520     fftwf_destroy_plan( plan );
    521 
    522     return out;
    523 }
    524 
    525 
    526 psVector *psVectorReal( psVector *out, const psVector* in )
     493    fftwf_execute(plan);
     494
     495    fftwf_destroy_plan(plan);
     496
     497    return out;
     498}
     499
     500psVector *psVectorReal(psVector * out, const psVector * in)
    527501{
    528502    psElemType type;
    529503    unsigned int numElements;
    530504
    531     if ( in == NULL ) {
    532         psFree( out );
     505    if (in == NULL) {
     506        psFree(out);
    533507        return NULL;
    534508    }
    … …  
    538512
    539513    /* if not a complex number, this is logically just a copy */
    540     if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     514    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    541515        // Warn user, as this is probably not expected
    542         psLogMsg( __func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
    543                   "Just a vector copy was performed." );
    544         out = psVectorRecycle( out, numElements, type );
     516        psLogMsg(__func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
     517                 "Just a vector copy was performed.");
     518        out = psVectorRecycle(out, numElements, type);
    545519        out->n = numElements;
    546         memcpy( out->data.V, in->data.V, numElements * PSELEMTYPE_SIZEOF( type ) );
     520        memcpy(out->data.V, in->data.V, numElements * PSELEMTYPE_SIZEOF(type));
    547521        return out;
    548522    }
    549523
    550     if ( type == PS_TYPE_C32 ) {
    551         psF32 * outVec;
    552         psC32* inVec = in->data.C32;
    553 
    554         out = psVectorRecycle( out, numElements, PS_TYPE_F32 );
     524    if (type == PS_TYPE_C32) {
     525        psF32 *outVec;
     526        psC32 *inVec = in->data.C32;
     527
     528        out = psVectorRecycle(out, numElements, PS_TYPE_F32);
    555529        out->n = numElements;
    556530        outVec = out->data.F32;
    557531
    558         for ( unsigned int i = 0;i < numElements;i++ ) {
    559             outVec[ i ] = crealf( inVec[ i ] );
    560         }
    561     } else {
    562         psError( __func__, "Can not extract real component from given vector type (%d).",
    563                  type );
    564         psFree( out );
    565         return NULL;
    566     }
    567 
    568     return out;
    569 }
    570 
    571 psVector *psVectorImaginary( psVector *out, const psVector* in )
     532        for (unsigned int i = 0; i < numElements; i++) {
     533            outVec[i] = crealf(inVec[i]);
     534        }
     535    } else {
     536        psError(__func__, "Can not extract real component from given vector type (%d).", type);
     537        psFree(out);
     538        return NULL;
     539    }
     540
     541    return out;
     542}
     543
     544psVector *psVectorImaginary(psVector * out, const psVector * in)
    572545{
    573546    psElemType type;
    574547    unsigned int numElements;
    575548
    576 
    577     if ( in == NULL ) {
    578         psFree( out );
     549    if (in == NULL) {
     550        psFree(out);
    579551        return NULL;
    580552    }
    … …  
    584556
    585557    /* if not a complex number, this is logically just zeroed image of same size */
    586     if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     558    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    587559        // Warn user, as this is probably not expected
    588         psLogMsg( __func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
    589                   "A zeroed vector was returned." );
    590         out = psVectorRecycle( out, numElements, type );
     560        psLogMsg(__func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
     561                 "A zeroed vector was returned.");
     562        out = psVectorRecycle(out, numElements, type);
    591563        out->n = numElements;
    592         memset( out->data.V, 0, PSELEMTYPE_SIZEOF( type ) * numElements );
     564        memset(out->data.V, 0, PSELEMTYPE_SIZEOF(type) * numElements);
    593565        return out;
    594566    }
    595567
    596     if ( type == PS_TYPE_C32 ) {
    597         psF32 * outVec;
    598         psC32* inVec = in->data.C32;
    599 
    600         out = psVectorRecycle( out, numElements, PS_TYPE_F32 );
     568    if (type == PS_TYPE_C32) {
     569        psF32 *outVec;
     570        psC32 *inVec = in->data.C32;
     571
     572        out = psVectorRecycle(out, numElements, PS_TYPE_F32);
    601573        out->n = numElements;
    602574        outVec = out->data.F32;
    603575
    604         for ( unsigned int i = 0;i < numElements;i++ ) {
    605             outVec[ i ] = cimagf( inVec[ i ] );
    606         }
    607     } else {
    608         psError( __func__, "Can not extract imaginary component from given vector type (%d).",
    609                  type );
    610         psFree( out );
    611         return NULL;
    612     }
    613 
    614     return out;
    615 }
    616 
    617 psVector *psVectorComplex( psVector* out, psVector *real, const psVector *imag )
     576        for (unsigned int i = 0; i < numElements; i++) {
     577            outVec[i] = cimagf(inVec[i]);
     578        }
     579    } else {
     580        psError(__func__, "Can not extract imaginary component from given vector type (%d).", type);
     581        psFree(out);
     582        return NULL;
     583    }
     584
     585    return out;
     586}
     587
     588psVector *psVectorComplex(psVector * out, psVector * real, const psVector * imag)
    618589{
    619590    psElemType type;
    620591    unsigned int numElements;
    621592
    622 
    623     if ( real == NULL || imag == NULL ) {
    624         psFree( out );
     593    if (real == NULL || imag == NULL) {
     594        psFree(out);
    625595        return NULL;
    626596    }
    627597
    628598    type = real->type.type;
    629     if ( real->n < imag->n ) {
     599    if (real->n < imag->n) {
    630600        numElements = real->n;
    631601    } else {
    … …  
    633603    }
    634604
    635     if ( imag->type.type != type ) {
    636         psError( __func__, "The inputs to psVectorComplex must be the same type." );
    637         psFree( out );
    638         return NULL;
    639     }
    640 
    641     if ( PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    642         psError( __func__, "The inputs to psVectorComplex can not be complex." );
    643         psFree( out );
    644         return NULL;
    645     }
    646 
    647     if ( type == PS_TYPE_F32 ) {
    648         psC32 * outVec;
    649         psF32* realVec = real->data.F32;
    650         psF32* imagVec = imag->data.F32;
    651 
    652         out = psVectorRecycle( out, numElements, PS_TYPE_C32 );
     605    if (imag->type.type != type) {
     606        psError(__func__, "The inputs to psVectorComplex must be the same type.");
     607        psFree(out);
     608        return NULL;
     609    }
     610
     611    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
     612        psError(__func__, "The inputs to psVectorComplex can not be complex.");
     613        psFree(out);
     614        return NULL;
     615    }
     616
     617    if (type == PS_TYPE_F32) {
     618        psC32 *outVec;
     619        psF32 *realVec = real->data.F32;
     620        psF32 *imagVec = imag->data.F32;
     621
     622        out = psVectorRecycle(out, numElements, PS_TYPE_C32);
    653623        out->n = numElements;
    654624        outVec = out->data.C32;
    655625
    656         for ( unsigned int i = 0;i < numElements;i++ ) {
    657             outVec[ i ] = realVec[ i ] + I * imagVec[ i ];
    658         }
    659     } else {
    660         psError( __func__, "Can not merge real and imaginary portions for given vector type (%d).",
    661                  type );
    662         psFree( out );
    663         return NULL;
    664     }
    665 
    666     return out;
    667 }
    668 
    669 psVector *psVectorConjugate( psVector *out, const psVector *in )
     626        for (unsigned int i = 0; i < numElements; i++) {
     627            outVec[i] = realVec[i] + I * imagVec[i];
     628        }
     629    } else {
     630        psError(__func__, "Can not merge real and imaginary portions for given vector type (%d).", type);
     631        psFree(out);
     632        return NULL;
     633    }
     634
     635    return out;
     636}
     637
     638psVector *psVectorConjugate(psVector * out, const psVector * in)
    670639{
    671640    psElemType type;
    672641    unsigned int numElements;
    673642
    674 
    675     if ( in == NULL ) {
    676         psFree( out );
     643    if (in == NULL) {
     644        psFree(out);
    677645        return NULL;
    678646    }
    … …  
    682650
    683651    /* if not a complex number, this is logically just a image copy */
    684     if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
     652    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    685653        // Warn user, as this is probably not expected
    686         psLogMsg( __func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
    687                   "Vector copy was performed instead." );
    688 
    689         out = psVectorRecycle( out, numElements, type );
     654        psLogMsg(__func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
     655                 "Vector copy was performed instead.");
     656
     657        out = psVectorRecycle(out, numElements, type);
    690658        out->n = numElements;
    691         memcpy( out->data.V, in->data.V, PSELEMTYPE_SIZEOF( type ) * numElements );
     659        memcpy(out->data.V, in->data.V, PSELEMTYPE_SIZEOF(type) * numElements);
    692660        return out;
    693661    }
    694662
    695     if ( type == PS_TYPE_C32 ) {
    696         psC32 * outVec;
    697         psC32* inVec = in->data.C32;
    698 
    699         out = psVectorRecycle( out, numElements, PS_TYPE_C32 );
     663    if (type == PS_TYPE_C32) {
     664        psC32 *outVec;
     665        psC32 *inVec = in->data.C32;
     666
     667        out = psVectorRecycle(out, numElements, PS_TYPE_C32);
    700668        out->n = numElements;
    701669        outVec = out->data.C32;
    702670
    703         for ( unsigned int i = 0;i < numElements;i++ ) {
    704             outVec[ i ] = crealf( inVec[ i ] ) - I * cimagf( inVec[ i ] );
    705         }
    706     } else {
    707         psError( __func__, "Can not compute complex conjugate for given vector type (%d).",
    708                  type );
    709         psFree( out );
    710         return NULL;
    711     }
    712 
    713     return out;
    714 }
    715 
    716 psVector *psVectorPowerSpectrum( psVector* out, const psVector* in )
     671        for (unsigned int i = 0; i < numElements; i++) {
     672            outVec[i] = crealf(inVec[i]) - I * cimagf(inVec[i]);
     673        }
     674    } else {
     675        psError(__func__, "Can not compute complex conjugate for given vector type (%d).", type);
     676        psFree(out);
     677        return NULL;
     678    }
     679
     680    return out;
     681}
     682
     683psVector *psVectorPowerSpectrum(psVector * out, const psVector * in)
    717684{
    718685    psElemType type;
    … …  
    722689    unsigned int inNumElementsSquared;
    723690
    724     if ( in == NULL ) {
    725         psFree( out );
     691    if (in == NULL) {
     692        psFree(out);
    726693        return NULL;
    727694    }
    … …  
    734701
    735702    /* if not a complex number, this is not implemented */
    736     if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    737         psError( __func__, "Power Spectrum for non-complex inputs is not implemented." );
    738         psFree( out );
    739         return NULL;
    740     }
    741 
    742     if ( type == PS_TYPE_C32 ) {
    743         psF32 * outVec;
    744         psC32* inVec = in->data.C32;
     703    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
     704        psError(__func__, "Power Spectrum for non-complex inputs is not implemented.");
     705        psFree(out);
     706        return NULL;
     707    }
     708
     709    if (type == PS_TYPE_C32) {
     710        psF32 *outVec;
     711        psC32 *inVec = in->data.C32;
    745712        psF32 inAbs1;
    746713        psF32 inAbs2;
    747714
    748         out = psVectorRecycle( out, outNumElements, PS_TYPE_F32 );
     715        out = psVectorRecycle(out, outNumElements, PS_TYPE_F32);
    749716        out->n = outNumElements;
    750717        outVec = out->data.F32;
    751718
    752719        // from ADD: P_0 = |C_0|^2/N^2
    753         inAbs1 = cabsf( inVec[ 0 ] );
    754         outVec[ 0 ] = inAbs1 * inAbs1 / inNumElementsSquared;
     720        inAbs1 = cabsf(inVec[0]);
     721        outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared;
    755722
    756723        // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
    757         for ( unsigned int i = 1;i < inHalfNumElements;i++ ) {
    758             inAbs1 = cabsf( inVec[ i ] );
    759             inAbs2 = cabsf( inVec[ inNumElements - i ] );
    760             outVec[ i ] = ( inAbs1 * inAbs1 + inAbs2 * inAbs2 ) / inNumElementsSquared;
     724        for (unsigned int i = 1; i < inHalfNumElements; i++) {
     725            inAbs1 = cabsf(inVec[i]);
     726            inAbs2 = cabsf(inVec[inNumElements - i]);
     727            outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared;
    761728        }
    762729
    763730        // from ADD: P_N/2 = |C_N/2|^2/N^2
    764         inAbs1 = cabsf( inVec[ inHalfNumElements ] );
    765         outVec[ inHalfNumElements ] = inAbs1 * inAbs1 / inNumElementsSquared;
    766     } else {
    767         psError( __func__, "Can not power spectrum for given vector type (%d).",
    768                  type );
    769         psFree( out );
    770         return NULL;
    771     }
    772 
    773     return out;
    774 
    775 }
     731        inAbs1 = cabsf(inVec[inHalfNumElements]);
     732        outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
     733    } else {
     734        psError(__func__, "Can not power spectrum for given vector type (%d).", type);
     735        psFree(out);
     736        return NULL;
     737    }
     738
     739    return out;
     740
     741}
Note: See TracChangeset for help on using the changeset viewer.