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

Timestamp:

Aug 28, 2008, 6:04:19 PM (18 years ago)

Author:

Paul Price

Message:

We want a bad/poor mask distinction in the mask convolution for stacking. Again, otherwise the number of bad pixels just explodes. Since convolving a mask with the PSF-matching kernel is something we want to do in both stack rejection and normal convolution, pulled the code that determines the size of the 'bad' region into a separate function in pmSubtraction.c. Then pmStackReject basically does the same thing that's happening in pmSubtractionConvolve/convolveRegion: go through the image by footprints, convolve the footprint and plug the result back in to the original image. Small API change for pmStackReject. This is very effective in reducing the amount of bad pixels in the stacked image.

File:

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

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -300 +300 @@
                                   psImage *subMask, // Subtraction mask
                                   const pmSubtractionKernels *kernels, // Kernels
-                                  psImage *polyValues, // Polynomial values
+                                  const psImage *polyValues, // Polynomial values
                                   float background, // Background value to apply
                                   psRegion region, // Region to convolve
@@ -340 +340 @@
     // 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;
-            }
-        }
-
+        int box = p_pmSubtractionBadRadius(*kernelImage, kernels, polyValues,
+                                           wantDual, poorFrac); // Size of bad box
         if (box > 0) {
             int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
@@ -415 +389 @@
     return;
 }
+
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -430 +406 @@
     psFree(conv);
     return convolved;
+}
+
+int p_pmSubtractionBadRadius(psKernel *preKernel, const pmSubtractionKernels *kernels,
+                             const psImage *polyValues, bool wantDual, float poorFrac)
+{
+    psKernel *kernel;                   // Kernel to use
+    if (!preKernel) {
+        kernel = solvedKernel(NULL, kernels, polyValues, wantDual);
+    } else {
+        kernel = psMemIncrRefCounter(preKernel);
+    }
+    PS_ASSERT_IMAGE_NON_NULL(polyValues, -1);
+
+    int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
+
+    // Determine the threshold between bad and poor
+    double sumKernel2 = 0.0;            // Sum of the kernel-squared
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            sumKernel2 += PS_SQR(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 -= PS_SQR(kernel->kernel[-box][x]) + PS_SQR(kernel->kernel[box][x]);
+        }
+        for (int y = -box + 1; y <= box - 1; y++) {
+            // Note: not doing corners
+            sumBox -= PS_SQR(kernel->kernel[y][-box]) + PS_SQR(kernel->kernel[y][box]);
+        }
+        if (sumBox < threshold) {
+            break;
+        }
+    }
+
+    psFree(kernel);
+
+    return box;
+}
+
+psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels,
+                                             int numCols, int numRows, int x, int y)
+{
+    assert(kernels);
+    assert(numCols > 0 && numRows > 0);
+
+    // Size to use when calculating normalised coordinates (different from actual size when convolving
+    // subimage)
+    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
+    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
+
+    // Normalised coordinates
+    float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize;
+    float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize;
+
+    return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm);
 }
 
@@ -847 +884 @@
     int xMin = region->x0, xMax = region->x1, yMin = region->y0, yMax = region->y1; // Bounds of patch
 
-    // Size to use when calculating normalised coordinates (different from actual size when convolving
-    // subimage)
-    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
-    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
-
-    // Normalised coordinates
-    float yNorm = 2.0 * (float)(yMin + y0 + size + 1 - yNormSize/2.0) / (float)yNormSize; // Normalised coord
-    float xNorm = 2.0 * (float)(xMin + x0 + size + 1 - xNormSize/2.0) / (float)xNormSize; // Normalised coord
-
     psKernel *kernelImage = NULL;       // Kernel for the images
     psKernel *kernelWeight = NULL;      // Kernel for the weight maps
@@ -861 +889 @@
     // Only generate polynomial values every kernel footprint, since we have already assumed
     // (with the stamps) that it does not vary rapidly on this scale.
-    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder,
-                                                    xNorm, yNorm); // Pre-calculated polynomial values
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
+                                                              xMin + x0 + size + 1,
+                                                              yMin + y0 + size + 1);
     float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term