IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 25, 2004, 11:10:09 AM (22 years ago)
Author:
desonia
Message:

extracted image FFT functions into a new file.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1624 r1625  
    1 /** @file  psFFT.c
     1/** @file  psVectorFFT.c
    22 *
    3  *  @brief Contains FFT transforms functions
     3 *  @brief Contains FFT transform related functions for psVector
    44 *
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2004-08-25 20:51:57 $
     7 *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2004-08-25 21:10:08 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2626
    2727static bool p_fftwWisdomImported = false;
    28 
    29 psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction)
    30 {
    31     unsigned int numCols;
    32     unsigned int numRows;
    33     psElemType type;
    34     fftwf_plan plan;
    35 
    36     /* 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)", type);
    46         psFree(out);
    47         return NULL;
    48     }
    49 
    50     if (type != PS_TYPE_C32 && direction == PS_FFT_REVERSE) {
    51         psError(__func__, "Input image must be complex image for reverse FFT (type=%d).", type);
    52         psFree(out);
    53         return NULL;
    54 
    55     }
    56 
    57     if (type != PS_TYPE_F32 && direction == PS_FFT_FORWARD) {
    58         psError(__func__, "Input image must be real image for forward FFT (type=%d).", type);
    59         psFree(out);
    60         return NULL;
    61     }
    62 
    63     /* make sure the system-level wisdom information is imported. */
    64     if (!p_fftwWisdomImported) {
    65         fftwf_import_system_wisdom();
    66         p_fftwWisdomImported = true;
    67     }
    68 
    69     numRows = in->numRows;
    70     numCols = in->numCols;
    71 
    72     out = psImageCopy(out, in, PS_TYPE_C32);
    73 
    74     plan = fftwf_plan_dft_2d(numCols, numRows,
    75                              (fftwf_complex *) out->data.C32[0],
    76                              (fftwf_complex *) out->data.C32[0], direction, P_FFTW_PLAN_RIGOR);
    77 
    78     /* check if a plan exists now */
    79     if (plan == NULL) {
    80         psError(__func__, "Failed to create FFTW plan.");
    81         psFree(out);
    82         return NULL;
    83     }
    84 
    85     /* finally, call FFTW with the plan made above */
    86     fftwf_execute(plan);
    87 
    88     fftwf_destroy_plan(plan);
    89 
    90     return out;
    91 
    92 }
    93 
    94 psImage* psImageReal(psImage* out, const psImage* in)
    95 {
    96     psElemType type;
    97     unsigned int numCols;
    98     unsigned int numRows;
    99 
    100     if (in == NULL) {
    101         psFree(out);
    102         return NULL;
    103     }
    104 
    105     type = in->type.type;
    106     numCols = in->numCols;
    107     numRows = in->numRows;
    108 
    109     /* 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).", type);
    145         psFree(out);
    146         return NULL;
    147     }
    148 
    149     return out;
    150 }
    151 
    152 psImage* psImageImaginary(psImage* out, const psImage* in)
    153 {
    154     psElemType type;
    155     unsigned int numCols;
    156     unsigned int numRows;
    157 
    158     if (in == NULL) {
    159         psFree(out);
    160         return NULL;
    161     }
    162 
    163     type = in->type.type;
    164     numCols = in->numCols;
    165     numRows = in->numRows;
    166 
    167     /* if not a complex number, this is logically just zeroed image of same size */
    168     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    169         // Warn user, as this is probably not expected
    170         psLogMsg(__func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
    171                  "A zero image was returned.");
    172         out = psImageRecycle(out, numCols, numRows, type);
    173         memset(out->data.V[0], 0, PSELEMTYPE_SIZEOF(type) * numCols * numRows);
    174         return out;
    175     }
    176 
    177     if (type == PS_TYPE_C32) {
    178         psF32* outRow;
    179         psC32* inRow;
    180 
    181         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
    182         for (unsigned int row = 0; row < numRows; row++) {
    183             outRow = out->data.F32[row];
    184             inRow = in->data.C32[row];
    185 
    186             for (unsigned int col = 0; col < numCols; col++) {
    187                 outRow[col] = cimagf(inRow[col]);
    188             }
    189         }
    190     } else if (type == PS_TYPE_C64) {
    191         psF64* outRow;
    192         psC64* inRow;
    193 
    194         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
    195         for (unsigned int row = 0; row < numRows; row++) {
    196             outRow = out->data.F64[row];
    197             inRow = in->data.C64[row];
    198 
    199             for (unsigned int col = 0; col < numCols; col++) {
    200                 outRow[col] = cimag(inRow[col]);
    201             }
    202         }
    203     } else {
    204         psError(__func__, "Can not extract imaginary component from given image type (%d).", type);
    205         psFree(out);
    206         return NULL;
    207     }
    208 
    209     return out;
    210 }
    211 
    212 psImage* psImageComplex(psImage* out, const psImage* real, const psImage* imag)
    213 {
    214     psElemType type;
    215     unsigned int numCols;
    216     unsigned int numRows;
    217 
    218     if (real == NULL || imag == NULL) {
    219         psFree(out);
    220         return NULL;
    221     }
    222 
    223     type = real->type.type;
    224     numCols = real->numCols;
    225     numRows = real->numRows;
    226 
    227     if (imag->type.type != type) {
    228         psError(__func__, "The inputs to psImageComplex must be the same type.");
    229         psFree(out);
    230         return NULL;
    231     }
    232 
    233     if (imag->numCols != numCols || imag->numRows != numRows) {
    234         psError(__func__, "The inputs to psImageComplex must be the same dimensions.");
    235         psFree(out);
    236         return NULL;
    237     }
    238 
    239     if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
    240         psError(__func__, "The inputs to psImageComplex can not be complex.");
    241         psFree(out);
    242         return NULL;
    243     }
    244 
    245     if (type != PS_TYPE_F32 && type != PS_TYPE_F64) {
    246         psError(__func__, "The input type to psImageComplex must be a floating point.");
    247         psFree(out);
    248         return NULL;
    249     }
    250 
    251     if (type == PS_TYPE_F32) {
    252         psC32* outRow;
    253         psF32* realRow;
    254         psF32* imagRow;
    255 
    256         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32);
    257 
    258         for (unsigned int row = 0; row < numRows; row++) {
    259             outRow = out->data.C32[row];
    260             realRow = real->data.F32[row];
    261             imagRow = imag->data.F32[row];
    262 
    263             for (unsigned int col = 0; col < numCols; col++) {
    264                 outRow[col] = realRow[col] + I * imagRow[col];
    265             }
    266         }
    267     } else if (type == PS_TYPE_F64) {
    268         psC64* outRow;
    269         psF64* realRow;
    270         psF64* imagRow;
    271 
    272         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64);
    273         for (unsigned int row = 0; row < numRows; row++) {
    274             outRow = out->data.C64[row];
    275             realRow = real->data.F64[row];
    276             imagRow = imag->data.F64[row];
    277 
    278             for (unsigned int col = 0; col < numCols; col++) {
    279                 outRow[col] = realRow[col] + I * imagRow[col];
    280             }
    281         }
    282     } else {
    283         psError(__func__, "Can not merge real and imaginary portions for given image type (%d).", type);
    284         psFree(out);
    285         return NULL;
    286     }
    287 
    288     return out;
    289 }
    290 
    291 psImage* psImageConjugate(psImage* out, const psImage* in)
    292 {
    293     psElemType type;
    294     unsigned int numCols;
    295     unsigned int numRows;
    296 
    297     if (in == NULL) {
    298         psFree(out);
    299         return NULL;
    300     }
    301 
    302     type = in->type.type;
    303     numCols = in->numCols;
    304     numRows = in->numRows;
    305 
    306     /* if not a complex number, this is logically just a image copy */
    307     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    308         // Warn user, as this is probably not expected
    309         psLogMsg(__func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
    310                  "Image copy was performed instead.");
    311         return psImageCopy(out, in, type);
    312     }
    313 
    314     if (type == PS_TYPE_C32) {
    315         psC32* outRow;
    316         psC32* inRow;
    317 
    318         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32);
    319         for (unsigned int row = 0; row < numRows; row++) {
    320             outRow = out->data.C32[row];
    321             inRow = in->data.C32[row];
    322 
    323             for (unsigned int col = 0; col < numCols; col++) {
    324                 outRow[col] = crealf(inRow[col]) - I * cimagf(inRow[col]);
    325             }
    326         }
    327     } else if (type == PS_TYPE_C64) {
    328         psC64* outRow;
    329         psC64* inRow;
    330 
    331         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64);
    332         for (unsigned int row = 0; row < numRows; row++) {
    333             outRow = out->data.C64[row];
    334             inRow = in->data.C64[row];
    335 
    336             for (unsigned int col = 0; col < numCols; col++) {
    337                 outRow[col] = creal(inRow[col]) - I * cimag(inRow[col]);
    338             }
    339         }
    340     } else {
    341         psError(__func__, "Can not compute complex conjugate for given image type (%d).", type);
    342         psFree(out);
    343         return NULL;
    344     }
    345 
    346     return out;
    347 }
    348 
    349 psImage* psImagePowerSpectrum(psImage* out, const psImage* in)
    350 {
    351     psElemType type;
    352     unsigned int numCols;
    353     unsigned int numRows;
    354     int numElementsSquared;
    355 
    356     if (in == NULL) {
    357         psFree(out);
    358         return NULL;
    359     }
    360 
    361     type = in->type.type;
    362     numCols = in->numCols;
    363     numRows = in->numRows;
    364     numElementsSquared = numCols * numCols * numRows * numRows;
    365 
    366     /* if not a complex number, this is not implemented */
    367     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    368         psError(__func__, "Power Spectrum for non-complex inputs is not implemented.");
    369         psFree(out);
    370         return NULL;
    371     }
    372 
    373     if (type == PS_TYPE_C32) {
    374         psF32* outRow;
    375         psC32* inRow;
    376         psF32 real;
    377         psF32 imag;
    378 
    379         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
    380         for (unsigned int row = 0; row < numRows; row++) {
    381             outRow = out->data.F32[row];
    382             inRow = in->data.C32[row];
    383 
    384             for (unsigned int col = 0; col < numCols; col++) {
    385                 real = crealf(inRow[col]);
    386                 imag = cimagf(inRow[col]);
    387                 outRow[col] = (real * real + imag * imag) / numElementsSquared;
    388             }
    389         }
    390     } else if (type == PS_TYPE_C64) {
    391         psF64* outRow;
    392         psC64* inRow;
    393         psF64 real;
    394         psF64 imag;
    395 
    396         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
    397         for (unsigned int row = 0; row < numRows; row++) {
    398             outRow = out->data.F64[row];
    399             inRow = in->data.C64[row];
    400 
    401             for (unsigned int col = 0; col < numCols; col++) {
    402                 real = crealf(inRow[col]);
    403                 imag = cimagf(inRow[col]);
    404                 outRow[col] = real * real + imag * imag / numElementsSquared;
    405             }
    406         }
    407     } else {
    408         psError(__func__, "Can not power spectrum for given image type (%d).", type);
    409         psFree(out);
    410         return NULL;
    411     }
    412 
    413     return out;
    414 
    415 }
    416 
    417 /************************************** Vector Functions ***************************************/
    41828
    41929psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction)
Note: See TracChangeset for help on using the changeset viewer.