Changeset 19164 for trunk/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Aug 22, 2008, 12:45:36 PM (18 years ago)

Author:

Paul Price

Message:

Reworking the masking so as to avoid exploding the number of bad pixels in the convolved image. Introduced the bad/poor distinction. Bad pixels that are convolved now have a large 'poor' halo in the mask, with a smaller 'bad' halo. Got rid of the 'FOOTPRINT' mask value (internal subtraction mask): there's no need to convolve the mask in order to find out if a sprinkling of footprints are bad (not sure this is so efficient when a list of sources isn't supplied), and I needed the mask bits to signify the difference between 'bad' and 'poor'. Tested on a single subtraction, which seems to work.

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (14 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -249 +249 @@
     psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, size, -size)); // Cut off the edges
 
-    if (threaded) {
-        psMutexLock(target);
-    }
-    psImage *subTarget = psImageSubset(target, region); // Target for this subregion
-    if (threaded) {
-        psMutexUnlock(target);
-    }
+    int xMin = region.x0, xMax = region.x1, yMin = region.y0, yMax = region.y1; // Bounds of region
     if (background != 0.0) {
-        psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32));
+        for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
+            for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) {
+                target->data.F32[yTarget][xTarget] = convolved->data.F32[ySource][xSource] + background;
+            }
+        }
     } else {
-        int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
-        for (int y = 0; y < subTarget->numRows; y++) {
-            memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes);
+        int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
+        for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
+            memcpy(&target->data.F32[yTarget][xMin], &convolved->data.F32[ySource][size], numBytes);
         }
     }
     psFree(subConv);
     psFree(convolved);
-    if (threaded) {
-        psMutexLock(target);
-    }
-    psFree(subTarget);
-    if (threaded) {
-        psMutexUnlock(target);
-    }
 
     return;
@@ -291 +282 @@
             for (int v = -size; v <= size; v++) {
                 for (int u = -size; u <= size; u++) {
-                    target->data.F32[y][x] += kernel->kernel[v][u] *
-                        image->data.F32[y - v][x - u];
+                    target->data.F32[y][x] += kernel->kernel[v][u] * image->data.F32[y - v][x - u];
                 }
             }
@@ -303 +293 @@
 static inline void convolveRegion(psImage *convImage, // Convolved image (output)
                                   psImage *convWeight, // Convolved weight map (output), or NULL
+                                  psImage *convMask, // Convolve mask (output), or NULL
                                   psKernel **kernelImage, // Convolution kernel for the image
                                   psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL
@@ -308 +299 @@
                                   psImage *weight, // Weight map to convolve, or NULL
                                   psImage *subMask, // Subtraction mask
-                                  psMaskType maskVal, // Mask value to apply in convolution
                                   const pmSubtractionKernels *kernels, // Kernels
                                   psImage *polyValues, // Polynomial values
                                   float background, // Background value to apply
                                   psRegion region, // Region to convolve
+                                  float poorFrac, // Fraction for "poor"
                                   bool useFFT,  // Use FFT to convolve?
                                   bool wantDual // Want the dual convolution?
@@ -318 +309 @@
 {
     *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual);
-    if (weight) {
+    if (weight || subMask) {
         *kernelWeight = varianceKernel(*kernelWeight, *kernelImage);
     }
 
+    psMaskType maskSource;              // Mask these values when
+    psMaskType maskTarget;              // Set this mask value when convolving
+    if (kernels->mode == PM_SUBTRACTION_MODE_1 || (kernels->mode == PM_SUBTRACTION_MODE_DUAL && !wantDual)) {
+        maskSource = PM_SUBTRACTION_MASK_BAD_1;
+        maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1;
+    } else {
+        maskSource = PM_SUBTRACTION_MASK_BAD_2;
+        maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_BAD_2;
+    }
+
+    // Convolve the image and weight
     if (useFFT) {
         // Use Fast Fourier Transform to do the convolution
         // This provides a big speed-up for large kernels
-
-        convolveFFT(convImage, image, subMask, maskVal, *kernelImage, region, background, kernels->size);
+        convolveFFT(convImage, image, subMask, maskSource, *kernelImage, region, background, kernels->size);
         if (weight) {
-            convolveFFT(convWeight, weight, subMask, maskVal, *kernelWeight, region, 0.0, kernels->size);
+            convolveFFT(convWeight, weight, subMask, maskSource, *kernelWeight, region, 0.0, kernels->size);
         }
     } else {
@@ -334 +335 @@
         if (weight) {
             convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size);
+        }
+    }
+
+    // Convolve the mask for bad pixels
+    if (subMask && convMask) {
+        psKernel *kernel = *kernelWeight; // Kernel of interest
+        int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
+
+        // Determine the thresholds
+        double sumKernel2 = 0.0;        // Sum of the kernel-squared
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                sumKernel2 += kernel->kernel[y][x];
+            }
+        }
+        float threshold = sumKernel2 * poorFrac; // Threshold between poor and bad
+
+        // Get bounds of threshold region
+        // Start with the entire kernel, and keep reducing the size of the box until it drops below threshold
+        int box = kernels->size;                    // Size of box with bad pixels
+        for (double sumBox = sumKernel2; box > 0; box--) {
+            for (int x = -box; x <= box; x++) {
+                sumBox -= kernel->kernel[-box][x] + kernel->kernel[box][x];
+            }
+            for (int y = -box + 1; y <= box - 1; y++) {
+                // Note: not doing corners
+                sumBox -= kernel->kernel[y][-box] + kernel->kernel[y][box];
+            }
+            if (sumBox < threshold) {
+                break;
+            }
+        }
+
+        if (box > 0) {
+            bool threaded = pmSubtractionThreaded(); // Are we running threaded?
+            if (threaded) {
+                psMutexLock(subMask);
+            }
+            psImage *image = psImageSubset(subMask, psRegionSet(region.x0 + xMin, region.x1 + xMax,
+                                                                region.y0 + yMin, region.y1 + yMax));
+            if (threaded) {
+                psMutexUnlock(subMask);
+            }
+
+            psImage *convolved = psImageConvolveMask(NULL, image, maskSource, maskTarget,
+                                                     -box, box, -box, box); // Convolved subtraction mask
+
+            if (threaded) {
+                psMutexLock(subMask);
+            }
+            psFree(image);
+            if (threaded) {
+                psMutexUnlock(subMask);
+            }
+
+            int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
+            int numBytes = (colMax - colMin) * PSELEMTYPE_SIZEOF(PS_TYPE_MASK); // Number of bytes to copy
+
+            // Since we're copying back into the source, we need to lock
+            if (threaded) {
+                psMutexLock(subMask);
+            }
+            for (int yTarget = rowMin, ySource = -yMin; yTarget < rowMax; yTarget++, ySource++) {
+                memcpy(&subMask->data.PS_TYPE_MASK_DATA[yTarget][colMin],
+                       &convolved->data.PS_TYPE_MASK_DATA[ySource][-xMin], numBytes);
+            }
+            if (threaded) {
+                psMutexUnlock(subMask);
+            }
+
+            // No need to lock: we own this
+            psFree(convolved);
         }
     }
