IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/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,
Note: See TracChangeset for help on using the changeset viewer.