IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 7, 2007, 6:17:58 PM (19 years ago)
Author:
Paul Price
Message:

Reworked image convolution (split direct and FFT methods into separate functions) and FFT functions. psVectorFFT needs a bit more work to clean up the N/2 issues from FFTW's r2c and c2r.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r11668 r11703  
    1 /** @file  psImageFFT.c
    2  *
    3  *  @brief Contains FFT transform related functions for psImage.
    4  *
    5  *  @author Robert DeSonia, MHPCC
    6  *
    7  *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-02-06 21:36:09 $
    9  *
    10  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    11  */
     1/// @file  psImageFFT.c
     2///
     3/// @brief Contains FFT transform related functions for psImage.
     4///
     5/// @author Paul Price, IfA
     6/// @author Robert DeSonia, MHPCC
     7///
     8/// @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     9/// @date $Date: 2007-02-08 04:17:58 $
     10///
     11/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     12///
    1213
    1314#ifdef HAVE_CONFIG_H
    … …  
    2021#include <fftw3.h>
    2122
    22 #include "psImageFFT.h"
     23#include "psAssert.h"
    2324#include "psError.h"
    2425#include "psMemory.h"
    2526#include "psLogMsg.h"
     27#include "psConstants.h"
    2628#include "psImageStructManip.h"
    27 
    28 
    29 
    30 #define PS_FFTW_PLAN_RIGOR FFTW_ESTIMATE
    31 
    32 static bool p_fftwWisdomImported = false;
    33 
    34 psImage* psImageFFT(psImage* out, const psImage* image, psFFTFlags direction)
    35 {
    36     psU32 numCols;
    37     psU32 numRows;
    38     psElemType type;
    39     fftwf_plan plan;
    40 
    41     /* got good image data? */
    42     if (image == NULL) {
    43         psFree(out);
     29#include "psImageFFT.h"
     30
     31
     32#define FFTW_PLAN_RIGOR FFTW_ESTIMATE   // How rigorous the FFTW planning is
     33
     34static psBool fftwWisdomImported = false; // Has the FFTW wisdom been imported yet?
     35
     36bool psImageForwardFFT(psImage **real, psImage **imag, const psImage *in)
     37{
     38    PS_ASSERT_IMAGE_NON_NULL(in, false);
     39    PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, false);
     40    PS_ASSERT_PTR_NON_NULL(real, false);
     41    PS_ASSERT_PTR_NON_NULL(imag, false);
     42
     43    // Make sure the system-level wisdom information is imported.
     44    if (!fftwWisdomImported) {
     45        fftwf_import_system_wisdom();
     46        fftwWisdomImported = true;
     47    }
     48
     49    int numCols = in->numCols;          // Number of columns
     50    int numRows = in->numRows;          // Number of rows
     51
     52    psF32 *input;                       // Input data, for FFTW
     53    if (!in->parent) {
     54        // No parent --- can just use the data buffer
     55        input = psMemIncrRefCounter(in->p_rawDataBuffer);
     56    } else {
     57        // Need to copy the data
     58        input = psAlloc(numCols * numRows);
     59        for (int y = 0; y < numRows; y++) {
     60            memcpy(&input[y * numRows], in->data.F32[y], numCols);
     61        }
     62    }
     63
     64    // Do the FFT
     65    fftwf_complex *out = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // Output data
     66    fftwf_plan plan = fftwf_plan_dft_r2c_2d(numRows, numCols, input, out, FFTW_PLAN_RIGOR);
     67    fftwf_execute(plan);
     68    fftwf_destroy_plan(plan);
     69    psFree(input);
     70
     71    // Pull the real and imaginary parts out
     72    numCols = numCols/2 + 1;            // x dimension is halved by FFTW
     73    *real = psImageRecycle(*real, numCols, numRows, PS_TYPE_F32);
     74    *imag = psImageRecycle(*imag, numCols, numRows, PS_TYPE_F32);
     75    for (int y = 0, index = 0; y < numRows; y++) {
     76        for (int x = 0; x < numCols; index++, x++) {
     77#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
     78            // C99 complex support
     79            (*real)->data.F32[y][x] = creal(out[index]);
     80            (*imag)->data.F32[y][x] = cimag(out[index]);
     81#else
     82            // FFTW's backup complex support
     83            (*real)->data.F32[y][x] = out[index][0];
     84            (*imag)->data.F32[y][x] = out[index][1];
     85#endif
     86        }
     87    }
     88    fftwf_free(out);
     89
     90    return true;
     91}
     92
     93bool psImageBackwardFFT(psImage **out, const psImage *real, const psImage *imag, int origCols)
     94{
     95    PS_ASSERT_IMAGE_NON_NULL(real, false);
     96    PS_ASSERT_IMAGE_TYPE(real, PS_TYPE_F32, false);
     97    PS_ASSERT_IMAGE_NON_NULL(imag, false);
     98    PS_ASSERT_IMAGE_TYPE(imag, PS_TYPE_F32, false);
     99    PS_ASSERT_IMAGES_SIZE_EQUAL(real, imag, false);
     100    PS_ASSERT_PTR_NON_NULL(out, false);
     101
     102    int numCols = real->numCols;        // Number of columns
     103    int numRows = real->numRows;        // Number of rows
     104
     105    // Because of the way FFT r2c and c2r work, need the number of columns in the target to be:
     106    // 2 * numCols = 2*(numCols/2 + 1) = origCols % 2 ? origCols + 1 : origCols + 2
     107    if (2 * numCols != ((origCols % 2) ? origCols + 1 : origCols + 2)) {
     108        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     109                "Number of columns in FFT-ed images (%d) should be origCols (%d) / 2 + 1 = %d",
     110                numCols, origCols, origCols/2 + 1);
     111        return false;
     112    }
     113
     114    if (*out && (*out)->parent) {
     115        // It has a parent, so we can't write directly into the buffer.
     116        // It had better be the correct size and type, because we don't want to resize a child image
     117        PS_ASSERT_IMAGE_SIZE(*out, origCols, numRows, false);
     118        PS_ASSERT_IMAGE_TYPE(*out, PS_TYPE_F32, false);
     119    }
     120
     121
     122    // Make sure the system-level wisdom information is imported.
     123    if (!fftwWisdomImported) {
     124        fftwf_import_system_wisdom();
     125        fftwWisdomImported = true;
     126    }
     127
     128    // Stuff the real and imaginary parts in
     129    psF32 *target = psAlloc(2 * numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Target for FFTW
     130    fftwf_complex *in = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // Input data
     131    for (int y = 0, index = 0; y < numRows; y++) {
     132        for (int x = 0; x < numCols; index++, x++) {
     133#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
     134            // C99 complex support
     135            in[index] = real->data.F32[y][x] + imag->data.F32[y][x] * I;
     136#else
     137            // FFTW's backup complex support
     138            in[index][0] = real->data.F32[y][x];
     139            in[index][1] = imag->data.F32[y][x];
     140#endif
     141        }
     142    }
     143
     144    // Do the FFT
     145    fftwf_plan plan = fftwf_plan_dft_c2r_2d(numRows, origCols, in, target, FFTW_PLAN_RIGOR);
     146    fftwf_execute(plan);
     147    fftwf_destroy_plan(plan);
     148    fftwf_free(in);
     149
     150    if (!(*out) || !(*out)->parent) {
     151        *out = psImageRecycle(*out, origCols, numRows, PS_TYPE_F32);
     152    }
     153
     154    // Copy the target pixels into the output
     155    // There's a slight offset --- target has the additional padding required by FFTW,
     156    // while the output doesn't.
     157    for (int y = 0, index = 0; y < numRows; y++, index += origCols) {
     158        memcpy(&(*out)->data.F32[y][0], &target[index], numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
     159    }
     160    psFree(target);
     161
     162    return true;
     163}
     164
     165psImage* psImagePowerSpectrum(psImage *out, const psImage *in)
     166{
     167    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
     168    PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, NULL);
     169
     170    psImage *real = NULL;              // Real component of FFT
     171    psImage *imag = NULL;              // Imaginary component of FFT
     172
     173    if (!psImageForwardFFT(&real, &imag, in)) {
     174        psError(PS_ERR_UNKNOWN, false, "Unable to perform forward FFT.");
    44175        return NULL;
    45176    }
    46177
    47     if ( ((direction & PS_FFT_FORWARD) != 0) ) {
    48         if ((direction & PS_FFT_REVERSE) != 0) {
    49             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    50                     _("Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options."));
    51             psFree(out);
    52             return NULL;
    53         }
    54         if ((direction & PS_FFT_REAL_RESULT) != 0) {
    55             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    56                     _("The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported."));
    57             psFree(out);
    58             return NULL;
    59         }
    60     } else if ((direction & PS_FFT_REVERSE) == 0) {
     178    int numCols = real->numCols;        // Number of columns
     179    int numRows = real->numRows;        // Number of rows
     180
     181    float norm = 1.0 / PS_SQR(in->numCols) / PS_SQR(in->numRows);
     182    for (int y = 0; y < numRows; y++) {
     183        for (int x = 0; x < numCols; x++) {
     184            // Power spectrum is the square of the complex modulus
     185            real->data.F32[y][x] = norm * (PS_SQR(real->data.F32[y][x]) + PS_SQR(imag->data.F32[y][x]));
     186        }
     187    }
     188    psFree(imag);
     189
     190    return real;
     191}
     192
     193
     194
     195bool psImageComplexMultiply(psImage **outReal, psImage **outImag,
     196                            const psImage *in1Real, const psImage *in1Imag,
     197                            const psImage *in2Real, const psImage *in2Imag)
     198{
     199    PS_ASSERT_IMAGE_NON_NULL(in1Real, false);
     200    PS_ASSERT_IMAGE_NON_NULL(in1Imag, false);
     201    PS_ASSERT_IMAGE_NON_NULL(in2Real, false);
     202    PS_ASSERT_IMAGE_NON_NULL(in2Imag, false);
     203    PS_ASSERT_IMAGE_TYPE(in1Real, PS_TYPE_F32, false);
     204    PS_ASSERT_IMAGE_TYPE(in1Imag, PS_TYPE_F32, false);
     205    PS_ASSERT_IMAGE_TYPE(in2Real, PS_TYPE_F32, false);
     206    PS_ASSERT_IMAGE_TYPE(in2Imag, PS_TYPE_F32, false);
     207    PS_ASSERT_IMAGES_SIZE_EQUAL(in1Imag, in1Real, false);
     208    PS_ASSERT_IMAGES_SIZE_EQUAL(in2Real, in1Real, false);
     209    PS_ASSERT_IMAGES_SIZE_EQUAL(in2Imag, in1Real, false);
     210    PS_ASSERT_PTR_NON_NULL(outReal, false);
     211    PS_ASSERT_PTR_NON_NULL(outImag, false);
     212    if (*outReal) {
     213        PS_ASSERT_IMAGE_NON_NULL(*outReal, false);
     214        PS_ASSERT_IMAGE_NON_NULL(*outImag, false);
     215    } else if (*outImag) {
    61216        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    62                 _("Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option."));
    63         psFree(out);
    64         return NULL;
    65     }
    66 
    67     type = image->type.type;
    68 
    69     /* make sure the system-level wisdom information is imported. */
    70     if (!p_fftwWisdomImported) {
    71         fftwf_import_system_wisdom();
    72         p_fftwWisdomImported = true;
    73     }
    74 
    75     numRows = image->numRows;
    76     numCols = image->numCols;
    77 
    78     // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place.
    79     psS32 sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD;
    80 
    81     fftwf_complex* outBuffer = fftwf_malloc(numRows*numCols*sizeof(fftwf_complex));
    82     p_psImageCopyToRawBuffer(outBuffer, image, PS_TYPE_C32);
    83 
    84     plan = fftwf_plan_dft_2d(numRows, numCols,
    85                              outBuffer,
    86                              outBuffer,
    87                              sign,
    88                              PS_FFTW_PLAN_RIGOR);
    89 
    90     /* check if a plan exists now -- if not, it is a real problem at this point */
    91     if (plan == NULL) {
    92         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    93                 _("Could not create a valid FFT plan to perform the transform."));
    94         psFree(out);
    95         return NULL;
    96     }
    97 
    98     /* finally, call FFTW with the plan made above */
    99     fftwf_execute(plan);
    100 
    101     fftwf_destroy_plan(plan);
    102 
    103     if (direction & PS_FFT_REAL_RESULT) {
    104         // n.b., we do this instead of using fftwf_plan_dft_c2r because that
    105         // plan requires a half-image, which would require a image reordering
    106         // that is not as simple as performing a normal complex transform and
    107         // then taking the real part of the result.  If performance here
    108         // becomes an issue, the use of fftwf_plan_dft_r2c should be considered
    109         // as well as fftwf_plan_dft_c2r.
    110         out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
    111         int index = 0;
    112         for (int row=0; row < numRows; row++) {
    113             psF32* outRow = out->data.F32[row];
    114             for (int col=0; col < numCols; col++) {
    115                 outRow[col] = crealf(outBuffer[index++]); // take just the real part
    116             }
    117         }
     217                "If the output real part is provided, the output imaginary part must also be provided.");
     218        return false;
     219    }
     220
     221    int numRows = in1Real->numRows;     // Number of rows
     222    int numCols = in1Real->numCols;     // Number of columns
     223
     224    // Need to worry if the outputs are children
     225    psImage *targetReal, *targetImag;   // Target real and imaginary parts
     226
     227    if ((*outReal)->parent) {
     228        // It had better be the correct size and type, because we don't want to resize a child image
     229        PS_ASSERT_IMAGE_TYPE(*outReal, PS_TYPE_F32, false);
     230        PS_ASSERT_IMAGE_SIZE(*outReal, numCols, numRows, false);
     231        targetReal = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    118232    } else {
    119         out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
    120         int index = 0;
    121         for (int row=0; row < numRows; row++) {
    122             psC32* outRow = out->data.C32[row];
    123             for (int col=0; col < numCols; col++) {
    124                 outRow[col] = outBuffer[index++]; // take just the real part
    125             }
    126         }
    127         //        memcpy(out->rawDataBuffer, outBuffer, numRows*numCols*sizeof(fftwf_complex));
    128     }
    129 
    130     fftwf_free(outBuffer);
    131 
    132     return out;
    133 
    134 }
    135 
    136 psImage* psImageReal(psImage* out, const psImage* in)
    137 {
    138     psElemType type;
    139     psU32 numCols;
    140     psU32 numRows;
    141 
    142     if (in == NULL) {
    143         psFree(out);
    144         return NULL;
    145     }
    146 
    147     type = in->type.type;
    148     numCols = in->numCols;
    149     numRows = in->numRows;
    150 
    151     /* if not a complex number, this is logically just a copy then */
    152     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    153         return psImageCopy(out, in, type);
    154     }
    155 
    156     if (type == PS_TYPE_C32) {
    157         psF32* outRow;
    158         psC32* inRow;
    159 
    160         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
    161         for (psU32 row = 0; row < numRows; row++) {
    162             outRow = out->data.F32[row];
    163             inRow = in->data.C32[row];
    164 
    165             for (psU32 col = 0; col < numCols; col++) {
    166                 outRow[col] = crealf(inRow[col]);
    167             }
    168         }
    169     } else if (type == PS_TYPE_C64) {
    170         psF64* outRow;
    171         psC64* inRow;
    172 
    173         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
    174         for (psU32 row = 0; row < numRows; row++) {
    175             outRow = out->data.F64[row];
    176             inRow = in->data.C64[row];
    177 
    178             for (psU32 col = 0; col < numCols; col++) {
    179                 outRow[col] = creal(inRow[col]);
    180             }
    181         }
     233        *outReal = psImageRecycle(*outReal, numCols, numRows, PS_TYPE_F32);
     234        targetReal = psMemIncrRefCounter(*outReal);
     235    }
     236    if ((*outImag)->parent) {
     237        // It had better be the correct size and type, because we don't want to resize a child image
     238        if ((*outReal)->numCols != numCols || (*outReal)->numRows != numRows ||
     239            (*outImag)->type.type != PS_TYPE_F32) {
     240            // Plug potential memory leak --- need to free targetReal if there's a problem
     241            psFree(targetReal);
     242        }
     243        PS_ASSERT_IMAGE_TYPE(*outReal, PS_TYPE_F32, false);
     244        PS_ASSERT_IMAGE_SIZE(*outReal, numCols, numRows, false);
     245        targetImag = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    182246    } else {
    183         char* typeStr;
    184         PS_TYPE_NAME(typeStr,type);
    185         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    186                 _("Specified psImage type, %s, is not supported."),
    187                 typeStr);
    188 
    189         psFree(out);
    190         return NULL;
    191     }
    192 
    193     return out;
    194 }
    195 
    196 psImage* psImageImaginary(psImage* out, const psImage* in)
    197 {
    198     psElemType type;
    199     psU32 numCols;
    200     psU32 numRows;
    201 
    202     if (in == NULL) {
    203         psFree(out);
    204         return NULL;
    205     }
    206 
    207     type = in->type.type;
    208     numCols = in->numCols;
    209     numRows = in->numRows;
    210 
    211     /* if not a complex image type, this is logically just zeroed image of same size */
    212     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    213         out = psImageRecycle(out, numCols, numRows, type);
    214         memset(out->data.V[0], 0, PSELEMTYPE_SIZEOF(type) * numCols * numRows);
    215         return out;
    216     }
    217 
    218     if (type == PS_TYPE_C32) {
    219         psF32* outRow;
    220         psC32* inRow;
    221 
    222         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
    223         for (psU32 row = 0; row < numRows; row++) {
    224             outRow = out->data.F32[row];
    225             inRow = in->data.C32[row];
    226 
    227             for (psU32 col = 0; col < numCols; col++) {
    228                 outRow[col] = cimagf(inRow[col]);
    229             }
    230         }
    231     } else if (type == PS_TYPE_C64) {
    232         psF64* outRow;
    233         psC64* inRow;
    234 
    235         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
    236         for (psU32 row = 0; row < numRows; row++) {
    237             outRow = out->data.F64[row];
    238             inRow = in->data.C64[row];
    239 
    240             for (psU32 col = 0; col < numCols; col++) {
    241                 outRow[col] = cimag(inRow[col]);
    242             }
    243         }
    244     } else {
    245         char* typeStr;
    246         PS_TYPE_NAME(typeStr,type);
    247         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    248                 _("Specified psImage type, %s, is not supported."),
    249                 typeStr);
    250         psFree(out);
    251         return NULL;
    252     }
    253 
    254     return out;
    255 }
    256 
    257 psImage* psImageComplex(psImage* out, const psImage* real, const psImage* imag)
    258 {
    259     psElemType type;
    260     psU32 numCols;
    261     psU32 numRows;
    262 
    263     if (real == NULL || imag == NULL) {
    264         psFree(out);
    265         return NULL;
    266     }
    267 
    268     type = real->type.type;
    269     numCols = real->numCols;
    270     numRows = real->numRows;
    271 
    272     if (imag->type.type != type) {
    273         char* typeStrReal;
    274         char* typeStrImag;
    275         PS_TYPE_NAME(typeStrReal,type);
    276         PS_TYPE_NAME(typeStrImag,imag->type.type);
    277         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    278                 _("Real psImage type (%s) and imaginary psImage type (%s) must be the same."),
    279                 typeStrReal,typeStrImag);
    280         psFree(out);
    281         return NULL;
    282     }
    283 
    284     if (imag->numCols != numCols || imag->numRows != numRows) {
    285         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    286                 _("Real psImage size (%dx%d) and imaginary psImage size (%dx%d) must be the same."),
    287                 numCols, numRows, imag->numCols, imag->numRows);
    288         psFree(out);
    289         return NULL;
    290     }
    291 
    292     if (type == PS_TYPE_F32) {
    293         psC32* outRow;
    294         psF32* realRow;
    295         psF32* imagRow;
    296 
    297         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32);
    298 
    299         for (psU32 row = 0; row < numRows; row++) {
    300             outRow = out->data.C32[row];
    301             realRow = real->data.F32[row];
    302             imagRow = imag->data.F32[row];
    303 
    304             for (psU32 col = 0; col < numCols; col++) {
    305                 outRow[col] = realRow[col] + I * imagRow[col];
    306             }
    307         }
    308     } else if (type == PS_TYPE_F64) {
    309         psC64* outRow;
    310         psF64* realRow;
    311         psF64* imagRow;
    312 
    313         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64);
    314         for (psU32 row = 0; row < numRows; row++) {
    315             outRow = out->data.C64[row];
    316             realRow = real->data.F64[row];
    317             imagRow = imag->data.F64[row];
    318 
    319             for (psU32 col = 0; col < numCols; col++) {
    320                 outRow[col] = realRow[col] + I * imagRow[col];
    321             }
    322         }
    323     } else {
    324         char* typeStr;
    325         PS_TYPE_NAME(typeStr,type);
    326         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    327                 _("Input psImage type, %s, is required to be either psF32 or psF64."),
    328                 typeStr);
    329         psFree(out);
    330         return NULL;
    331     }
    332 
    333 
    334     return out;
    335 }
    336 
    337 psImage* psImageConjugate(psImage* out, const psImage* in)
    338 {
    339     psElemType type;
    340     psU32 numCols;
    341     psU32 numRows;
    342 
    343     if (in == NULL) {
    344         psFree(out);
    345         return NULL;
    346     }
    347 
    348     type = in->type.type;
    349     numCols = in->numCols;
    350     numRows = in->numRows;
    351 
    352     /* if not a complex image, this is logically just a image copy */
    353     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    354         return psImageCopy(out, in, type);
    355     }
    356 
    357     if (type == PS_TYPE_C32) {
    358         psC32* outRow;
    359         psC32* inRow;
    360 
    361         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32);
    362         for (psU32 row = 0; row < numRows; row++) {
    363             outRow = out->data.C32[row];
    364             inRow = in->data.C32[row];
    365 
    366             for (psU32 col = 0; col < numCols; col++) {
    367                 outRow[col] = crealf(inRow[col]) - I * cimagf(inRow[col]);
    368             }
    369         }
    370     } else if (type == PS_TYPE_C64) {
    371         psC64* outRow;
    372         psC64* inRow;
    373 
    374         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64);
    375         for (psU32 row = 0; row < numRows; row++) {
    376             outRow = out->data.C64[row];
    377             inRow = in->data.C64[row];
    378 
    379             for (psU32 col = 0; col < numCols; col++) {
    380                 outRow[col] = creal(inRow[col]) - I * cimag(inRow[col]);
    381             }
    382         }
    383     } else {
    384         char* typeStr;
    385         PS_TYPE_NAME(typeStr,type);
    386         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    387                 _("Input psImage type, %s, is required to be either psC32 or psC64."),
    388                 typeStr);
    389         psFree(out);
    390         return NULL;
    391     }
    392 
    393     return out;
    394 }
    395 
    396 psImage* psImagePowerSpectrum(psImage* out, const psImage* in)
    397 {
    398     psElemType type;
    399     psU32 numCols;
    400     psU32 numRows;
    401 
    402     if (in == NULL) {
    403         psFree(out);
    404         return NULL;
    405     }
    406 
    407     type = in->type.type;
    408     numCols = in->numCols;
    409     numRows = in->numRows;
    410 
    411     if (type == PS_TYPE_C32) {
    412         psF32* outRow;
    413         psC32* inRow;
    414         psF32 real;
    415         psF32 imag;
    416         psF32 fNumCols = numCols;
    417         psF32 fNumRows = numRows;
    418         psF32 numElementsSquared = fNumCols * fNumCols * fNumRows * fNumRows;
    419 
    420         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
    421         for (psU32 row = 0; row < numRows; row++) {
    422             outRow = out->data.F32[row];
    423             inRow = in->data.C32[row];
    424 
    425             for (psU32 col = 0; col < numCols; col++) {
    426                 real = crealf(inRow[col]);
    427                 imag = cimagf(inRow[col]);
    428                 outRow[col] = (real * real + imag * imag) / numElementsSquared;
    429             }
    430         }
    431     } else if (type == PS_TYPE_C64) {
    432         psF64* outRow;
    433         psC64* inRow;
    434         psF64 real;
    435         psF64 imag;
    436         psF64 numElementsSquared = numCols * numCols * numRows * numRows;
    437 
    438         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
    439         for (psU32 row = 0; row < numRows; row++) {
    440             outRow = out->data.F64[row];
    441             inRow = in->data.C64[row];
    442 
    443             for (psU32 col = 0; col < numCols; col++) {
    444                 real = creal(inRow[col]);
    445                 imag = cimag(inRow[col]);
    446                 outRow[col] = (real * real + imag * imag) / numElementsSquared;
    447             }
    448         }
    449     } else {
    450         char* typeStr;
    451         PS_TYPE_NAME(typeStr,type);
    452         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    453                 _("Input psImage type, %s, is required to be either psC32 or psC64."),
    454                 typeStr);
    455         psFree(out);
    456         return NULL;
    457     }
    458 
    459     return out;
    460 
    461 }
     247        *outImag = psImageRecycle(*outImag, numCols, numRows, PS_TYPE_F32);
     248        targetImag = psMemIncrRefCounter(*outImag);
     249    }
     250
     251    for (int y = 0; y < numRows; y++) {
     252        for (int x = 0; x < numCols; x++) {
     253            // (a + bi) * (c + di) = (ac - bd) + (bc + ad)i
     254            float real = in1Real->data.F32[y][x] * in2Real->data.F32[y][x] -
     255                in1Imag->data.F32[y][x] * in2Imag->data.F32[y][x];
     256            float imag = in1Imag->data.F32[y][x] * in2Real->data.F32[y][x] +
     257                in1Real->data.F32[y][x] * in2Imag->data.F32[y][x];
     258            targetReal->data.F32[y][x] = real;
     259            targetImag->data.F32[y][x] = imag;
     260        }
     261    }
     262
     263    if ((*outReal)->parent) {
     264        *outReal = psImageCopy(*outReal, targetReal, PS_TYPE_F32);
     265    }
     266    if ((*outImag)->parent) {
     267        *outImag = psImageCopy(*outImag, targetImag, PS_TYPE_F32);
     268    }
     269
     270    psFree(targetReal);
     271    psFree(targetImag);
     272
     273    return true;
     274}
Note: See TracChangeset for help on using the changeset viewer.