@@ -753 +826 @@
 
 // XXX Put kernelImage, kernelWeight and polyValues on thread-dependent data
-static bool subtractionConvolvePatch(int numCols, int numRows, int x0, int y0,
-                                     pmReadout *out1, pmReadout *out2, psImage *convMask,
-                                     const pmReadout *ro1, const pmReadout *ro2, psImage *subMask,
-                                     psMaskType blank, psMaskType maskSource, psMaskType maskTarget,
-                                     const psRegion *region, const pmSubtractionKernels *kernels,
-                                     bool doBG, bool useFFT)
+static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image
+                                     int x0, int y0, // Offsets for image
+                                     pmReadout *out1, pmReadout *out2, // Output readouts
+                                     psImage *convMask, // Output convolved mask
+                                     const pmReadout *ro1, const pmReadout *ro2, // Input readouts
+                                     psImage *subMask, // Input subtraction mask
+                                     psMaskType maskBad, // Mask value to give bad pixels
+                                     psMaskType maskPoor, // Mask value to give poor pixels
+                                     float poorFrac, // Fraction for "poor"
+                                     const psRegion *region, // Patch to convolve
+                                     const pmSubtractionKernels *kernels, // Kernels
+                                     bool doBG, // Add in background when convolving?
+                                     bool useFFT // Use FFT to do the convolution?
+    )
 {
     int size = kernels->size;           // Half-size of kernel
@@ -782 +863 @@
 
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        convolveRegion(out1->image, out1->weight, &kernelImage, &kernelWeight, ro1->image, ro1->weight,
-                       subMask, maskSource, kernels, polyValues, background, *region, useFFT,
-                       false);
+        convolveRegion(out1->image, out1->weight, convMask, &kernelImage, &kernelWeight,
+                       ro1->image, ro1->weight, subMask, kernels, polyValues, background,
+                       *region, poorFrac, useFFT, false);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        convolveRegion(out2->image, out2->weight, &kernelImage, &kernelWeight, ro2->image, ro2->weight,
-                       subMask, maskSource, kernels, polyValues, background, *region, useFFT,
-                       kernels->mode == PM_SUBTRACTION_MODE_DUAL);
+        convolveRegion(out2->image, out2->weight, convMask, &kernelImage, &kernelWeight,
+                       ro2->image, ro2->weight, subMask, kernels, polyValues, background,
+                       *region, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL);
     }
 
@@ -796 +877 @@
     psFree(polyValues);
 
