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

Timestamp:

Oct 3, 2008, 4:37:51 PM (18 years ago)

Author:

Paul Price

Message:

Adding psImageSmoothMask: a faster version of psImageSmoothMaskF32 (by about 30% when running optimised on ipp002).

File:

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

trunk/psLib/src/imageops/psImageConvolve.c

@@ -7 +7 @@
 /// @author Eugene Magnier, IfA
 ///
-/// @version $Revision: 1.76 $ $Name: not supported by cvs2svn $
-/// @date $Date: 2008-08-30 02:24:21 $
+/// @version $Revision: 1.77 $ $Name: not supported by cvs2svn $
+/// @date $Date: 2008-10-04 02:37:51 $
 ///
 /// Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -33 +33 @@
 
 #include "psImageConvolve.h"
+
+#define MIN_GAUSS_FRAC 0.25             // Minimum Gaussian fraction to accept when smoothing
+
 
 static bool threaded = false;           // Run image convolution threaded?
@@ -642 +645 @@
 }
 
+// Smooth an image with masked pixels
+// The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
+// [col][row] instead of [row][col]) to optimise the iteration over rows; such instances are marked "BW"
+#define IMAGESMOOTH_MASK_CASE(TYPE) \
+case PS_TYPE_##TYPE: { \
+    psImage *calculation = psImageAlloc(numRows, numCols, PS_TYPE_##TYPE); /* Calculation image; BW */ \
+    psImage *calcMask = psImageAlloc(numRows, numCols, PS_TYPE_MASK); /* Mask for calculation image; BW */ \
+    \
+    /** Smooth in X direction **/ \
+    for (int j = 0; j < numRows; j++) { \
+        for (int i = 0; i < numCols; i++) { \
+            int xMin = PS_MAX(i - size, 0); \
+            int xMax = PS_MIN(i + size, xLast); \
+            const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[j][xMin]; \
+            const ps##TYPE *imageData = &image->data.TYPE[j][xMin]; \
+            int uMin = - PS_MIN(i, size); /* Minimum kernel index */ \
+            const psF32 *gaussData = &gauss[uMin]; \
+            double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */ \
+            for (int x = xMin; x <= xMax; x++, maskData++, imageData++, gaussData++) { \
+                if (*maskData & maskVal) { \
+                    continue; \
+                } \
+                sumIG += *imageData * *gaussData; \
+                sumG += *gaussData; \
+            } \
+            if (sumG > MIN_GAUSS_FRAC) { \
+                /* BW */ \
+                calculation->data.TYPE[i][j] = sumIG / sumG; \
+                calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0; \
+            } else { \
+                /* BW */ \
+                calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0xFF; \
+            } \
+        } \
+    } \
+    \
+    output = psImageRecycle(output, numCols, numRows, PS_TYPE_##TYPE); \
+    \
+    /** Smooth in Y direction  **/ \
+    for (int i = 0; i < numCols; i++) { \
+        for (int j = 0; j < numRows; j++) { \
+            int yMin = PS_MAX(j - size, 0); \
+            int yMax = PS_MIN(j + size, yLast); \
+            const psMaskType *maskData = &calcMask->data.PS_TYPE_MASK_DATA[i][yMin]; /* BW */ \
+            const ps##TYPE *imageData = &calculation->data.TYPE[i][yMin]; /* BW */ \
+            int vMin = - PS_MIN(j, size); /* Minimum kernel index */ \
+            const psF32 *gaussData = &gauss[vMin]; \
+            double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */ \
+            for (int y = yMin; y <= yMax; y++, maskData++, imageData++, gaussData++) { \
+                if (*maskData) { \
+                    continue; \
+                } \
+                sumIG += *imageData * *gaussData; \
+                sumG += *gaussData; \
+            } \
+            output->data.TYPE[j][i] = (sumG > MIN_GAUSS_FRAC) ? sumIG / sumG : NAN; \
+        } \
+    } \
+    \
+    psFree(calculation); \
+    psFree(calcMask); \
+    psFree(gaussNorm); \
+    \
+    return output; \
+}
+
+
+
+psImage *psImageSmoothMask(psImage *output,
+                           const psImage *image,
+                           const psImage *mask,
+                           psMaskType maskVal,
+                           double sigma,
+                           double numSigma)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
+    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL);
+
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int xLast = numCols - 1, yLast = numRows - 1; // Last index
+
+    // relevant terms
+    int size = sigma * numSigma + 0.5;  // Half-size of kernel
+    int fullSize = 2*size + 1;          // Total number of pixels in convolution kernel
+
+    // Generate normalized gaussian
+    psVector *gaussNorm = psVectorAlloc(fullSize + 1, PS_TYPE_F32); // Normalised Gaussian
+    {
+        double sum = 0.0;
+        double norm = -0.5/(sigma*sigma);
+        for (int i = 0, u = -size; i <= fullSize; i++, u++) {
+            float value = expf(norm*PS_SQR(u));
+            gaussNorm->data.F32[i] = value;
+            sum += value;
+        }
+        sum = 1.0 / sum;
+        for (int i = 0; i <= fullSize; i++) {
+            gaussNorm->data.F32[i] *= sum;
+        }
+    }
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+
+    switch (image->type.type) {
+        IMAGESMOOTH_MASK_CASE(F32);
+        IMAGESMOOTH_MASK_CASE(F64);
+      default:
+        psFree(gaussNorm);
+        char *typeStr;
+        PS_TYPE_NAME(typeStr,image->type.type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                _("Specified psImage type, %s, is not supported."),
+                typeStr);
+        return false;
+    }
+
+    psAbort("Should never reach here.");
+    return false;
+}
+
+
 bool psImageSmoothMaskF32 (psImage *image,
                            psImage *mask,
@@ -695 +820 @@
                 g += *sg;
             }
-            *vo = (g > 0.25) ? s / g : NAN;
+            *vo = (g > MIN_GAUSS_FRAC) ? s / g : NAN;
         }
         memcpy(image->data.F32[j], calculation->data.F32, Nx*sizeof(psF32));
@@ -740 +865 @@
         psF32 *vs = outsum->data.F32;
         for (int i = 0; i < Nx; i++, vo++, vs++) {
-            *vo = (*vs > 0.25) ? *vo / *vs : NAN;
+            *vo = (*vs > MIN_GAUSS_FRAC) ? *vo / *vs : NAN;
         }