Changeset 17320 for trunk/psLib/src/imageops/psImageConvolve.c

Timestamp:

Apr 4, 2008, 12:44:56 PM (18 years ago)

Author:

Paul Price

Message:

Optimising psImageConvolveFFT. Made a couple of optimisations: 1. Moved code into psImageFFT.c, so I can use FFTW directly. FFTW, through its 'advanced interface' allows two FFTs to be performed at the same time and 'the resulting plans can often be faster than calling FFTW multiple times for the individual transforms' (FFTW 3.1.2 manual, p30). 2. The convolved image needs to be normalised because FFTW doesn't take out the sqrt(N) factors in the FFT. Instead of multiplying the entire convolved image, we can multiply the much smaller number of pixels in the kernel. Tested this code with tap_psImageConvolve2 and it seems to work. Running convolutionBench unoptimised on mithrandir produces times that are comparable to what I got a while back running optimised. Running it optimized on alala produces speed-ups ranging from 25 to 350%.

File:

• trunk/psLib/src/imageops/psImageConvolve.c (modified) (3 diffs)

trunk/psLib/src/imageops/psImageConvolve.c

@@ -7 +7 @@
 /// @author Eugene Magnier, IfA
 ///
-/// @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
-/// @date $Date: 2008-04-03 03:00:25 $
+/// @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
+/// @date $Date: 2008-04-04 22:44:56 $
 ///
 /// Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -386 +386 @@
 
 
-psImage *psImageConvolveFFT(const psImage *in,
-                            const psImage *mask,
-                            psMaskType maskVal,
-                            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;
-    }
-    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
-    if (mask && maskVal) {
-        // Need to replace non-finite (assumed masked) pixels, since they propagate everywhere during FFT
-        for (int y = 0; y < numRows; y++) {
-            for (int x = 0; x < numCols; x++) {
-                paddedImage->data.F32[y][x] = (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? pad :
-                    in->data.F32[y][x];
-            }
-        }
-    } else {
-        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;
-    }
-    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
-        if (xMin < 0) {
-            // x is negative
-            memcpy(&paddedKernel->data.F32[paddedRows + y][paddedCols + xMin], &kernel->kernel[y][xMin],
-                   (PS_MIN(0, xMax) - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
-        }
-        if (xMax >= 0) {
-            // x is positive
-            int min = PS_MAX(0, xMin);  // Minimum value of x when positive
-            memcpy(&paddedKernel->data.F32[paddedRows + y][min], &kernel->kernel[y][min],
-                   (xMax - min + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
-        }
-    }
-    for (int y = PS_MAX(0, yMin); y <= PS_MAX(0, yMax); y++) {
-        // y is positive
-        if (xMin < 0) {
-            // x is negative
-            memcpy(&paddedKernel->data.F32[y][paddedCols + xMin], &kernel->kernel[y][xMin],
-                   (PS_MIN(0, xMax) - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
-        }
-        if (xMax >= 0) {
-            // x is positive
-            int min = PS_MAX(0, xMin);  // Minimum value of x when positive
-            memcpy(&paddedKernel->data.F32[y][min], &kernel->kernel[y][min],
-                   (xMax - min + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
-        }
-    }
-
-    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;
-    }
-    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;
-    }
-    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;
-    }
-    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));
-    psImage *out = (psImage*)psBinaryOp(NULL, convolved, "*",
-                                        psScalarAlloc(1.0 / paddedCols / paddedRows, PS_TYPE_F32));
-    psFree(convolved);
-    psFree(paddedConvolved);
-
-    return out;
-}
-
-
 psImage *psImageConvolveMaskFFT(psImage *out, const psImage *mask, psMaskType maskVal,
                                 psMaskType setVal, int xMin, int xMax, int yMin, int yMax, float thresh)
@@ -573 +434 @@
     psKernel *kernel = psKernelAlloc(xMin, xMax, yMin, yMax);
     psImageInit(kernel->image, 1.0);
-    psImage *convolved = psImageConvolveFFT(onoff, NULL, 0, kernel, 0.0);
+    psImage *convolved = psImageConvolveFFT(NULL, onoff, NULL, 0, kernel);
     psFree(onoff);
     psFree(kernel);