IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 11703


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.

Location:
trunk/psLib
Files:
1 added
9 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}

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

    r11248 r11703  
    1 /* @file  psImageFFT.h
    2  * @brief Contains FFT transform related functions for psImage
    3  *
    4  * @author Robert DeSonia, MHPCC
    5  *
    6  * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2007-01-23 22:47:22 $
    8  * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    9  */
     1/// @file  psImageFFT.h
     2/// @brief Contains FFT transform related functions for psImage
     3///
     4/// @author Paul Price, IfA
     5/// @author Robert DeSonia, MHPCC
     6///
     7/// @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     8/// @date $Date: 2007-02-08 04:17:58 $
     9/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     10///
    1011
    1112#ifndef PS_IMAGE_FFT_H
    … …  
    1314
    1415#include "psImage.h"
    15 #include "psVectorFFT.h"               // for psFFTFlags
    1616
    1717/// @addtogroup MathOps Mathematical Operations
    1818/// @{
    1919
    20 /** Forward and reverse FFT calculations.
    21  *
    22  *  This takes as input the image of interest (in) and the direction
    23  *  (direction), which is specified by an enumerated type psFftDirection.
    24  *  The input image may be of type psF32 or psC32, the result is always
    25  *  psC32. If the input vector is psF32, the direction must be forward.
    26  *
    27  *  @return psImage* the FFT transformation result
    28  */
    29 psImage* psImageFFT(
    30     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    31     const psImage* image,              ///< the psImage to apply transform to
    32     psFFTFlags direction               ///< the direction of the transform
    33 );
     20/// Forward FFT of an image
     21///
     22/// Applies a forward FFT (exponent -1), with the result returned in both real and imaginary parts.  The FFT
     23/// is not normalised (a forward followed by a reverse is the original scaled by the total number of pixels).
     24/// The FFT takes advantage of the fact that the input is purely real; hence the output number of columns is
     25/// numCols/2 + 1 (with division rounding down).  Only implemented for F32 input.
     26bool psImageForwardFFT(psImage **real,  ///< Real part of FFT
     27                       psImage **imag,  ///< Imaginary part of FFT
     28                       const psImage *in///< Input image (F32)
     29    );
    3430
    35 /** extract the real portion of a complex image
    36  *
    37  *  @return psImage*   real portion of the input image.
    38  */
    39 psImage* psImageReal(
    40     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    41     const psImage* in                  ///< the psImage to extract real portion from
    42 );
     31/// Backward FFT of an image
     32///
     33/// Applies a backward FFT (exponent +1) from the real and imaginary parts with the (purely) real result
     34/// returned.  The FFT is not normalised (a forward followed by a reverse is the original scaled by the total
     35/// number of pixels).  The FFT takes advantage of the fact that the output will be purely real; hence the
     36/// input number of columns is numCols/2 + 1 (with division rounding down); to manage the redundancy (is the
     37/// original number of columns even or odd?), we need the original number of columns (the number of columns of
     38/// the image that was input to psImageForwardFFT) to be provided.  Only implemented for F32 input.
     39bool psImageBackwardFFT(psImage **out,///< Output image
     40                        const psImage *real, ///< Real input (F32)
     41                        const psImage *imag, ///< Imaginary input (F32)
     42                        int origCols    ///< Original number of columns
     43    );
    4344
    44 /** extract the imaginary portion of a complex image
    45  *
    46  *  @return psImage*   imaginary portion of the input image.
    47  */
    48 psImage* psImageImaginary(
    49     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    50     const psImage* in                  ///< the psImage to extract imaginary portion from
    51 );
    5245
    53 /** creates a complex image from separate real and imaginary plane images
    54  *
    55  *  @return psImage*   resulting complex image
    56  */
    57 psImage* psImageComplex(
    58     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    59     const psImage* real,               ///< the real plane image
    60     const psImage* imag                ///< the imaginary plane image
    61 );
     46/// Power spectrum of an image
     47///
     48/// Generates the power spectrum of an image.  Only implemented for F32 input.
     49psImage* psImagePowerSpectrum(psImage *out, const psImage *in);
    6250
    63 /** computes the complex conjugate of an image
    64  *
    65  *  @return psImage*   the complex conjugate of the 'in' image
    66  */
    67 psImage* psImageConjugate(
    68     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    69     const psImage* in                  ///< the psImage to compute conjugate of
    70 );
    71 
    72 /** computes the power spectrum of an image
    73  *
    74  *  @return psImage*   the power spectrum of the 'in' image
    75  */
    76 psImage* psImagePowerSpectrum(
    77     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    78     const psImage* in                  ///< the psImage to power spectrum of
    79 );
     51/// Multiply complex images
     52///
     53/// The input images are the real and imaginary parts of each of two images.  The real and imaginary parts
     54/// of the output image are returned.  Only implemented for F32 input.
     55bool psImageComplexMultiply(psImage **outReal, ///< Real part of output
     56                            psimage **outImag, ///< Imaginary part of output
     57                            const psImage *in1Real, ///< Real part of input 1
     58                            const psImage *in1Imag, ///< Imaginary part of input 1
     59                            const psImage *in2Real, ///< Real part of input 2
     60                            const psImage *in2Imag ///< Imaginary part of input 2
     61    );
    8062
    8163/// @}

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

    r11668 r11703  
    33 *  @brief Contains FFT transform related functions for psVector
    44 *
     5 *  @author Paul Price, IfA
    56 *  @author Robert DeSonia, MHPCC
    67 *
    7  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-02-06 21:36:09 $
     8 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-02-08 04:17:58 $
    910 *
    1011 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    1213
    1314#ifdef HAVE_CONFIG_H
    14 # include "config.h"
     15#include "config.h"
    1516#endif
    1617
    … …  
    2122#include <fftw3.h>
    2223
    23 #include "psVectorFFT.h"
     24#include "psAssert.h"
    2425#include "psError.h"
    2526#include "psMemory.h"
    2627#include "psLogMsg.h"
     28#include "psConstants.h"
     29#include "psVectorFFT.h"
     30
     31#define FFTW_PLAN_RIGOR FFTW_ESTIMATE   // How rigorous the FFTW planning is
     32
     33static psBool fftwWisdomImported = false; // Has the system wisdom been imported?
     34
     35bool psVectorForwardFFT(psVector **real, psVector **imag, const psVector *in)
     36{
     37    PS_ASSERT_VECTOR_NON_NULL(in, false);
     38    PS_ASSERT_VECTOR_TYPE(in, PS_TYPE_F32, false);
     39    PS_ASSERT_PTR_NON_NULL(real, false);
     40    PS_ASSERT_PTR_NON_NULL(imag, false);
     41
     42    // Make sure the system-level wisdom information is imported.
     43    if (!fftwWisdomImported) {
     44        fftwf_import_system_wisdom();
     45        fftwWisdomImported = true;
     46    }
     47
     48    long num = in->n;                   // Number of elements
     49
     50    // Do the FFT
     51    fftwf_complex *out = fftwf_malloc(num * sizeof(fftwf_complex)); // Output data
     52    fftwf_plan plan = fftwf_plan_dft_r2c_1d(num, in->data.F32, out, FFTW_PLAN_RIGOR);
     53    fftwf_execute(plan);
     54    fftwf_destroy_plan(plan);
     55
     56    // Pull the real and imaginary parts out
     57    *real = psVectorRecycle(*real, num/2 + 1, PS_TYPE_F32);
     58    *imag = psVectorRecycle(*imag, num/2 + 1, PS_TYPE_F32);
     59    for (int i = 0; i < num; i++) {
     60#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
     61        // C99 complex support
     62        (*real)->data.F32[i] = creal(out[i]);
     63        (*imag)->data.F32[i] = cimag(out[i]);
     64#else
     65        // FFTW's backup complex support
     66        (*real)->data.F32[i] = out[i][0];
     67        (*imag)->data.F32[i] = out[i][1];
     68#endif
     69    }
     70    fftwf_free(out);
     71
     72    return true;
     73}
     74
     75bool psVectorBackwardFFT(psVector **out, const psVector *real, const psVector *imag)
     76{
     77    PS_ASSERT_VECTOR_NON_NULL(real, false);
     78    PS_ASSERT_VECTOR_TYPE(real, PS_TYPE_F32, false);
     79    PS_ASSERT_VECTOR_NON_NULL(imag, false);
     80    PS_ASSERT_VECTOR_TYPE(imag, PS_TYPE_F32, false);
     81    PS_ASSERT_VECTORS_SIZE_EQUAL(real, imag, false);
     82    PS_ASSERT_PTR_NON_NULL(out, false);
     83
     84    // Make sure the system-level wisdom information is imported.
     85    if (!fftwWisdomImported) {
     86        fftwf_import_system_wisdom();
     87        fftwWisdomImported = true;
     88    }
     89
     90    long num = real->n;                 // Number of elements
     91
     92    // Stuff the real and imaginary parts in
     93    fftwf_complex *in = fftwf_malloc(num * sizeof(fftwf_complex)); // Input data
     94    for (int i = 0; i < num; i++) {
     95#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
     96        // C99 complex support
     97        in[i] = real->data.F32[i] + imag->data.F32[i] * I;
     98#else
     99        // FFTW's backup complex support
     100        in[i][0] = real->data.F32[i];
     101        in[i][1] = imag->data.F32[i];
     102#endif
     103    }
    27104
    28105
     106    // Do the FFT
     107    *out = psVectorRecycle(*out, 2 * num, PS_TYPE_F32);
     108    fftwf_plan plan = fftwf_plan_dft_c2r_1d(2 * num, in, (*out)->data.F32, FFTW_PLAN_RIGOR);
     109    fftwf_execute(plan);
     110    fftwf_destroy_plan(plan);
    29111
    30 #define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE
     112    fftwf_free(in);
    31113
    32 static bool p_fftwWisdomImported = false;
     114    return true;
     115}
    33116
    34 psVector* psVectorFFT(psVector* out, const psVector* in, psFFTFlags direction)
     117psVector *psVectorPowerSpectrum(psVector* out, const psVector* in)
    35118{
    36     psU32 numElements;
    37     psElemType type;
    38     fftwf_plan plan;
     119    PS_ASSERT_VECTOR_NON_NULL(in, NULL);
     120    PS_ASSERT_VECTOR_TYPE(in, PS_TYPE_F32, NULL);
    39121
    40     /* got good image data? */
    41     if (in == NULL) {
    42         psFree(out);
     122    psVector *real = NULL;              // Real component of FFT
     123    psVector *imag = NULL;              // Imaginary component of FFT
     124
     125    if (!psVectorForwardFFT(&real, &imag, in)) {
     126        psError(PS_ERR_UNKNOWN, false, "Unable to perform forward FFT.");
    43127        return NULL;
    44128    }
    45129
    46     type = in->type.type;
     130    int num = real->n;                  // Number of elements
    47131
    48     /* make sure the system-level wisdom information is imported. */
    49     if (!p_fftwWisdomImported) {
    50         fftwf_import_system_wisdom();
    51         p_fftwWisdomImported = true;
     132    float norm = 1.0 / PS_SQR(in->n);
     133    for (int i = 0; i < num; i++) {
     134        // Power spectrum is the square of the complex modulus
     135        real->data.F32[i] = norm * (PS_SQR(real->data.F32[i]) + PS_SQR(imag->data.F32[i]));
     136    }
     137    psFree(imag);
     138
     139    return real;
     140}
     141
     142bool psVectorComplexMultiply(psVector **outReal, psVector **outImag,
     143                             const psVector *in1Real, const psVector *in1Imag,
     144                             const psVector *in2Real, const psVector *in2Imag)
     145{
     146    PS_ASSERT_VECTOR_NON_NULL(in1Real, false);
     147    PS_ASSERT_VECTOR_NON_NULL(in1Imag, false);
     148    PS_ASSERT_VECTOR_NON_NULL(in2Real, false);
     149    PS_ASSERT_VECTOR_NON_NULL(in2Imag, false);
     150    PS_ASSERT_VECTOR_TYPE(in1Real, PS_TYPE_F32, false);
     151    PS_ASSERT_VECTOR_TYPE(in1Imag, PS_TYPE_F32, false);
     152    PS_ASSERT_VECTOR_TYPE(in2Real, PS_TYPE_F32, false);
     153    PS_ASSERT_VECTOR_TYPE(in2Imag, PS_TYPE_F32, false);
     154    PS_ASSERT_VECTORS_SIZE_EQUAL(in1Imag, in1Real, false);
     155    PS_ASSERT_VECTORS_SIZE_EQUAL(in2Real, in1Real, false);
     156    PS_ASSERT_VECTORS_SIZE_EQUAL(in2Imag, in1Real, false);
     157    PS_ASSERT_PTR_NON_NULL(outReal, false);
     158    PS_ASSERT_PTR_NON_NULL(outImag, false);
     159    if (*outReal) {
     160        PS_ASSERT_VECTOR_NON_NULL(*outReal, false);
     161        PS_ASSERT_VECTOR_NON_NULL(*outImag, false);
     162    } else if (*outImag) {
     163        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     164                "If the output real part is provided, the output imaginary part must also be provided.");
     165        return false;
    52166    }
    53167
    54     numElements = in->n;
     168    int num = in1Real->n;               // Number of elements
    55169
    56     out = psVectorCopy(out, in, PS_TYPE_C32);
    57     out->n = numElements;
     170    *outReal = psVectorRecycle(*outReal, num, PS_TYPE_F32);
     171    *outImag = psVectorRecycle(*outImag, num, PS_TYPE_F32);
    58172
    59     if ((direction & PS_FFT_FORWARD) != 0) {
    60         plan = fftwf_plan_dft_1d(numElements,
    61                                  (fftwf_complex *) out->data.C32,
    62                                  (fftwf_complex *) out->data.C32, FFTW_FORWARD, P_FFTW_PLAN_RIGOR);
    63     } else if ((direction & PS_FFT_REVERSE) != 0) {
    64         plan = fftwf_plan_dft_1d(numElements,
    65                                  (fftwf_complex *) out->data.C32,
    66                                  (fftwf_complex *) out->data.C32, FFTW_BACKWARD, P_FFTW_PLAN_RIGOR);
    67     } else {
    68         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    69                 _("Must specify the direction as either PS_FFT_FORWARD or PS_FFT_REVERSE."));
    70         psFree(out);
    71         return NULL;
     173    for (int i = 0; i < num; i++) {
     174        // (a + bi) * (c + di) = (ac - bd) + (bc + ad)i
     175        (*outReal)->data.F32[i] = in1Real->data.F32[i] * in2Real->data.F32[i] -
     176            in1Imag->data.F32[i] * in2Imag->data.F32[i];
     177        (*outImag)->data.F32[i] = in1Imag->data.F32[i] * in2Real->data.F32[i] +
     178            in1Real->data.F32[i] * in2Imag->data.F32[i];
    72179    }
    73180
    74     /* check if a plan exists now */
    75     if (plan == NULL) {
    76         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    77                 _("Could not create a valid FFT plan to perform the transform."));
    78         psFree(out);
    79         return NULL;
    80     }
    81 
    82     /* finally, call FFTW with the plan made above */
    83     fftwf_execute(plan);
    84 
    85     fftwf_destroy_plan(plan);
    86 
    87     if ((direction & PS_FFT_REAL_RESULT) != 0) {
    88         for (psS32 i = 0; i < numElements; i++) {
    89             out->data.F32[i] = out->data.C32[i];
    90         }
    91         out->type.type = PS_TYPE_F32;
    92         out->data.U8 = psRealloc(out->data.U8,PSELEMTYPE_SIZEOF(PS_TYPE_F32)*out->nalloc);
    93     }
    94 
    95     return out;
     181    return true;
    96182}
    97 
    98 psVector* psVectorReal(psVector* out, const psVector* in)
    99 {
    100     psElemType type;
    101     psU32 numElements;
    102 
    103     if (in == NULL) {
    104         psFree(out);
    105         return NULL;
    106     }
    107 
    108     type = in->type.type;
    109     numElements = in->n;
    110 
    111     /* if not a complex number, this is logically just a copy */
    112     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    113         // Warn user, as this is probably not expected
    114         psLogMsg(__func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
    115                  "Just a vector copy was performed.");
    116         out = psVectorRecycle(out, numElements, type);
    117         out->n = numElements;
    118         memcpy(out->data.U8, in->data.U8, numElements * PSELEMTYPE_SIZEOF(type));
    119         return out;
    120     }
    121 
    122     if (type == PS_TYPE_C32) {
    123         psF32* outVec;
    124         psC32* inVec = in->data.C32;
    125 
    126         out = psVectorRecycle(out, numElements, PS_TYPE_F32);
    127         out->n = numElements;
    128         outVec = out->data.F32;
    129 
    130         for (psU32 i = 0; i < numElements; i++) {
    131             outVec[i] = crealf(inVec[i]);
    132         }
    133     } else if (type == PS_TYPE_C64) {
    134         psF64* outVec;
    135         psC64* inVec = in->data.C64;
    136 
    137         out = psVectorRecycle(out, numElements, PS_TYPE_F64);
    138         out->n = numElements;
    139         outVec = out->data.F64;
    140 
    141         for (psU32 i = 0; i < numElements; i++) {
    142             outVec[i] = creal(inVec[i]);
    143         }
    144     } else {
    145         char* typeStr;
    146         PS_TYPE_NAME(typeStr,type);
    147         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    148                 _("Specified psVector type, %s, is not supported."),
    149                 typeStr);
    150         psFree(out);
    151         return NULL;
    152     }
    153 
    154     return out;
    155 }
    156 
    157 psVector* psVectorImaginary(psVector* out, const psVector* in)
    158 {
    159     psElemType type;
    160     psU32 numElements;
    161 
    162     if (in == NULL) {
    163         psFree(out);
    164         return NULL;
    165     }
    166 
    167     type = in->type.type;
    168     numElements = in->n;
    169 
    170     /* if not a complex number, this is logically just zeroed image of same size */
    171     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    172         // Warn user, as this is probably not expected
    173         psLogMsg(__func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
    174                  "A zeroed vector was returned.");
    175         out = psVectorRecycle(out, numElements, type);
    176         out->n = numElements;
    177         memset(out->data.U8, 0, PSELEMTYPE_SIZEOF(type) * numElements);
    178         return out;
    179     }
    180 
    181     if (type == PS_TYPE_C32) {
    182         psF32* outVec;
    183         psC32* inVec = in->data.C32;
    184 
    185         out = psVectorRecycle(out, numElements, PS_TYPE_F32);
    186         out->n = numElements;
    187         outVec = out->data.F32;
    188 
    189         for (psU32 i = 0; i < numElements; i++) {
    190             outVec[i] = cimagf(inVec[i]);
    191         }
    192     } else if (type == PS_TYPE_C64) {
    193         psF64* outVec;
    194         psC64* inVec = in->data.C64;
    195 
    196         out = psVectorRecycle(out, numElements, PS_TYPE_F64);
    197         out->n = numElements;
    198         outVec = out->data.F64;
    199 
    200         for (psU32 i = 0; i < numElements; i++) {
    201             outVec[i] = cimag(inVec[i]);
    202         }
    203     } else {
    204         char* typeStr;
    205         PS_TYPE_NAME(typeStr,type);
    206         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    207                 _("Specified psVector type, %s, is not supported."),
    208                 typeStr);
    209         psFree(out);
    210         return NULL;
    211     }
    212 
    213     return out;
    214 }
    215 
    216 psVector* psVectorComplex(psVector* out, const psVector* real, const psVector* imag)
    217 {
    218     psElemType type;
    219     psU32 numElements;
    220 
    221     if (real == NULL || imag == NULL) {
    222         psFree(out);
    223         return NULL;
    224     }
    225 
    226     type = real->type.type;
    227     if (real->n < imag->n) {
    228         numElements = real->n;
    229     } else {
    230         numElements = imag->n;
    231     }
    232 
    233     if (imag->type.type != type) {
    234         char* typeStrReal;
    235         char* typeStrImag;
    236         PS_TYPE_NAME(typeStrReal,type);
    237         PS_TYPE_NAME(typeStrImag,imag->type.type);
    238         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    239                 _("Real psVector type, %s, and imaginary psVector type, %s, must be the same."),
    240                 typeStrReal,typeStrImag);
    241         psFree(out);
    242         return NULL;
    243     }
    244 
    245     if (type == PS_TYPE_F32) {
    246         psC32* outVec;
    247         psF32* realVec = real->data.F32;
    248         psF32* imagVec = imag->data.F32;
    249 
    250         out = psVectorRecycle(out, numElements, PS_TYPE_C32);
    251         out->n = numElements;
    252         outVec = out->data.C32;
    253 
    254         for (psU32 i = 0; i < numElements; i++) {
    255             outVec[i] = realVec[i] + I * imagVec[i];
    256         }
    257     } else if (type == PS_TYPE_F64) {
    258         psC64* outVec;
    259         psF64* realVec = real->data.F64;
    260         psF64* imagVec = imag->data.F64;
    261 
    262         out = psVectorRecycle(out, numElements, PS_TYPE_C64);
    263         out->n = numElements;
    264         outVec = out->data.C64;
    265 
    266         for (psU32 i = 0; i < numElements; i++) {
    267             outVec[i] = realVec[i] + I * imagVec[i];
    268         }
    269     } else {
    270         char* typeStr;
    271         PS_TYPE_NAME(typeStr,type);
    272         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    273                 _("Input psVector type, %s, is required to be either psF32 or psF64."),
    274                 typeStr);
    275         psFree(out);
    276         return NULL;
    277     }
    278 
    279     return out;
    280 }
    281 
    282 psVector* psVectorConjugate(psVector* out, const psVector* in)
    283 {
    284     psElemType type;
    285     psU32 numElements;
    286 
    287     if (in == NULL) {
    288         psFree(out);
    289         return NULL;
    290     }
    291 
    292     type = in->type.type;
    293     numElements = in->n;
    294 
    295     /* if not a complex number, this is logically just a image copy */
    296     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    297         // Warn user, as this is probably not expected
    298         psLogMsg(__func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
    299                  "Vector copy was performed instead.");
    300 
    301         out = psVectorRecycle(out, numElements, type);
    302         out->n = numElements;
    303         memcpy(out->data.U8, in->data.U8, PSELEMTYPE_SIZEOF(type) * numElements);
    304         return out;
    305     }
    306 
    307     if (type == PS_TYPE_C32) {
    308         psC32* outVec;
    309         psC32* inVec = in->data.C32;
    310 
    311         out = psVectorRecycle(out, numElements, PS_TYPE_C32);
    312         out->n = numElements;
    313         outVec = out->data.C32;
    314 
    315         for (psU32 i = 0; i < numElements; i++) {
    316             outVec[i] = crealf(inVec[i]) - I * cimagf(inVec[i]);
    317         }
    318     } else if (type == PS_TYPE_C64) {
    319         psC64* outVec;
    320         psC64* inVec = in->data.C64;
    321 
    322         out = psVectorRecycle(out, numElements, PS_TYPE_C64);
    323         out->n = numElements;
    324         outVec = out->data.C64;
    325 
    326         for (psU32 i = 0; i < numElements; i++) {
    327             outVec[i] = creal(inVec[i]) - I * cimag(inVec[i]);
    328         }
    329     } else {
    330         char* typeStr;
    331         PS_TYPE_NAME(typeStr,type);
    332         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    333                 _("Input psVector type, %s, is required to be either psC32 or psC64."),
    334                 typeStr);
    335         psFree(out);
    336         return NULL;
    337     }
    338 
    339     return out;
    340 }
    341 
    342 psVector* psVectorPowerSpectrum(psVector* out, const psVector* in)
    343 {
    344     psElemType type;
    345     psU32 outNumElements;
    346     psU32 inNumElements;
    347     psU32 inHalfNumElements;
    348     psU32 inNumElementsSquared;
    349 
    350     if (in == NULL) {
    351         psFree(out);
    352         return NULL;
    353     }
    354 
    355     type = in->type.type;
    356     inNumElements = in->n;
    357     inNumElementsSquared = inNumElements * inNumElements;
    358     inHalfNumElements = inNumElements / 2;
    359     outNumElements = inHalfNumElements + 1;
    360 
    361     if (type == PS_TYPE_C32) {
    362         psF32* outVec;
    363         psC32* inVec = in->data.C32;
    364         psF32 inAbs1;
    365         psF32 inAbs2;
    366 
    367         out = psVectorRecycle(out, outNumElements, PS_TYPE_F32);
    368         out->n = outNumElements;
    369         outVec = out->data.F32;
    370 
    371         // from ADD: P_0 = |C_0|^2/N^2
    372         inAbs1 = cabsf(inVec[0]);
    373         outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared;
    374 
    375         // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
    376         for (psU32 i = 1; i < inHalfNumElements; i++) {
    377             inAbs1 = cabsf(inVec[i]);
    378             inAbs2 = cabsf(inVec[inNumElements - i]);
    379             outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared;
    380         }
    381 
    382         // from ADD: P_N/2 = |C_N/2|^2/N^2
    383         inAbs1 = cabsf(inVec[inHalfNumElements]);
    384         outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
    385     } else if (type == PS_TYPE_C64) {
    386         psF64* outVec;
    387         psC64* inVec = in->data.C64;
    388         psF64 inAbs1;
    389         psF64 inAbs2;
    390 
    391         out = psVectorRecycle(out, outNumElements, PS_TYPE_F64);
    392         out->n = outNumElements;
    393         outVec = out->data.F64;
    394 
    395         // from ADD: P_0 = |C_0|^2/N^2
    396         inAbs1 = cabs(inVec[0]);
    397         outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared;
    398 
    399         // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
    400         for (psU32 i = 1; i < inHalfNumElements; i++) {
    401             inAbs1 = cabs(inVec[i]);
    402             inAbs2 = cabs(inVec[inNumElements - i]);
    403             outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared;
    404         }
    405 
    406         // from ADD: P_N/2 = |C_N/2|^2/N^2
    407         inAbs1 = cabs(inVec[inHalfNumElements]);
    408         outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
    409     } else {
    410         char* typeStr;
    411         PS_TYPE_NAME(typeStr,type);
    412         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    413                 _("Input psVector type, %s, is required to be either psC32 or psC64."),
    414                 typeStr);
    415         psFree(out);
    416         return NULL;
    417     }
    418 
    419     return out;
    420 
    421 }

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

    r11248 r11703  
    1 /* @file  psVectorFFT.h
    2  * @brief Contains FFT transform related functions for psVector
    3  *
    4  * @author Robert DeSonia, MHPCC
    5  *
    6  * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2007-01-23 22:47:22 $
    8  * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    9  */
     1/// @file  psVectorFFT.h
     2/// @brief Contains FFT transform related functions for psVector
     3///
     4/// @author Paul Price, IfA
     5/// @author Robert DeSonia, MHPCC
     6///
     7/// @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
     8/// @date $Date: 2007-02-08 04:17:58 $
     9/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     10///
    1011
    1112#ifndef PS_VECTOR_FFT_H
    … …  
    1718/// @{
    1819
    19 /** Specify direction of FFT */
    20 typedef enum {
    21     /// psImageFFT/psVectorFFT should perform a forward FFT.
    22     PS_FFT_FORWARD = 1,
     20/// Forward FFT of a vector
     21///
     22/// Applies a forward FFT (exponent -1), with the result returned in both real and imaginary parts.  The FFT
     23/// is not normalised (a forward followed by a reverse is the original scaled by the number of elements).  The
     24/// FFT takes advantage of the fact that the input is purely real; hence the output size is N/2 + 1 (with
     25/// division rounding down).  Only implemented for F32 input.
     26bool psVectorForwardFFT(psVector **real,///< Real part of FFT
     27                        psVector **imag,///< Imaginary part of FFT
     28                        const psVector *in ///< Input vector (F32)
     29    );
    2330
    24     /// psImageFFT/psVectorFFT should perform a reverse FFT.
    25     PS_FFT_REVERSE = 2,
     31/// Backward FFT of a vector
     32///
     33/// Applies a backward FFT (exponent +1) from the real and imaginary parts with the (purely) real result
     34/// returned.  The FFT is not normalised (a forward followed by a reverse is the original scaled by the number
     35/// of elements).  The FFT takes advantage of the fact that the output will be purely real; hence the input
     36/// size is N/2 + 1 (with division rounding down); to manage the redundancy (is the original size even or
     37/// odd?), we need the original size (the size of the array that was input to psVectorForwardFFT) to be
     38/// provided.  Only implemented for F32 input.
     39bool psVectorBackwardFFT(psVector **out,///< Output vector
     40                         const psVector *real, ///< Real input (F32)
     41                         const psVector *imag, ///< Imaginary input (F32)
     42                         int origN      ///< Original number of elements
     43    );
    2644
    27     /// psImageFFT/psVectorFFT should perform a reverse FFT with a real result.
    28     PS_FFT_REAL_RESULT = 4
    29 } psFFTFlags;
     45/// Power spectrum of a vector
     46///
     47/// Generates the power spectrum of a vector.  Only implemented for F32 input.
     48psVector *psVectorPowerSpectrum(psVector *out, ///< Output power spectrum, or NULL
     49                                const psVector* in ///< Input vector (F32)
     50    );
    3051
    31 
    32 /** Forward and reverse FFT calculations.
    33  *
    34  *  This takes as input the vector of interest (in) and the direction
    35  *  (direction), which is specified by an enumerated type psFftDirection.
    36  *  The input vector may be of type psF32 or psC32, the result is always
    37  *  psC32. If the input vector is psF32, the direction must be forward.
    38  *
    39  *  @return psVector* the FFT transformation result
    40  */
    41 psVector* psVectorFFT(
    42     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    43     const psVector* in,                ///< the vector to apply transform to
    44     psFFTFlags direction               ///< the direction of the transform
    45 );
    46 
    47 /** extract the real portion of a complex vector
    48  *
    49  *  @return psVector*   real portion of the input vector.
    50  */
    51 psVector* psVectorReal(
    52     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    53     const psVector* in                ///< the psVector to extract real portion from
    54 );
    55 
    56 /** extract the imaginary portion of a complex vector
    57  *
    58  *  @return psVector*   imaginary portion of the input vector.
    59  */
    60 psVector* psVectorImaginary(
    61     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    62     const psVector* in                 ///< the psVector to extract imaginary portion from
    63 );
    64 
    65 /** creates a complex vector from separate real and imaginary vectors
    66  *
    67  *  @return psVector*   resulting complex vector
    68  */
    69 psVector* psVectorComplex(
    70     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    71     const psVector* real,              ///< the real vector
    72     const psVector* imag               ///< the imaginary vector
    73 );
    74 
    75 /** computes the complex conjugate of a vector
    76  *
    77  *  @return psVector*   the complex conjugate of the 'in' vector
    78  */
    79 psVector* psVectorConjugate(
    80     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    81     const psVector* in                 ///< the psVector to compute conjugate of
    82 );
    83 
    84 /** computes the power spectrum of a vector
    85  *
    86  *  @return psVector*   the power spectrum of the 'in' vector
    87  */
    88 psVector* psVectorPowerSpectrum(
    89     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    90     const psVector* in                 ///< the psVector to power spectrum of
    91 );
     52/// Multiply complex vectors
     53///
     54/// The input vectors are the real and imaginary parts of each of two vectors.  The real and imaginary parts
     55/// of the output vector are returned.  Only implemented for F32 input.
     56bool psVectorComplexMultiply(psVector **outReal, ///< Real part of output
     57                             psVector **outImag, ///< Imaginary part of output
     58                             const psVector *in1Real, ///< Real part of input 1
     59                             const psVector *in1Imag, ///< Imaginary part of input 1
     60                             const psVector *in2Real, ///< Real part of input 2
     61                             const psVector *in2Imag ///< Imaginary part of input 2
     62    );
    9263
    9364/// @}

  • trunk/psLib/src/imageops/psImageConvolve.c

    r11153 r11703  
    1 /*  @file  psImageConvolve.c
    2  *
    3  *  @brief Contains FFT transform related functions for psImage.
    4  *
    5  *  @author Robert DeSonia, MHPCC
    6  *
    7  *  @version $Revision: 1.42 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-01-19 04:30:33 $
    9  *
    10  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    11  */
     1/// @file  psImageConvolve.c
     2///
     3/// @brief Contains FFT transform related functions for psImage.
     4///
     5/// @author Robert DeSonia, MHPCC
     6/// @author Paul Price, IfA
     7/// @author Eugene Magnier, IfA
     8///
     9/// @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
     10/// @date $Date: 2007-02-08 04:17:58 $
     11///
     12/// Copyright 2004-2007 Institute for Astronomy, University of Hawaii
     13///
    1214
    1315#ifdef HAVE_CONFIG_H
    14 # include "config.h"
     16#include "config.h"
    1517#endif
    1618
    1719#include <string.h>
    1820#include <math.h>
    19 #include "psImageConvolve.h"
    20 #include "psImageFFT.h"
    21 #include "psImageStructManip.h"
    22 #include "psBinaryOp.h"
     21#include "psAbort.h"
    2322#include "psMemory.h"
    2423#include "psLogMsg.h"
    2524#include "psError.h"
    2625#include "psAssert.h"
    27 
    28 
    29 
    30 #define FOURIER_PADDING 32 /* padding amount in every side of the image for fourier convolution */
    31 
    32 static void freeKernel(psKernel* ptr);
    33 
    34 psKernel* psKernelAlloc(int xMin, int xMax, int yMin, int yMax)
    35 {
    36     psKernel* result;
    37     psS32 numRows;
    38     psS32 numCols;
    39 
    40     // following is explicitly spelled out in the SDRS as a requirement
     26#include "psScalar.h"
     27#include "psBinaryOp.h"
     28#include "psImageFFT.h"
     29#include "psImageStructManip.h"
     30#include "psImagePixelManip.h"
     31
     32#include "psImageConvolve.h"
     33
     34static void kernelFree(psKernel *kernel)
     35{
     36    if (kernel) {
     37        psFree(kernel->image);
     38        psFree(kernel->p_kernelRows);
     39    }
     40    return;
     41}
     42
     43psKernel *psKernelAlloc(int xMin, int xMax, int yMin, int yMax)
     44{
     45    // Check the inputs to make sure max > min; if not, switch.
    4146    if (yMin > yMax) {
    4247        psLogMsg(__func__, PS_LOG_WARN,
    … …  
    4449                 yMin, yMax);
    4550
    46         psS32 temp = yMin;
     51        int temp = yMin;
    4752        yMin = yMax;
    4853        yMax = temp;
    … …  
    5560                 xMin, xMax);
    5661
    57         psS32 temp = xMin;
     62        int temp = xMin;
    5863        xMin = xMax;
    5964        xMax = temp;
    6065    }
    6166
    62     numRows = yMax - yMin + 1;
    63     numCols = xMax - xMin + 1;
    64 
    65     result = psAlloc(sizeof(psKernel));
    66     result->xMin = xMin;
    67     result->xMax = xMax;
    68     result->yMin = yMin;
    69     result->yMax = yMax;
    70     result->image = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL);
    71     memset((result->image->p_rawDataBuffer),0,numCols*numRows*PSELEMTYPE_SIZEOF(PS_TYPE_KERNEL));
    72     result->p_kernelRows = psAlloc(sizeof(float*)*numRows);
    73 
    74     float** kernelRows = result->p_kernelRows;
    75     float** imageRows = result->image->data.PS_TYPE_KERNEL_DATA;
    76     for (psS32 i = 0; i < numRows; i++) {
    77         kernelRows[i] = imageRows[i] - xMin;
    78     }
    79     result->kernel = kernelRows - yMin;
    80 
    81     psMemSetDeallocator(result,(psFreeFunc)freeKernel);
    82 
    83     return result;
    84 }
    85 
    86 void freeKernel(psKernel* ptr)
    87 {
    88     if (ptr != NULL) {
    89         psFree(ptr->image);
    90         psFree(ptr->p_kernelRows);
    91     }
    92 }
     67    int numRows = yMax - yMin + 1;      // Number of rows for kernel image
     68    int numCols = xMax - xMin + 1;      // Number of columns for kernel image
     69
     70    psKernel *kernel = psAlloc(sizeof(psKernel)); // The kernel, to be returned
     71    psMemSetDeallocator(kernel,(psFreeFunc)kernelFree);
     72
     73    kernel->xMin = xMin;
     74    kernel->xMax = xMax;
     75    kernel->yMin = yMin;
     76    kernel->yMax = yMax;
     77    kernel->image = psImageAlloc(numCols, numRows, PS_TYPE_KERNEL);
     78    psImageInit(kernel->image, 0.0);
     79
     80    // Set up indirections, so we can refer to kernel->kernel[-1][-3] for the (-1,-1) element, instead of
     81    // kernel->image[kernel->yMax - kernel->yMin + 1][kernel->yMax - kernel->yMin + 3] (yuk!).
     82    kernel->p_kernelRows = psAlloc(sizeof(float*)*numRows);
     83    for (int i = 0; i < numRows; i++) {
     84        kernel->p_kernelRows[i] = kernel->image->data.PS_TYPE_KERNEL_DATA[i] - xMin;
     85    }
     86    kernel->kernel = kernel->p_kernelRows - yMin;
     87
     88    return kernel;
     89}
     90
    9391
    9492
    … …  
    9694{
    9795    PS_ASSERT_PTR(ptr, false);
    98     return ( psMemGetDeallocator(ptr) == (psFreeFunc)freeKernel );
    99 }
    100 
    101 
    102 psKernel* psKernelGenerate(const psVector* tShifts,
    103                            const psVector* xShifts,
    104                            const psVector* yShifts,
    105                            bool relative)
    106 {
    107     psS32 lastX;
    108     psS32 lastY;
    109     psS32 lastT;
    110     psS32 x;
    111     psS32 y;
    112     psS32 t;
    113     psS32 xMin = 0;
    114     psS32 xMax = 0;
    115     psS32 yMin = 0;
    116     psS32 yMax = 0;
    117     psS32 length = 0;
    118     float normalizeTime = 1.0;  // fraction of total time for each shift clock
    119     psKernel* result = NULL;
    120     float** kernel = NULL;
    121 
    122     // got non-NULL vectors?
    123     if (tShifts == NULL || xShifts == NULL || yShifts == NULL) {
    124         psError(PS_ERR_BAD_PARAMETER_NULL, true,
    125                 _("Specified shift vectors can not be NULL."));
     96    return ( psMemGetDeallocator(ptr) == (psFreeFunc)kernelFree );
     97}
     98
     99
     100psKernel *psKernelGenerate(const psVector *tShifts,
     101                           const psVector *xShifts,
     102                           const psVector *yShifts,
     103                           bool tRelative,
     104                           bool xyRelative)
     105{
     106    PS_ASSERT_VECTOR_NON_NULL(tShifts, NULL);
     107    PS_ASSERT_VECTOR_NON_NULL(xShifts, NULL);
     108    PS_ASSERT_VECTOR_NON_NULL(yShifts, NULL);
     109    PS_ASSERT_VECTORS_SIZE_EQUAL(tShifts, xShifts, NULL);
     110    PS_ASSERT_VECTORS_SIZE_EQUAL(tShifts, yShifts, NULL);
     111    PS_ASSERT_VECTOR_TYPE(tShifts, PS_TYPE_F32, NULL);
     112    PS_ASSERT_VECTOR_TYPE(xShifts, PS_TYPE_S32, NULL);
     113    PS_ASSERT_VECTOR_TYPE(yShifts, PS_TYPE_S32, NULL);
     114
     115    // If there are no shifts, the kernel is just a 1 at 0,0
     116    long num = tShifts->n;              // Number of shifts
     117    if (num == 0) {
     118        psKernel *kernel = psKernelAlloc(0,0,0,0);
     119        kernel->kernel[0][0] = 1;
     120        return kernel;
     121    }
     122
     123    // Get dimensions and scaling
     124    int xMin, xMax, yMin, yMax;         // Range of values for kernel
     125    int xLast, yLast;                   // Last location, for relative shifts
     126    xLast = xMin = xMax = xShifts->data.S32[0];
     127    yLast = yMin = yMax = yShifts->data.S32[0];
     128    float tSum = tShifts->data.F32[0];   // Sum of the times
     129    for (long i = 1; i < num; i++) {
     130        int x = xShifts->data.S32[i];    // x position in kernel
     131        int y = yShifts->data.S32[i];    // y position in kernel
     132        if (xyRelative) {
     133            x += xLast;
     134            y += yLast;
     135            xLast = x;
     136            yLast = y;
     137        }
     138        if (x < xMin) {
     139            xMin = x;
     140        }
     141        if (x > xMax) {
     142            xMax = x;
     143        }
     144        if (y < yMin) {
     145            yMin = y;
     146        }
     147        if (y > yMax) {
     148            yMax = y;
     149        }
     150
     151        if (tRelative) {
     152            tSum += tShifts->data.F32[i];
     153        }
     154    }
     155
     156    if (!tRelative) {
     157        // Then the total time is simply the final value
     158        // NB: We assume the counter starts at zero!
     159        tSum = tShifts->data.F32[tShifts->n - 1];
     160    }
     161
     162    // One more pass through to set the kernel
     163    psKernel *kernel = psKernelAlloc(xMin, xMax, yMin, yMax); // The kernel
     164    xLast = xShifts->data.S32[0];
     165    yLast = yShifts->data.S32[0];
     166    float tLast = 0.0;                  // Last value for t
     167    for (int i = 0; i < num; i++) {
     168        int x = xShifts->data.S32[i];    // x position in kernel
     169        int y = yShifts->data.S32[i];    // y position in kernel
     170        if (xyRelative) {
     171            x += xLast;
     172            y += yLast;
     173            xLast = x;
     174            yLast = y;
     175        }
     176        float t = tShifts->data.F32[i];
     177        if (tRelative) {
     178            t -= tLast;
     179            tLast = tShifts->data.F32[i];
     180        }
     181
     182        kernel->kernel[y][x] += t;
     183    }
     184
     185    // Normalise the kernel by the total time (kernel sum should be unity)
     186    psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(1.0 / tSum, PS_TYPE_F32));
     187
     188    return kernel;
     189}
     190
     191psImage *psImageConvolveDirect(psImage *out,
     192                               const psImage *in,
     193                               const psKernel *kernel)
     194{
     195    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
     196    PS_ASSERT_IMAGE_TYPE_F32_OR_F64(in, NULL);
     197    PS_ASSERT_PTR_NON_NULL(kernel, NULL);
     198    PS_ASSERT_PTR_NON_NULL(kernel->kernel, NULL);
     199
     200    // Pull out kernel parameters, for convenience
     201    int xMin = kernel->xMin;
     202    int xMax = kernel->xMax;
     203    int yMin = kernel->yMin;
     204    int yMax = kernel->yMax;
     205    float **kernelData = kernel->kernel;
     206
     207    int numRows = in->numRows;          // Number of rows
     208    int numCols = in->numCols;          // Number of columns
     209
     210#define SPATIAL_CONVOLVE_CASE(TYPE) \
     211    case PS_TYPE_##TYPE: { \
     212        ps##TYPE **inData = in->data.TYPE; \
     213        out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \
     214        for (int row = 0; row < numRows; row++) { \
     215            ps##TYPE *outRow = out->data.TYPE[row]; \
     216            for (int col = 0; col < numCols; col++) { \
     217                ps##TYPE pixel = 0.0; \
     218                for (int kRow = PS_MAX(yMin, -row); kRow <= PS_MIN(yMax, numRows - row - 1); kRow++) { \
     219                    for (int kCol = PS_MAX(xMin, -col); kCol <= PS_MIN(xMax, numCols - col - 1); kCol++) { \
     220                        pixel += kernelData[kRow][kCol] * inData[row + kRow][col + kCol]; \
     221                    } \
     222                } \
     223                outRow[col] = pixel; \
     224            } \
     225        } \
     226    } \
     227    break;
     228
     229    switch (in->type.type) {
     230        SPATIAL_CONVOLVE_CASE(F32);
     231        SPATIAL_CONVOLVE_CASE(F64);
     232      default:
     233        psAbort(PS_FILE_LINE, "Should never get here: bad type that was asserted on previously.");
     234    }
     235
     236    return out;
     237}
     238
     239psImage *psImageConvolveFFT(psImage *out,
     240                            const psImage *in,
     241                            const psKernel *kernel,
     242                            float pad)
     243{
     244    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
     245    PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, NULL);
     246    PS_ASSERT_PTR_NON_NULL(kernel, NULL);
     247    PS_ASSERT_IMAGE_NON_NULL(kernel->image, NULL);
     248
     249    // Pull out kernel parameters, for convenience
     250    int xMin = kernel->xMin;
     251    int xMax = kernel->xMax;
     252    int yMin = kernel->yMin;
     253    int yMax = kernel->yMax;
     254
     255    int numRows = in->numRows;          // Number of rows in input image
     256    int numCols = in->numCols;          // Number of columns in input image
     257
     258    // Need to pad the input image to protect from wrap-around effects
     259    if (xMax - xMin > numCols || yMax - yMin > numRows) {
     260        // Cannot pad the image if the kernel is larger.
     261        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     262                _("Kernel cannot extend further than input image size (%dx%d vs %dx%d)."),
     263                xMax, yMax, numCols, numRows);
    126264        return NULL;
    127265    }
    128 
    129     // types match?
    130     if (xShifts->type.type != yShifts->type.type ||
    131             tShifts->type.type != xShifts->type.type) {
    132         char* typeXStr;
    133         char* typeYStr;
    134         char* typeTStr;
    135         PS_TYPE_NAME(typeXStr,xShifts->type.type);
    136         PS_TYPE_NAME(typeYStr,yShifts->type.type);
    137         PS_TYPE_NAME(typeTStr,tShifts->type.type);
    138         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    139                 _("Input t-, x-, and y-shift vector types (%s/%s/%s) must match."),
    140                 typeTStr, typeXStr, typeYStr);
     266    int paddedCols = numCols + PS_MAX(-xMin, xMax); // Number of columns in padded image
     267    int paddedRows = numRows + PS_MAX(-yMin, yMax); // Number of rows in padded image
     268
     269    // Generate padded image
     270    psImage *paddedImage = psImageAlloc(paddedCols,paddedRows,in->type.type); // Padded input image
     271    psImageOverlaySection(paddedImage, in, 0, 0, "=");
     272    for (int y = 0; y < numRows; y++) {
     273        for (int x = numCols; x < paddedCols; x++) {
     274            paddedImage->data.F32[y][x] = pad;
     275        }
     276    }
     277    for (int y = numRows; y < paddedRows; y++) {
     278        for (int x = 0; x < paddedCols; x++) {
     279            paddedImage->data.F32[y][x] = pad;
     280        }
     281    }
     282
     283    // Result of FFT
     284    psImage *inRealFFT = NULL, *inImagFFT = NULL;
     285    if (!psImageForwardFFT(&inRealFFT, &inImagFFT, paddedImage)) {
     286        psError(PS_ERR_UNKNOWN, false, _("Failed to fourier transform input image."));
     287        psFree(paddedImage);
    141288        return NULL;
    142289    }
    143 
    144     // sizes match?
    145     length = xShifts->n;
    146     if (length != yShifts->n ||
    147             length != tShifts->n) {
    148         psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    149                 "Shift vectors can not be of different sizes.");
     290    psFree(paddedImage);
     291
     292    // Generate padded kernel image
     293    psImage *paddedKernel = psImageAlloc(paddedCols, paddedRows, PS_TYPE_F32);
     294    psImageInit(paddedKernel, 0.0);
     295    for (int y = PS_MIN(-1, yMin); y <= PS_MIN(-1, yMax); y++) {
     296        // y is negative
     297        for (int x = PS_MIN(-1, xMin); x <= PS_MIN(-1, xMax); x++) {
     298            // x is negative
     299            paddedKernel->data.F32[paddedRows + y][paddedCols + x] = kernel->kernel[y][x];
     300        }
     301        for (int x = PS_MAX(0, xMin); x <= PS_MAX(0, xMax); x++) {
     302            // x is positive
     303            paddedKernel->data.F32[paddedRows + y][x] = kernel->kernel[y][x];
     304        }
     305    }
     306    for (int y = PS_MAX(0, yMin); y <= PS_MAX(0, yMax); y++) {
     307        // y is positive
     308        for (int x = PS_MIN(-1, xMin); x <= PS_MIN(-1, xMax); x++) {
     309            // x is negative
     310            paddedKernel->data.F32[y][paddedCols + x] = kernel->kernel[y][x];
     311        }
     312        for (int x = PS_MAX(0, xMin); x <= PS_MAX(0, xMax); x++) {
     313            // x is positive
     314            paddedKernel->data.F32[y][x] = kernel->kernel[y][x];
     315        }
     316    }
     317
     318    psImage *kernelRealFFT = NULL, *kernelImagFFT = NULL;
     319    if (!psImageForwardFFT(&kernelRealFFT, &kernelImagFFT, paddedKernel)) {
     320        psError(PS_ERR_UNKNOWN, false, _("Failed to fourier transform kernel."));
     321        psFree(inRealFFT);
     322        psFree(inImagFFT);
     323        psFree(paddedKernel);
    150324        return NULL;
    151325    }
    152 
    153     // if no shifts, the kernel is just a 1 at 0,0
    154     if (length < 1) {
    155         result = psKernelAlloc(0,0,0,0);
    156         result->kernel[0][0] = 1;
    157         return result;
    158     }
    159 
    160     #define KERNEL_GENERATE_CASE(TYPE) \
    161 case PS_TYPE_##TYPE: { \
    162         ps##TYPE *tShiftData = tShifts->data.TYPE; \
    163         ps##TYPE *xShiftData = xShifts->data.TYPE; \
    164         ps##TYPE *yShiftData = yShifts->data.TYPE; \
    165         lastX = xShiftData[length-1]; \
    166         lastY = yShiftData[length-1]; \
    167         lastT = tShiftData[length-1]; \
    168         \
    169         for (int lcv = 0; lcv < length; lcv++) { \
    170             x = lastX - xShiftData[lcv]; \
    171             y = lastY - yShiftData[lcv]; \
    172             \
    173             if (x < xMin) { \
    174                 xMin = x; \
    175             } else if (x > xMax) { \
    176                 xMax = x; \
    177             } \
    178             if (y < yMin) { \
    179                 yMin = y; \
    180             } else if (y > yMax) { \
    181                 yMax = y; \
    182             } \
    183         } \
    184         \
    185         normalizeTime = 1.0 / (float)(tShiftData[length-1]); \
    186         result = psKernelAlloc(xMin,xMax,yMin,yMax); \
    187         kernel = result->kernel; \
    188         \
    189         psS32 prevT = 0; \
    190         for (int i = 0; i < length; i++) { \
    191             t = tShiftData[i] - prevT; \
    192             x = lastX - xShiftData[i]; \
    193             y = lastY - yShiftData[i]; \
    194             \
    195             kernel[y][x] += (float)t / (float)lastT; \
    196             prevT = tShiftData[i]; \
    197         } \
    198         break; \
    199     }
    200 
    201     #define RELATIVE_KERNEL_GENERATE_CASE(TYPE) \
    202 case PS_TYPE_##TYPE: { \
    203         ps##TYPE *tShiftData = tShifts->data.TYPE; \
    204         ps##TYPE *xShiftData = xShifts->data.TYPE; \
    205         ps##TYPE *yShiftData = yShifts->data.TYPE; \
    206         \
    207         x = 0; \
    208         y = 0; \
    209         t = 0; \
    210         \
    211         for (int lcv = length-1; lcv >= 0; lcv--) { \
    212             t += tShiftData[lcv]; \
    213             \
    214             if (x < xMin) { \
    215                 xMin = x; \
    216             } else if (x > xMax) { \
    217                 xMax = x; \
    218             } \
    219             if (y < yMin) { \
    220                 yMin = y; \
    221             } else if (y > yMax) { \
    222                 yMax = y; \
    223             } \
    224             x -= xShiftData[lcv]; \
    225             y -= yShiftData[lcv]; \
    226             \
    227         } \
    228         result = psKernelAlloc(xMin,xMax,yMin,yMax); \
    229         kernel = result->kernel; \
    230         \
    231         normalizeTime = 1.0 / (float)t; \
    232         x = 0; \
    233         y = 0; \
    234         for (psS32 i = length-1; i >= 0; i--) { \
    235             kernel[y][x] += (float)(tShiftData[i]) * normalizeTime; \
    236             x -= xShiftData[i]; \
    237             y -= yShiftData[i]; \
    238             \
    239         } \
    240         break; \
    241     }
    242 
    243     if (relative) {
    244         switch (xShifts->type.type) {
    245             RELATIVE_KERNEL_GENERATE_CASE(U8);
    246             RELATIVE_KERNEL_GENERATE_CASE(U16);
    247             RELATIVE_KERNEL_GENERATE_CASE(U32);
    248             RELATIVE_KERNEL_GENERATE_CASE(U64);
    249             RELATIVE_KERNEL_GENERATE_CASE(S8);
    250             RELATIVE_KERNEL_GENERATE_CASE(S16);
    251             RELATIVE_KERNEL_GENERATE_CASE(S32);
    252             RELATIVE_KERNEL_GENERATE_CASE(S64);
    253             RELATIVE_KERNEL_GENERATE_CASE(F32);
    254             RELATIVE_KERNEL_GENERATE_CASE(F64);
    255             RELATIVE_KERNEL_GENERATE_CASE(C32);
    256             RELATIVE_KERNEL_GENERATE_CASE(C64);
    257 
    258         default: {
    259                 char* typeStr;
    260                 PS_TYPE_NAME(typeStr,xShifts->type.type);
    261                 psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    262                         _("Specified psImage type, %s, is not supported."),
    263                         typeStr);
    264             }
    265         }
    266     } else {
    267         switch (xShifts->type.type) {
    268             KERNEL_GENERATE_CASE(U8);
    269             KERNEL_GENERATE_CASE(U16);
    270             KERNEL_GENERATE_CASE(U32);
    271             KERNEL_GENERATE_CASE(U64);
    272             KERNEL_GENERATE_CASE(S8);
    273             KERNEL_GENERATE_CASE(S16);
    274             KERNEL_GENERATE_CASE(S32);
    275             KERNEL_GENERATE_CASE(S64);
    276             KERNEL_GENERATE_CASE(F32);
    277             KERNEL_GENERATE_CASE(F64);
    278             KERNEL_GENERATE_CASE(C32);
    279             KERNEL_GENERATE_CASE(C64);
    280 
    281         default: {
    282                 char* typeStr;
    283                 PS_TYPE_NAME(typeStr,xShifts->type.type);
    284                 psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    285                         _("Specified psImage type, %s, is not supported."),
    286                         typeStr);
    287             }
    288         }
    289     }
    290 
    291     return result;
    292 }
    293 
    294 psImage* psImageConvolve(psImage* out,
    295                          const psImage* in,
    296                          const psKernel* kernel,
    297                          bool direct)
    298 {
    299     if (in == NULL) {
    300         psFree(out);
     326    psFree(paddedKernel);
     327
     328    // Convolution in fourier domain is just a pixel-wise multiplication
     329    if (!psImageComplexMultiply(&inRealFFT, &inImagFFT, inRealFFT, inImagFFT, kernelRealFFT, kernelImagFFT)) {
     330        psError(PS_ERR_UNKNOWN, false, _("Unable to multiply fourier transformts."));
     331        psFree(inRealFFT);
     332        psFree(inImagFFT);
     333        psFree(kernelRealFFT);
     334        psFree(kernelImagFFT);
    301335        return NULL;
    302336    }
    303 
    304     if (kernel == NULL) {
    305         psError(PS_ERR_BAD_PARAMETER_NULL, true,
    306                 _("Specified psKernel can not be NULL."));
    307         psFree(out);
     337    psFree(kernelRealFFT);
     338    psFree(kernelImagFFT);
     339
     340    psImage *paddedConvolved = NULL; // Padded convolved image
     341    if (!psImageBackwardFFT(&paddedConvolved, inRealFFT, inImagFFT, paddedCols)) {
     342        psError(PS_ERR_UNKNOWN, false, _("Failed to invert fourier transform of convolution image."));
     343        psFree(inRealFFT);
     344        psFree(inImagFFT);
    308345        return NULL;
    309346    }
    310     psS32 xMin = kernel->xMin;
    311     psS32 xMax = kernel->xMax;
    312     psS32 yMin = kernel->yMin;
    313     psS32 yMax = kernel->yMax;
    314     float** kData = kernel->kernel;
    315 
    316     // make the output image to the proper size and type
    317     psS32 numRows = in->numRows;
    318     psS32 numCols = in->numCols;
    319 
    320 
    321 
    322     if (direct) {
    323         // spatial convolution
    324 
    325         #define SPATIAL_CONVOLVE_CASE(TYPE) \
    326     case PS_TYPE_##TYPE: { \
    327             ps##TYPE** inData = in->data.TYPE; \
    328             out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \
    329             for (psS32 row=0;row<numRows;row++) { \
    330                 ps##TYPE* outRow = out->data.TYPE[row]; \
    331                 for (psS32 col=0;col<numCols;col++) { \
    332                     ps##TYPE pixel = 0.0; \
    333                     for (psS32 kRow = yMin; kRow < yMax; kRow++) { \
    334                         if (row-kRow >= 0 && row-kRow < numRows) { \
    335                             for (psS32 kCol = xMin; kCol < xMax; kCol++) { \
    336                                 if (col-kCol >= 0 && col-kCol < numCols) { \
    337                                     pixel += kData[kRow][kCol] * inData[row-kRow][col-kCol]; \
    338                                 } \
    339                             } \
    340                         } \
    341                     } \
    342                     outRow[col] = pixel; \
    343                 } \
    344             } \
    345         } \
    346         break;
    347 
    348         switch (in->type.type) {
    349             SPATIAL_CONVOLVE_CASE(U8)
    350             SPATIAL_CONVOLVE_CASE(U16)
    351             SPATIAL_CONVOLVE_CASE(U32)
    352             SPATIAL_CONVOLVE_CASE(U64)
    353             SPATIAL_CONVOLVE_CASE(S8)
    354             SPATIAL_CONVOLVE_CASE(S16)
    355             SPATIAL_CONVOLVE_CASE(S32)
    356             SPATIAL_CONVOLVE_CASE(S64)
    357             SPATIAL_CONVOLVE_CASE(F32)
    358             SPATIAL_CONVOLVE_CASE(F64)
    359             SPATIAL_CONVOLVE_CASE(C32)
    360             SPATIAL_CONVOLVE_CASE(C64)
    361 
    362         default: {
    363                 char* typeStr;
    364                 PS_TYPE_NAME(typeStr,in->type.type);
    365                 psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    366                         _("Specified psImage type, %s, is not supported."),
    367                         typeStr);
    368                 psFree(out);
    369                 return NULL;
    370 
    371             }
    372         }
    373 
    374 
    375     } else {
    376         // fourier convolution
    377         psS32 paddedCols = numCols+2*FOURIER_PADDING;
    378         psS32 paddedRows = numRows+2*FOURIER_PADDING;
    379 
    380         // check to see if kernel is smaller, otherwise padding it up will fail.
    381         psS32 kRows = kernel->image->numRows;
    382         psS32 kCols = kernel->image->numCols;
    383         if (kRows >= numRows || kCols >= numCols) {
    384             psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    385                     _("Specified psKernel size, %dx%d, can not be larger than input psImage size, %dx%d."),
    386                     kCols,kRows,
    387                     numCols, numRows);
    388             psFree(out);
    389             return NULL;
    390         }
    391 
    392         // pad the image
    393         psImage* paddedImage = psImageAlloc(paddedCols,paddedRows,in->type.type);
    394         psS32 elementSize = PSELEMTYPE_SIZEOF(in->type.type);
    395 
    396         // zero out padded area on top and bottom
    397         memset(paddedImage->data.U8[0],0,FOURIER_PADDING*paddedCols*elementSize);
    398         memset(paddedImage->data.U8[FOURIER_PADDING+numRows-1],0,FOURIER_PADDING*paddedCols*elementSize);
    399 
    400         // fill in the image-containing rows.
    401         psS32 sidePaddingSize = FOURIER_PADDING*elementSize;
    402         psS32 imageRowSize = numCols*elementSize;
    403         psU8* paddedData = paddedImage->data.U8[FOURIER_PADDING];
    404         for (psS32 row=0;row<numRows;row++) {
    405             // zero out padded area on left edge.
    406             memset(paddedData,0,sidePaddingSize);
    407             paddedData += sidePaddingSize;
    408             memcpy(paddedData,in->data.U8[row],imageRowSize);
    409             paddedData += imageRowSize;
    410             // zero out padded area on right edge.
    411             memset(paddedData,0,sidePaddingSize);
    412             paddedData += sidePaddingSize;
    413         }
    414 
    415         // pad the kernel to the same size of paddedImage
    416         psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL);
    417         memset(paddedKernel->data.U8[0],0,sizeof(float)*numCols*numRows); // zero-out image
    418         psS32 yMax = kernel->yMax;
    419         psS32 xMax = kernel->xMax;
    420         for (psS32 row = kernel->yMin; row <= yMax;row++) {
    421             psS32 padRow = row;
    422             if (padRow < 0) {
    423                 padRow += paddedRows;
    424             }
    425             float* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow];
    426             float* kernelRow = kernel->kernel[row];
    427             for (psS32 col = kernel->xMin; col <= xMax; col++) {
    428                 if (col < 0) {
    429                     padData[col+paddedCols] = kernelRow[col];
    430                 } else {
    431                     padData[col] = kernelRow[col];
    432                 }
    433             }
    434         }
    435 
    436         psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD);
    437         if (kernelFourier == NULL) {
    438             psError(PS_ERR_UNKNOWN, false,
    439                     _("Failed to perform a fourier transform of kernel."));
    440             psFree(out);
    441             return NULL;
    442         }
    443 
    444         psImage* inFourier = psImageFFT(NULL, paddedImage, PS_FFT_FORWARD);
    445         if (inFourier == NULL) {
    446             psError(PS_ERR_UNKNOWN, false,
    447                     _("Failed to perform a fourier transform of input image."));
    448             psFree(out);
    449             return NULL;
    450         }
    451 
    452         // convolution in fourier domain is just a pixel-wise multiplication
    453         for (int row = 0; row < paddedRows; row++) {
    454             psC32* inRow = inFourier->data.C32[row];
    455             psC32* kRow = kernelFourier->data.C32[row];
    456             for (int col = 0; col < paddedCols; col++) {
    457                 inRow[col] *= kRow[col];
    458             }
    459         }
    460 
    461         psImage* complexOut = psImageFFT(NULL, inFourier,
    462                                          PS_FFT_REVERSE);
    463         if (complexOut == NULL) {
    464             psError(PS_ERR_UNKNOWN, false,
    465                     _("Failed to perform a fourier transform of input image."));
    466             psFree(out);
    467             return NULL;
    468         }
    469 
    470         // subset out the padded area now.
    471         psImage* complexOutSansPad = psImageSubset(complexOut,
    472                                      psRegionSet(FOURIER_PADDING, FOURIER_PADDING+numCols,FOURIER_PADDING,FOURIER_PADDING+numRows));
    473 
    474         out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
    475         float factor = 1.0f/(float)paddedCols/(float)paddedRows;
    476         for (psS32 row = 0; row < numRows; row++) {
    477             psF32* outRow = out->data.F32[row];
    478             psC32* resultRow = complexOutSansPad->data.C32[row];
    479             for (psS32 col = 0; col < numCols; col++) {
    480                 outRow[col] = crealf(resultRow[col])*factor;
    481             }
    482         }
    483 
    484         psFree(complexOut); // frees complexOutSansPad, as it is a child.
    485         psFree(kernelFourier);
    486         psFree(inFourier);
    487         psFree(paddedImage);
    488         psFree(paddedKernel);
    489 
    490     }
     347    psFree(inRealFFT);
     348    psFree(inImagFFT);
     349
     350    // Trim off the padding, then renormalise (which also does a copy, so there's no parent for the output)
     351    psImage *convolved = psImageSubset(paddedConvolved, psRegionSet(0, numCols, 0, numRows));
     352    out = (psImage*)psBinaryOp(out, convolved, "*",
     353                               psScalarAlloc(1.0 / paddedCols / paddedRows, PS_TYPE_F32));
     354    psFree(convolved);
     355    psFree(paddedConvolved);
    491356
    492357    return out;
    … …  
    661526        IMAGESMOOTH_CASE(F64);
    662527    default: {
    663             char* typeStr;
     528            char *typeStr;
    664529            PS_TYPE_NAME(typeStr,image->type.type);
    665530            psError(PS_ERR_BAD_PARAMETER_TYPE, true,

  • trunk/psLib/src/imageops/psImageConvolve.h

    r11248 r11703  
    55 * @author Robert DeSonia, MHPCC
    66 *
    7  * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2007-01-23 22:47:23 $
     7 * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2007-02-08 04:17:58 $
    99 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    1010 */
    … …  
    2020#include "psType.h"
    2121
    22 #define PS_TYPE_KERNEL PS_TYPE_F32     /**< the data member to use for kernel image */
    23 #define PS_TYPE_KERNEL_DATA F32        /**< the data member to use for kernel image */
    24 #define PS_TYPE_KERNEL_NAME "psF32"    /**< the data type for kernel as a string */
     22#define PS_TYPE_KERNEL PS_TYPE_F32     ///< the data member to use for kernel image */
     23#define PS_TYPE_KERNEL_DATA F32        ///< the data member to use for kernel image */
     24#define PS_TYPE_KERNEL_NAME "psF32"    ///< the data type for kernel as a string */
    2525
    26 /** Kernel Type
    27  *
    28  *  A floating-point data type used for storing kernel data.
    29  *
    30  */
    31 //typedef float psKernelType;
    32 
    33 /** A convolution kernel */
     26/// A convolution kernel
    3427typedef struct
    3528{
    36     psImage* image;                    ///< Kernel data, in the form of an image
     29    psImage *image;                    ///< Kernel data, in the form of an image
    3730    int xMin;                          ///< Most negative x index
    3831    int yMin;                          ///< Most negative y index
    3932    int xMax;                          ///< Most positive x index
    4033    int yMax;                          ///< Most positive y index
    41     float** kernel;                    ///< Pointer to the kernel data
    42     float** p_kernelRows;              ///< Pointer to the rows of the kernel data; not intended for user use.
     34    float **kernel;                    ///< Pointer to the kernel data
     35    float **p_kernelRows;              ///< Pointer to the rows of the kernel data; not intended for user use.
    4336}
    4437psKernel;
    4538
    46 /** Allocates a convolution kernel of the given range
    47  *
    48  * In order to perform a convolution, we need to define the convolution
    49  * kernel. We need a more general object than a psImage so that we can
    50  * incorporate the offset from the (0, 0) pixel to the (0, 0) value of the
    51  * kernel. It might be convenient to allow both positive and negative
    52  * indices to convey the positive and negative shifts. One might consider
    53  * setting the x0 and y0 members of a psImage to the appropriate offsets,
    54  * but this is not the purpose of these members, and doing so may affect the
    55  * behavior of other psImage operations.
    56  *
    57  * This construction allows the kernel member to use negative indices, while
    58  * preserving the location of psMemBlocks relative to allocated memory.
    59  *
    60  * The maximum extent of the kernel shifts shall be defined by the xMin,
    61  * xMax, yMin and yMax members. Note that xMin and yMin, under normal
    62  * circumstances, should be negative numbers. That is,
    63  * myKernel->kernel[-3][-2] may be defined if yMin and xMin are equal to or
    64  * more negative than -3 and -2, respectively.
    65  *
    66  * In the event that one of the minimum values is greater than the
    67  * corresponding maximum value, the function shall generate a warning, and
    68  * the offending values shall be exchanged.
    69  *
    70  * @return psKernel*          A new kernel object
    71  */
    72 psKernel* psKernelAlloc(
     39/// Allocates a convolution kernel of the given range
     40///
     41/// In order to perform a convolution, we need to define the convolution
     42/// kernel. We need a more general object than a psImage so that we can
     43/// incorporate the offset from the (0, 0) pixel to the (0, 0) value of the
     44/// kernel. It might be convenient to allow both positive and negative
     45/// indices to convey the positive and negative shifts. One might consider
     46/// setting the x0 and y0 members of a psImage to the appropriate offsets,
     47/// but this is not the purpose of these members, and doing so may affect the
     48/// behavior of other psImage operations.
     49///
     50/// This construction allows the kernel member to use negative indices, while
     51/// preserving the location of psMemBlocks relative to allocated memory.
     52///
     53/// The maximum extent of the kernel shifts shall be defined by the xMin,
     54/// xMax, yMin and yMax members. Note that xMin and yMin, under normal
     55/// circumstances, should be negative numbers. That is,
     56/// myKernel->kernel[-3][-2] may be defined if yMin and xMin are equal to or
     57/// more negative than -3 and -2, respectively.
     58///
     59/// In the event that one of the minimum values is greater than the
     60/// corresponding maximum value, the function shall generate a warning, and
     61/// the offending values shall be exchanged.
     62///
     63/// @return psKernel*          A new kernel object
     64///
     65psKernel *psKernelAlloc(
    7366    int xMin,                          ///< Most negative x index
    7467    int xMax,                          ///< Most positive x index
    … …  
    7770);
    7871
    79 /** Checks the type of a particular pointer.
    80  *
    81  * Uses the appropriate deallocation function in psMemBlock to check the ptr datatype.
    82  *
    83  * @return bool:       True if the pointer matches a psKernel structure, false otherwise.
    84  */
     72/// Checks the type of a particular pointer.
     73///
     74/// Uses the appropriate deallocation function in psMemBlock to check the ptr datatype.
     75///
     76/// @return bool:       True if the pointer matches a psKernel structure, false otherwise.
     77///
    8578bool psMemCheckKernel(
    8679    psPtr ptr                          ///< the pointer whose type to check
    … …  
    8881
    8982
    90 /** Generates a kernel given a list of shift values
    91  *
    92  *  Given a list of values (e.g., shifts made in the course of OT guiding),
    93  *  psKernelGenerate shall return the appropriate kernel.  The vectors xShifts
    94  *  and yShifts, which are a list of shifts relative to some starting point,
    95  *  will be supplied by the user. The elements of the vectors should be of an
    96  *  integer type; otherwise the values shall be truncated to integers. The
    97  *  output kernel shall be normalized such that the sum over the kernel is
    98  *  unity.
    99  *
    100  *  If the vectors are not of the same number of elements, then the function
    101  *  shall generate a warning shall be generated, following which, the longer
    102  *  vector trimmed to the length of the shorter, and the function shall continue.
    103  *
    104  *  @return psKernel*    new Kernel object
    105  */
    106 psKernel* psKernelGenerate(
    107     const psVector* tShifts,           ///< list of time shifts
    108     const psVector* xShifts,           ///< list of x-axis shifts
    109     const psVector* yShifts,           ///< list of y-axis shifts
    110     bool relative
    111     /**< specifies the starting point for the shifts; true=relative to previous shift
    112      *  false = relative to some other starting point.  */
     83/// Generates a kernel given a list of shift values
     84///
     85/// Given a list of values (e.g., shifts made in the course of OT guiding),
     86/// psKernelGenerate shall return the appropriate kernel.  The vectors xShifts
     87/// and yShifts, which are a list of shifts relative to some starting point,
     88/// will be supplied by the user. The elements of the vectors should be of an
     89/// integer type; otherwise the values shall be truncated to integers. The
     90/// output kernel shall be normalized such that the sum over the kernel is
     91/// unity.
     92///
     93/// If the vectors are not of the same number of elements, then the function
     94/// shall generate a warning shall be generated, following which, the longer
     95/// vector trimmed to the length of the shorter, and the function shall continue.
     96///
     97/// @return psKernel*    new Kernel object
     98///
     99psKernel *psKernelGenerate(
     100    const psVector *tShifts,           ///< list of time shifts (F32)
     101    const psVector *xShifts,           ///< list of x-axis shifts (S32)
     102    const psVector *yShifts,           ///< list of y-axis shifts (S32)
     103    bool tRelative,                    ///< Are times relative (durations) or absolute?
     104    bool xyRelative                    ///< Are x,y positions relative (shifts) or absolute?
    113105);
    114106
    115 /** convolve an image with a kernel
    116  *
    117  *  Given an input image and the convolution kernel, psImageConvolve shall
    118  *  convolve the input image, in, with the kernel, kernel and return the
    119  *  convolved image, out.
    120  *
    121  *  Two methods shall be available for the convolution: if direct is true,
    122  *  then the convolution shall be performed in real space (appropriate for
    123  *  small kernels); otherwise, the convolution shall be performed using Fast
    124  *  Fourier Transforms (FFTs; appropriate for larger kernels). The latter
    125  *  option involves padding the input image, copying the kernel into an image
    126  *  of the same size as the padded input image, performing an FFT on each,
    127  *  multiplying the FFTs, and performing an inverse FFT before trimming the
    128  *  image back to the original size.
    129  *
    130  *  @return psImage*  resulting image
    131  */
    132 psImage* psImageConvolve(
    133     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    134     const psImage* in,                 ///< the psImage to convolve
    135     const psKernel* kernel,            ///< kernel to colvolve with
    136     bool direct                        ///< specifies method, true=direct convolution, false=fourier
     107/// Convolve an image with a kernel, using a direct convolution
     108///
     109/// This is appropriate for small kernels, where there is no time saving to use FFT method.
     110///
     111/// @return psImage*  resulting image
     112///
     113psImage *psImageConvolveDirect(
     114    psImage *out,                       ///< a psImage to recycle.  If NULL, a new psImage is made.
     115    const psImage *in,                  ///< the psImage to convolve
     116    const psKernel *kernel              ///< kernel to colvolve with
    137117);
    138118
    139 /** Smooths an image by parts using 1D Gaussian independently in x and y.
    140  *
    141  *  Applies a circularly symmetric Gaussian smoothing first in x and then in y
    142  *  directions with just a vector.  This process is 2N faster than 2D convolutions (in general).
    143  *
    144  *  @return bool        TRUE if successful, otherwise FALSE
    145  */
     119/// Convolve an image with a kernel, using the FFT
     120///
     121/// This is appropriate for larger kernels, where the direct convolution is slow.  The input image and kernel
     122/// are suitably padded to avoid wrap-around effects.
     123psImage *psImageConvolveFFT(
     124    psImage *out,                       ///< a psImage to recycle.  If NULL, a new psImage is made.
     125    const psImage *in,                  ///< the psImage to convolve
     126    const psKernel *kernel,             ///< kernel to colvolve with
     127    float pad                           ///< Value to use to pad the input image
     128);
     129
     130/// Smooths an image by parts using 1D Gaussian independently in x and y.
     131///
     132/// Applies a circularly symmetric Gaussian smoothing first in x and then in y
     133/// directions with just a vector.  This process is 2N faster than 2D convolutions (in general).
     134///
     135/// @return bool        TRUE if successful, otherwise FALSE
     136///
    146137bool psImageSmooth(
    147138    psImage *image,                    ///< the image to be smoothed

  • trunk/psLib/test/imageops

    • Property svn:ignore

      •  

        old new  
        2323tap_psImageShift
        2424tap_psImageShiftKernel
         25tap_psImageConvolve
         26tap_psImageConvolve2

  • trunk/psLib/test/imageops/.cvsignore

    r9926 r11703  
    2323tap_psImageShift
    2424tap_psImageShiftKernel
     25tap_psImageConvolve
     26tap_psImageConvolve2

  • trunk/psLib/test/imageops/Makefile.am

    r11691 r11703  
    2121        tap_psImageStructManip \
    2222        tap_psImageConvolve \
     23        tap_psImageConvolve2 \
    2324        tap_psImagePixelExtract
    2425
Note: See TracChangeset for help on using the changeset viewer.