-    // Propagate the mask
+    // Propagate the subtraction mask
     if (subMask) {
         psMaskType **target = convMask->data.PS_TYPE_MASK_DATA; // Target mask
         psMaskType **source = subMask->data.PS_TYPE_MASK_DATA; // Source mask
 
+        psMaskType poor, bad;           // Mask value for "poor" and "bad" pixels in subtraction mask
+        switch (kernels->mode) {
+          case PM_SUBTRACTION_MODE_1:
+            poor = PM_SUBTRACTION_MASK_CONVOLVE_1;
+            bad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
+            break;
+          case PM_SUBTRACTION_MODE_2:
+            poor = PM_SUBTRACTION_MASK_CONVOLVE_2;
+            bad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
+            break;
+          case PM_SUBTRACTION_MODE_DUAL:
+            poor = PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_1;
+            bad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 |
+                PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
+            break;
+          default:
+            psAbort("Unrecognised subtraction mode: %x", kernels->mode);
+        }
+
         for (int y = yMin; y < yMax; y++) {
             for (int x = xMin; x < xMax; x++) {
-                if (source[y][x] & maskTarget) {
-                    target[y][x] |= blank;
+                // Pixels marked as "bad" shouldn't also be marked as "poor"
+                if (source[y][x] & bad) {
+                    target[y][x] |= maskBad;
+                } else if (source[y][x] & poor){
+                    target[y][x] |= maskPoor;
                 }
             }
@@ -848 +951 @@
     const pmReadout *ro2 = args->data[8]; // Input readout 2
     psImage *subMask = args->data[9]; // Subtraction mask
-    psMaskType blank = PS_SCALAR_VALUE(args->data[10], U8); // Output mask value
-    psMaskType maskSource = PS_SCALAR_VALUE(args->data[11], U8); // Mask value corresponding to source image
-    psMaskType maskTarget = PS_SCALAR_VALUE(args->data[12], U8); // Mask value corresponding to target image
+    psMaskType maskBad = PS_SCALAR_VALUE(args->data[10], U8); // Output mask value for bad pixels
+    psMaskType maskPoor = PS_SCALAR_VALUE(args->data[11], U8); // Output mask value for poor pixels
+    float poorFrac = PS_SCALAR_VALUE(args->data[12], F32); // Fraction for "poor"
     const psRegion *region = args->data[13]; // Region to convolve
     const pmSubtractionKernels *kernels = args->data[14]; // Kernels
@@ -856 +959 @@
     bool useFFT = PS_SCALAR_VALUE(args->data[16], U8); // Use FFT for convolution?
 
-    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask,
-                                    ro1, ro2, subMask, blank, maskSource, maskTarget,
-                                    region, kernels, doBG, useFFT);
+    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, subMask,
+                                    maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT);
 }
 
 bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2,
-                           psImage *subMask, psMaskType blank, const psRegion *region,
-                           const pmSubtractionKernels *kernels, bool doBG, bool useFFT)
+                           psImage *subMask, psMaskType maskBad, psMaskType maskPoor, float poorFrac,
+                           const psRegion *region, const pmSubtractionKernels *kernels,
+                           bool doBG, bool useFFT)
 {
     int numCols = 0, numRows = 0;       // Image dimensions
@@ -998 +1101 @@
         yMin = PS_MAX(region->y0, yMin);
         yMax = PS_MIN(region->y1, yMax);
-    }
-
-    psMaskType maskSource;              // Mask these pixels when convolving
-    psMaskType maskTarget;              // Mark these pixels as bad when propagating the subtractionMask
-    switch (kernels->mode) {
-      case PM_SUBTRACTION_MODE_1:
-        maskSource = PM_SUBTRACTION_MASK_BAD_1;
-        maskTarget = PM_SUBTRACTION_MASK_BAD_2 | PM_SUBTRACTION_MASK_CONVOLVE_1;
-        break;
-      case PM_SUBTRACTION_MODE_2:
-        maskSource = PM_SUBTRACTION_MASK_BAD_2;
-        maskTarget = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
-        break;
-      case PM_SUBTRACTION_MODE_DUAL:
-        maskSource = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
-        maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
-        break;
-      default:
-        psAbort("Unsupported subtraction mode: %x", kernels->mode);
     }
 
@@ -1045 +1129 @@
                 psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const
                 psArrayAdd(args, 1, (psImage*)subMask); // Casting away const
-                PS_ARRAY_ADD_SCALAR(args, blank, PS_TYPE_U8);
-                PS_ARRAY_ADD_SCALAR(args, maskSource, PS_TYPE_U8);
-                PS_ARRAY_ADD_SCALAR(args, maskTarget, PS_TYPE_U8);
+                PS_ARRAY_ADD_SCALAR(args, maskBad, PS_TYPE_U8);
+                PS_ARRAY_ADD_SCALAR(args, maskPoor, PS_TYPE_U8);
+                PS_ARRAY_ADD_SCALAR(args, poorFrac, PS_TYPE_F32);
                 psArrayAdd(args, 1, subRegion);
                 psArrayAdd(args, 1, (pmSubtractionKernels*)kernels); // Casting away const
@@ -1060 +1144 @@
             } else {
                 subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, subMask,
-                                         blank, maskSource, maskTarget, subRegion, kernels, doBG, useFFT);
+                                         maskBad, maskPoor, poorFrac, subRegion, kernels, doBG, useFFT);
             }
             psFree(subRegion);