Index: trunk/psLib/src/imageops/psImageConvolve.c =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.c (revision 11153) +++ trunk/psLib/src/imageops/psImageConvolve.c (revision 11703) @@ -1,42 +1,47 @@ -/* @file psImageConvolve.c - * - * @brief Contains FFT transform related functions for psImage. - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.42 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-01-19 04:30:33 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ +/// @file psImageConvolve.c +/// +/// @brief Contains FFT transform related functions for psImage. +/// +/// @author Robert DeSonia, MHPCC +/// @author Paul Price, IfA +/// @author Eugene Magnier, IfA +/// +/// @version $Revision: 1.43 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2007-02-08 04:17:58 $ +/// +/// Copyright 2004-2007 Institute for Astronomy, University of Hawaii +/// #ifdef HAVE_CONFIG_H -# include "config.h" +#include "config.h" #endif #include #include -#include "psImageConvolve.h" -#include "psImageFFT.h" -#include "psImageStructManip.h" -#include "psBinaryOp.h" +#include "psAbort.h" #include "psMemory.h" #include "psLogMsg.h" #include "psError.h" #include "psAssert.h" - - - -#define FOURIER_PADDING 32 /* padding amount in every side of the image for fourier convolution */ - -static void freeKernel(psKernel* ptr); - -psKernel* psKernelAlloc(int xMin, int xMax, int yMin, int yMax) -{ - psKernel* result; - psS32 numRows; - psS32 numCols; - - // following is explicitly spelled out in the SDRS as a requirement +#include "psScalar.h" +#include "psBinaryOp.h" +#include "psImageFFT.h" +#include "psImageStructManip.h" +#include "psImagePixelManip.h" + +#include "psImageConvolve.h" + +static void kernelFree(psKernel *kernel) +{ + if (kernel) { + psFree(kernel->image); + psFree(kernel->p_kernelRows); + } + return; +} + +psKernel *psKernelAlloc(int xMin, int xMax, int yMin, int yMax) +{ + // Check the inputs to make sure max > min; if not, switch. if (yMin > yMax) { psLogMsg(__func__, PS_LOG_WARN, @@ -44,5 +49,5 @@ yMin, yMax); - psS32 temp = yMin; + int temp = yMin; yMin = yMax; yMax = temp; @@ -55,40 +60,33 @@ xMin, xMax); - psS32 temp = xMin; + int temp = xMin; xMin = xMax; xMax = temp; } - numRows = yMax - yMin + 1; - numCols = xMax - xMin + 1; - - result = psAlloc(sizeof(psKernel)); - result->xMin = xMin; - result->xMax = xMax; - result->yMin = yMin; - result->yMax = yMax; - result->image = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL); - memset((result->image->p_rawDataBuffer),0,numCols*numRows*PSELEMTYPE_SIZEOF(PS_TYPE_KERNEL)); - result->p_kernelRows = psAlloc(sizeof(float*)*numRows); - - float** kernelRows = result->p_kernelRows; - float** imageRows = result->image->data.PS_TYPE_KERNEL_DATA; - for (psS32 i = 0; i < numRows; i++) { - kernelRows[i] = imageRows[i] - xMin; - } - result->kernel = kernelRows - yMin; - - psMemSetDeallocator(result,(psFreeFunc)freeKernel); - - return result; -} - -void freeKernel(psKernel* ptr) -{ - if (ptr != NULL) { - psFree(ptr->image); - psFree(ptr->p_kernelRows); - } -} + int numRows = yMax - yMin + 1; // Number of rows for kernel image + int numCols = xMax - xMin + 1; // Number of columns for kernel image + + psKernel *kernel = psAlloc(sizeof(psKernel)); // The kernel, to be returned + psMemSetDeallocator(kernel,(psFreeFunc)kernelFree); + + kernel->xMin = xMin; + kernel->xMax = xMax; + kernel->yMin = yMin; + kernel->yMax = yMax; + kernel->image = psImageAlloc(numCols, numRows, PS_TYPE_KERNEL); + psImageInit(kernel->image, 0.0); + + // Set up indirections, so we can refer to kernel->kernel[-1][-3] for the (-1,-1) element, instead of + // kernel->image[kernel->yMax - kernel->yMin + 1][kernel->yMax - kernel->yMin + 3] (yuk!). + kernel->p_kernelRows = psAlloc(sizeof(float*)*numRows); + for (int i = 0; i < numRows; i++) { + kernel->p_kernelRows[i] = kernel->image->data.PS_TYPE_KERNEL_DATA[i] - xMin; + } + kernel->kernel = kernel->p_kernelRows - yMin; + + return kernel; +} + @@ -96,397 +94,264 @@ { PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc)freeKernel ); -} - - -psKernel* psKernelGenerate(const psVector* tShifts, - const psVector* xShifts, - const psVector* yShifts, - bool relative) -{ - psS32 lastX; - psS32 lastY; - psS32 lastT; - psS32 x; - psS32 y; - psS32 t; - psS32 xMin = 0; - psS32 xMax = 0; - psS32 yMin = 0; - psS32 yMax = 0; - psS32 length = 0; - float normalizeTime = 1.0; // fraction of total time for each shift clock - psKernel* result = NULL; - float** kernel = NULL; - - // got non-NULL vectors? - if (tShifts == NULL || xShifts == NULL || yShifts == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Specified shift vectors can not be NULL.")); + return ( psMemGetDeallocator(ptr) == (psFreeFunc)kernelFree ); +} + + +psKernel *psKernelGenerate(const psVector *tShifts, + const psVector *xShifts, + const psVector *yShifts, + bool tRelative, + bool xyRelative) +{ + PS_ASSERT_VECTOR_NON_NULL(tShifts, NULL); + PS_ASSERT_VECTOR_NON_NULL(xShifts, NULL); + PS_ASSERT_VECTOR_NON_NULL(yShifts, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(tShifts, xShifts, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(tShifts, yShifts, NULL); + PS_ASSERT_VECTOR_TYPE(tShifts, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(xShifts, PS_TYPE_S32, NULL); + PS_ASSERT_VECTOR_TYPE(yShifts, PS_TYPE_S32, NULL); + + // If there are no shifts, the kernel is just a 1 at 0,0 + long num = tShifts->n; // Number of shifts + if (num == 0) { + psKernel *kernel = psKernelAlloc(0,0,0,0); + kernel->kernel[0][0] = 1; + return kernel; + } + + // Get dimensions and scaling + int xMin, xMax, yMin, yMax; // Range of values for kernel + int xLast, yLast; // Last location, for relative shifts + xLast = xMin = xMax = xShifts->data.S32[0]; + yLast = yMin = yMax = yShifts->data.S32[0]; + float tSum = tShifts->data.F32[0]; // Sum of the times + for (long i = 1; i < num; i++) { + int x = xShifts->data.S32[i]; // x position in kernel + int y = yShifts->data.S32[i]; // y position in kernel + if (xyRelative) { + x += xLast; + y += yLast; + xLast = x; + yLast = y; + } + if (x < xMin) { + xMin = x; + } + if (x > xMax) { + xMax = x; + } + if (y < yMin) { + yMin = y; + } + if (y > yMax) { + yMax = y; + } + + if (tRelative) { + tSum += tShifts->data.F32[i]; + } + } + + if (!tRelative) { + // Then the total time is simply the final value + // NB: We assume the counter starts at zero! + tSum = tShifts->data.F32[tShifts->n - 1]; + } + + // One more pass through to set the kernel + psKernel *kernel = psKernelAlloc(xMin, xMax, yMin, yMax); // The kernel + xLast = xShifts->data.S32[0]; + yLast = yShifts->data.S32[0]; + float tLast = 0.0; // Last value for t + for (int i = 0; i < num; i++) { + int x = xShifts->data.S32[i]; // x position in kernel + int y = yShifts->data.S32[i]; // y position in kernel + if (xyRelative) { + x += xLast; + y += yLast; + xLast = x; + yLast = y; + } + float t = tShifts->data.F32[i]; + if (tRelative) { + t -= tLast; + tLast = tShifts->data.F32[i]; + } + + kernel->kernel[y][x] += t; + } + + // Normalise the kernel by the total time (kernel sum should be unity) + psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(1.0 / tSum, PS_TYPE_F32)); + + return kernel; +} + +psImage *psImageConvolveDirect(psImage *out, + const psImage *in, + const psKernel *kernel) +{ + PS_ASSERT_IMAGE_NON_NULL(in, NULL); + PS_ASSERT_IMAGE_TYPE_F32_OR_F64(in, NULL); + PS_ASSERT_PTR_NON_NULL(kernel, NULL); + PS_ASSERT_PTR_NON_NULL(kernel->kernel, NULL); + + // Pull out kernel parameters, for convenience + int xMin = kernel->xMin; + int xMax = kernel->xMax; + int yMin = kernel->yMin; + int yMax = kernel->yMax; + float **kernelData = kernel->kernel; + + int numRows = in->numRows; // Number of rows + int numCols = in->numCols; // Number of columns + +#define SPATIAL_CONVOLVE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + ps##TYPE **inData = in->data.TYPE; \ + out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \ + for (int row = 0; row < numRows; row++) { \ + ps##TYPE *outRow = out->data.TYPE[row]; \ + for (int col = 0; col < numCols; col++) { \ + ps##TYPE pixel = 0.0; \ + for (int kRow = PS_MAX(yMin, -row); kRow <= PS_MIN(yMax, numRows - row - 1); kRow++) { \ + for (int kCol = PS_MAX(xMin, -col); kCol <= PS_MIN(xMax, numCols - col - 1); kCol++) { \ + pixel += kernelData[kRow][kCol] * inData[row + kRow][col + kCol]; \ + } \ + } \ + outRow[col] = pixel; \ + } \ + } \ + } \ + break; + + switch (in->type.type) { + SPATIAL_CONVOLVE_CASE(F32); + SPATIAL_CONVOLVE_CASE(F64); + default: + psAbort(PS_FILE_LINE, "Should never get here: bad type that was asserted on previously."); + } + + return out; +} + +psImage *psImageConvolveFFT(psImage *out, + const psImage *in, + const psKernel *kernel, + float pad) +{ + PS_ASSERT_IMAGE_NON_NULL(in, NULL); + PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, NULL); + PS_ASSERT_PTR_NON_NULL(kernel, NULL); + PS_ASSERT_IMAGE_NON_NULL(kernel->image, NULL); + + // Pull out kernel parameters, for convenience + int xMin = kernel->xMin; + int xMax = kernel->xMax; + int yMin = kernel->yMin; + int yMax = kernel->yMax; + + int numRows = in->numRows; // Number of rows in input image + int numCols = in->numCols; // Number of columns in input image + + // Need to pad the input image to protect from wrap-around effects + if (xMax - xMin > numCols || yMax - yMin > numRows) { + // Cannot pad the image if the kernel is larger. + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + _("Kernel cannot extend further than input image size (%dx%d vs %dx%d)."), + xMax, yMax, numCols, numRows); return NULL; } - - // types match? - if (xShifts->type.type != yShifts->type.type || - tShifts->type.type != xShifts->type.type) { - char* typeXStr; - char* typeYStr; - char* typeTStr; - PS_TYPE_NAME(typeXStr,xShifts->type.type); - PS_TYPE_NAME(typeYStr,yShifts->type.type); - PS_TYPE_NAME(typeTStr,tShifts->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input t-, x-, and y-shift vector types (%s/%s/%s) must match."), - typeTStr, typeXStr, typeYStr); + int paddedCols = numCols + PS_MAX(-xMin, xMax); // Number of columns in padded image + int paddedRows = numRows + PS_MAX(-yMin, yMax); // Number of rows in padded image + + // Generate padded image + psImage *paddedImage = psImageAlloc(paddedCols,paddedRows,in->type.type); // Padded input image + psImageOverlaySection(paddedImage, in, 0, 0, "="); + for (int y = 0; y < numRows; y++) { + for (int x = numCols; x < paddedCols; x++) { + paddedImage->data.F32[y][x] = pad; + } + } + for (int y = numRows; y < paddedRows; y++) { + for (int x = 0; x < paddedCols; x++) { + paddedImage->data.F32[y][x] = pad; + } + } + + // Result of FFT + psImage *inRealFFT = NULL, *inImagFFT = NULL; + if (!psImageForwardFFT(&inRealFFT, &inImagFFT, paddedImage)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to fourier transform input image.")); + psFree(paddedImage); return NULL; } - - // sizes match? - length = xShifts->n; - if (length != yShifts->n || - length != tShifts->n) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "Shift vectors can not be of different sizes."); + psFree(paddedImage); + + // Generate padded kernel image + psImage *paddedKernel = psImageAlloc(paddedCols, paddedRows, PS_TYPE_F32); + psImageInit(paddedKernel, 0.0); + for (int y = PS_MIN(-1, yMin); y <= PS_MIN(-1, yMax); y++) { + // y is negative + for (int x = PS_MIN(-1, xMin); x <= PS_MIN(-1, xMax); x++) { + // x is negative + paddedKernel->data.F32[paddedRows + y][paddedCols + x] = kernel->kernel[y][x]; + } + for (int x = PS_MAX(0, xMin); x <= PS_MAX(0, xMax); x++) { + // x is positive + paddedKernel->data.F32[paddedRows + y][x] = kernel->kernel[y][x]; + } + } + for (int y = PS_MAX(0, yMin); y <= PS_MAX(0, yMax); y++) { + // y is positive + for (int x = PS_MIN(-1, xMin); x <= PS_MIN(-1, xMax); x++) { + // x is negative + paddedKernel->data.F32[y][paddedCols + x] = kernel->kernel[y][x]; + } + for (int x = PS_MAX(0, xMin); x <= PS_MAX(0, xMax); x++) { + // x is positive + paddedKernel->data.F32[y][x] = kernel->kernel[y][x]; + } + } + + psImage *kernelRealFFT = NULL, *kernelImagFFT = NULL; + if (!psImageForwardFFT(&kernelRealFFT, &kernelImagFFT, paddedKernel)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to fourier transform kernel.")); + psFree(inRealFFT); + psFree(inImagFFT); + psFree(paddedKernel); return NULL; } - - // if no shifts, the kernel is just a 1 at 0,0 - if (length < 1) { - result = psKernelAlloc(0,0,0,0); - result->kernel[0][0] = 1; - return result; - } - - #define KERNEL_GENERATE_CASE(TYPE) \ -case PS_TYPE_##TYPE: { \ - ps##TYPE *tShiftData = tShifts->data.TYPE; \ - ps##TYPE *xShiftData = xShifts->data.TYPE; \ - ps##TYPE *yShiftData = yShifts->data.TYPE; \ - lastX = xShiftData[length-1]; \ - lastY = yShiftData[length-1]; \ - lastT = tShiftData[length-1]; \ - \ - for (int lcv = 0; lcv < length; lcv++) { \ - x = lastX - xShiftData[lcv]; \ - y = lastY - yShiftData[lcv]; \ - \ - if (x < xMin) { \ - xMin = x; \ - } else if (x > xMax) { \ - xMax = x; \ - } \ - if (y < yMin) { \ - yMin = y; \ - } else if (y > yMax) { \ - yMax = y; \ - } \ - } \ - \ - normalizeTime = 1.0 / (float)(tShiftData[length-1]); \ - result = psKernelAlloc(xMin,xMax,yMin,yMax); \ - kernel = result->kernel; \ - \ - psS32 prevT = 0; \ - for (int i = 0; i < length; i++) { \ - t = tShiftData[i] - prevT; \ - x = lastX - xShiftData[i]; \ - y = lastY - yShiftData[i]; \ - \ - kernel[y][x] += (float)t / (float)lastT; \ - prevT = tShiftData[i]; \ - } \ - break; \ - } - - #define RELATIVE_KERNEL_GENERATE_CASE(TYPE) \ -case PS_TYPE_##TYPE: { \ - ps##TYPE *tShiftData = tShifts->data.TYPE; \ - ps##TYPE *xShiftData = xShifts->data.TYPE; \ - ps##TYPE *yShiftData = yShifts->data.TYPE; \ - \ - x = 0; \ - y = 0; \ - t = 0; \ - \ - for (int lcv = length-1; lcv >= 0; lcv--) { \ - t += tShiftData[lcv]; \ - \ - if (x < xMin) { \ - xMin = x; \ - } else if (x > xMax) { \ - xMax = x; \ - } \ - if (y < yMin) { \ - yMin = y; \ - } else if (y > yMax) { \ - yMax = y; \ - } \ - x -= xShiftData[lcv]; \ - y -= yShiftData[lcv]; \ - \ - } \ - result = psKernelAlloc(xMin,xMax,yMin,yMax); \ - kernel = result->kernel; \ - \ - normalizeTime = 1.0 / (float)t; \ - x = 0; \ - y = 0; \ - for (psS32 i = length-1; i >= 0; i--) { \ - kernel[y][x] += (float)(tShiftData[i]) * normalizeTime; \ - x -= xShiftData[i]; \ - y -= yShiftData[i]; \ - \ - } \ - break; \ - } - - if (relative) { - switch (xShifts->type.type) { - RELATIVE_KERNEL_GENERATE_CASE(U8); - RELATIVE_KERNEL_GENERATE_CASE(U16); - RELATIVE_KERNEL_GENERATE_CASE(U32); - RELATIVE_KERNEL_GENERATE_CASE(U64); - RELATIVE_KERNEL_GENERATE_CASE(S8); - RELATIVE_KERNEL_GENERATE_CASE(S16); - RELATIVE_KERNEL_GENERATE_CASE(S32); - RELATIVE_KERNEL_GENERATE_CASE(S64); - RELATIVE_KERNEL_GENERATE_CASE(F32); - RELATIVE_KERNEL_GENERATE_CASE(F64); - RELATIVE_KERNEL_GENERATE_CASE(C32); - RELATIVE_KERNEL_GENERATE_CASE(C64); - - default: { - char* typeStr; - PS_TYPE_NAME(typeStr,xShifts->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - } - } - } else { - switch (xShifts->type.type) { - KERNEL_GENERATE_CASE(U8); - KERNEL_GENERATE_CASE(U16); - KERNEL_GENERATE_CASE(U32); - KERNEL_GENERATE_CASE(U64); - KERNEL_GENERATE_CASE(S8); - KERNEL_GENERATE_CASE(S16); - KERNEL_GENERATE_CASE(S32); - KERNEL_GENERATE_CASE(S64); - KERNEL_GENERATE_CASE(F32); - KERNEL_GENERATE_CASE(F64); - KERNEL_GENERATE_CASE(C32); - KERNEL_GENERATE_CASE(C64); - - default: { - char* typeStr; - PS_TYPE_NAME(typeStr,xShifts->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - } - } - } - - return result; -} - -psImage* psImageConvolve(psImage* out, - const psImage* in, - const psKernel* kernel, - bool direct) -{ - if (in == NULL) { - psFree(out); + psFree(paddedKernel); + + // Convolution in fourier domain is just a pixel-wise multiplication + if (!psImageComplexMultiply(&inRealFFT, &inImagFFT, inRealFFT, inImagFFT, kernelRealFFT, kernelImagFFT)) { + psError(PS_ERR_UNKNOWN, false, _("Unable to multiply fourier transformts.")); + psFree(inRealFFT); + psFree(inImagFFT); + psFree(kernelRealFFT); + psFree(kernelImagFFT); return NULL; } - - if (kernel == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Specified psKernel can not be NULL.")); - psFree(out); + psFree(kernelRealFFT); + psFree(kernelImagFFT); + + psImage *paddedConvolved = NULL; // Padded convolved image + if (!psImageBackwardFFT(&paddedConvolved, inRealFFT, inImagFFT, paddedCols)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to invert fourier transform of convolution image.")); + psFree(inRealFFT); + psFree(inImagFFT); return NULL; } - psS32 xMin = kernel->xMin; - psS32 xMax = kernel->xMax; - psS32 yMin = kernel->yMin; - psS32 yMax = kernel->yMax; - float** kData = kernel->kernel; - - // make the output image to the proper size and type - psS32 numRows = in->numRows; - psS32 numCols = in->numCols; - - - - if (direct) { - // spatial convolution - - #define SPATIAL_CONVOLVE_CASE(TYPE) \ - case PS_TYPE_##TYPE: { \ - ps##TYPE** inData = in->data.TYPE; \ - out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \ - for (psS32 row=0;rowdata.TYPE[row]; \ - for (psS32 col=0;col= 0 && row-kRow < numRows) { \ - for (psS32 kCol = xMin; kCol < xMax; kCol++) { \ - if (col-kCol >= 0 && col-kCol < numCols) { \ - pixel += kData[kRow][kCol] * inData[row-kRow][col-kCol]; \ - } \ - } \ - } \ - } \ - outRow[col] = pixel; \ - } \ - } \ - } \ - break; - - switch (in->type.type) { - SPATIAL_CONVOLVE_CASE(U8) - SPATIAL_CONVOLVE_CASE(U16) - SPATIAL_CONVOLVE_CASE(U32) - SPATIAL_CONVOLVE_CASE(U64) - SPATIAL_CONVOLVE_CASE(S8) - SPATIAL_CONVOLVE_CASE(S16) - SPATIAL_CONVOLVE_CASE(S32) - SPATIAL_CONVOLVE_CASE(S64) - SPATIAL_CONVOLVE_CASE(F32) - SPATIAL_CONVOLVE_CASE(F64) - SPATIAL_CONVOLVE_CASE(C32) - SPATIAL_CONVOLVE_CASE(C64) - - default: { - char* typeStr; - PS_TYPE_NAME(typeStr,in->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - psFree(out); - return NULL; - - } - } - - - } else { - // fourier convolution - psS32 paddedCols = numCols+2*FOURIER_PADDING; - psS32 paddedRows = numRows+2*FOURIER_PADDING; - - // check to see if kernel is smaller, otherwise padding it up will fail. - psS32 kRows = kernel->image->numRows; - psS32 kCols = kernel->image->numCols; - if (kRows >= numRows || kCols >= numCols) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - _("Specified psKernel size, %dx%d, can not be larger than input psImage size, %dx%d."), - kCols,kRows, - numCols, numRows); - psFree(out); - return NULL; - } - - // pad the image - psImage* paddedImage = psImageAlloc(paddedCols,paddedRows,in->type.type); - psS32 elementSize = PSELEMTYPE_SIZEOF(in->type.type); - - // zero out padded area on top and bottom - memset(paddedImage->data.U8[0],0,FOURIER_PADDING*paddedCols*elementSize); - memset(paddedImage->data.U8[FOURIER_PADDING+numRows-1],0,FOURIER_PADDING*paddedCols*elementSize); - - // fill in the image-containing rows. - psS32 sidePaddingSize = FOURIER_PADDING*elementSize; - psS32 imageRowSize = numCols*elementSize; - psU8* paddedData = paddedImage->data.U8[FOURIER_PADDING]; - for (psS32 row=0;rowdata.U8[row],imageRowSize); - paddedData += imageRowSize; - // zero out padded area on right edge. - memset(paddedData,0,sidePaddingSize); - paddedData += sidePaddingSize; - } - - // pad the kernel to the same size of paddedImage - psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL); - memset(paddedKernel->data.U8[0],0,sizeof(float)*numCols*numRows); // zero-out image - psS32 yMax = kernel->yMax; - psS32 xMax = kernel->xMax; - for (psS32 row = kernel->yMin; row <= yMax;row++) { - psS32 padRow = row; - if (padRow < 0) { - padRow += paddedRows; - } - float* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow]; - float* kernelRow = kernel->kernel[row]; - for (psS32 col = kernel->xMin; col <= xMax; col++) { - if (col < 0) { - padData[col+paddedCols] = kernelRow[col]; - } else { - padData[col] = kernelRow[col]; - } - } - } - - psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD); - if (kernelFourier == NULL) { - psError(PS_ERR_UNKNOWN, false, - _("Failed to perform a fourier transform of kernel.")); - psFree(out); - return NULL; - } - - psImage* inFourier = psImageFFT(NULL, paddedImage, PS_FFT_FORWARD); - if (inFourier == NULL) { - psError(PS_ERR_UNKNOWN, false, - _("Failed to perform a fourier transform of input image.")); - psFree(out); - return NULL; - } - - // convolution in fourier domain is just a pixel-wise multiplication - for (int row = 0; row < paddedRows; row++) { - psC32* inRow = inFourier->data.C32[row]; - psC32* kRow = kernelFourier->data.C32[row]; - for (int col = 0; col < paddedCols; col++) { - inRow[col] *= kRow[col]; - } - } - - psImage* complexOut = psImageFFT(NULL, inFourier, - PS_FFT_REVERSE); - if (complexOut == NULL) { - psError(PS_ERR_UNKNOWN, false, - _("Failed to perform a fourier transform of input image.")); - psFree(out); - return NULL; - } - - // subset out the padded area now. - psImage* complexOutSansPad = psImageSubset(complexOut, - psRegionSet(FOURIER_PADDING, FOURIER_PADDING+numCols,FOURIER_PADDING,FOURIER_PADDING+numRows)); - - out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); - float factor = 1.0f/(float)paddedCols/(float)paddedRows; - for (psS32 row = 0; row < numRows; row++) { - psF32* outRow = out->data.F32[row]; - psC32* resultRow = complexOutSansPad->data.C32[row]; - for (psS32 col = 0; col < numCols; col++) { - outRow[col] = crealf(resultRow[col])*factor; - } - } - - psFree(complexOut); // frees complexOutSansPad, as it is a child. - psFree(kernelFourier); - psFree(inFourier); - psFree(paddedImage); - psFree(paddedKernel); - - } + psFree(inRealFFT); + psFree(inImagFFT); + + // Trim off the padding, then renormalise (which also does a copy, so there's no parent for the output) + psImage *convolved = psImageSubset(paddedConvolved, psRegionSet(0, numCols, 0, numRows)); + out = (psImage*)psBinaryOp(out, convolved, "*", + psScalarAlloc(1.0 / paddedCols / paddedRows, PS_TYPE_F32)); + psFree(convolved); + psFree(paddedConvolved); return out; @@ -661,5 +526,5 @@ IMAGESMOOTH_CASE(F64); default: { - char* typeStr; + char *typeStr; PS_TYPE_NAME(typeStr,image->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, true,