Changeset 19282 for trunk/psModules/src/imcombine/pmStackReject.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/pmStackReject.c (modified) (9 diffs)

trunk/psModules/src/imcombine/pmStackReject.c

@@ -4 +4 @@
 
 #include <stdio.h>
+#include <string.h>
 #include <pslib.h>
 
@@ -14 +15 @@
 
 
-psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold,
-                        const psArray *subRegions, const psArray *kernels)
+psPixels *pmStackReject(const psPixels *in, int numCols, int numRows, float threshold, float poorFrac,
+                        const psArray *subRegions, const psArray *subKernels)
 {
     PS_ASSERT_PIXELS_NON_NULL(in, NULL);
@@ -21 +22 @@
     PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(threshold, 1.0, NULL);
     PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL);
-    PS_ASSERT_ARRAY_NON_NULL(kernels, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL);
+    PS_ASSERT_ARRAY_NON_NULL(subKernels, NULL);
+    PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, subKernels, NULL);
 
     // Trivial case
@@ -29 +30 @@
     }
 
-    // Get the original image size
-    int numRegions = subRegions->n;        // Number of regions
-    int numCols = 0, numRows = 0;       // Size of original image
-    int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image
+    // Check consistency of kernels
+    int numRegions = subRegions->n;     // Number of regions
     int size = 0;                       // Size of kernel
     for (int i = 0; i < numRegions; i++) {
-        psRegion *subRegion = subRegions->data[i]; // Region of interest
-        if (subRegion->x0 < minCols) {
-            minCols = subRegion->x0;
-        }
-        if (subRegion->y0 < minRows) {
-            minRows = subRegion->y0;
-        }
-        if (subRegion->x1 > numCols) {
-            numCols = subRegion->x1;
-        }
-        if (subRegion->y1 > numRows) {
-            numRows = subRegion->y1;
-        }
-
-        pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
+        pmSubtractionKernels *kernels = subKernels->data[i]; // Kernel of interest
         if (size == 0) {
-            size = kernel->size;
-        } else if (kernel->size != size) {
+            size = kernels->size;
+        } else if (kernels->size != size) {
             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Kernel sizes are not identical: %d vs %d",
-                    size, kernel->size);
+                    size, kernels->size);
             return NULL;
         }
     }
 
-    // Adjust the size for the size of the subimage
-    if (valid) {
-        minCols = PS_MAX(valid->x0, minCols);
-        minRows = PS_MAX(valid->y0, minRows);
-        numCols = PS_MIN(valid->x1, numCols);
-        numRows = PS_MIN(valid->y1, numRows);
-    }
-    psTrace("psModules.imcombine", 1, "Rejecting [%d:%d,%d:%d]\n", minCols, numCols, minRows, numRows);
-
-    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1),
-                                   0x01); // Mask
+    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols - 1, 0, numRows - 1), 0x01); // Mask
     psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve
     psFree(mask);
@@ -77 +52 @@
     pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
     inRO->image = image;
-    inRO->col0 = minCols;
-    inRO->row0 = minRows;
     for (int i = 0; i < numRegions; i++) {
         psRegion *region = subRegions->data[i]; // Region of interest
-        if (valid && (region->x0 > valid->x1 || region->x1 < valid->x0 ||
-                      region->y0 > valid->y1 || region->y1 < valid->y0)) {
-            // Region is outside of our sub-image
-            continue;
-        }
-        pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
-        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, 0, 1.0, region, kernel, false, true)) {
+        pmSubtractionKernels *kernels = subKernels->data[i]; // Kernel of interest
+        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, 0, 1.0, region, kernels, false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
             psFree(convRO);
@@ -98 +66 @@
 
         // Image of the kernel at the centre of the region
-        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernel->numCols/2.0) /
-            (float)kernel->numCols;
-        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernel->numRows/2.0) /
-            (float)kernel->numRows;
-        psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false);
+        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernels->numCols/2.0) /
+            (float)kernels->numCols;
+        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernels->numRows/2.0) /
+            (float)kernels->numRows;
+        psImage *image = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
         if (!image) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
@@ -118 +86 @@
 
         // Range for normalisation
-        int yMin = PS_MAX(minRows, region->y0) - inRO->row0;
-        int yMax = PS_MIN(numRows - 1, region->y1) - inRO->row0;
-        int xMin = PS_MAX(minCols, region->x0) - inRO->col0;
-        int xMax = PS_MIN(numCols - 1, region->x1) - inRO->col0;
+        int xMin = PS_MAX(0, region->x0), xMax = PS_MIN(numCols - 1, region->x1);
+        int yMin = PS_MAX(0, region->y0), yMax = PS_MIN(numRows - 1, region->y1);
         psTrace("psModules.imcombine", 2, "Normalising convolved mask image by %f over %d:%d,%d:%d\n",
                 sum, xMin, xMax, yMin, yMax);
+        sum = 1.0 / sum;
         for (int y = yMin; y <= yMax; y++) {
             for (int x = xMin; x <= xMax; x++) {
-                convRO->image->data.F32[y][x] /= sum;
+                convRO->image->data.F32[y][x] *= sum;
             }
         }
@@ -152 +119 @@
         for (int x = size; x < convolved->numCols - size; x++) {
             if (convolved->data.F32[y][x] > threshold) {
-                // Pixel coordinates in "bad" correspond to the full image
-                bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x + minCols, y + minRows);
+                bad = psPixelsAdd(bad, bad->nalloc, x, y);
             }
         }
@@ -160 +126 @@
 
     // Now, grow the mask to include everything that touches a bad pixel in the convolution
-    int x0 = minCols, y0 = minRows;     // Offset for mask image
-    mask = psPixelsToMask(NULL, bad, psRegionSet(x0, numCols - 1, y0, numRows - 1), 0xff);
-    for (int i = 0; i < bad->n; i++) {
-        int xPix = bad->data[i].x - x0, yPix = bad->data[i].y - y0; // Coordinates in frame of mask image
-        // Convolution limits
-        int xMin = PS_MAX(xPix - size, 0);
-        int xMax = PS_MIN(xPix + size, mask->numCols - 1);
-        int yMin = PS_MAX(yPix - size, 0);
-        int yMax = PS_MIN(yPix + size, mask->numRows - 1);
-        for (int y = yMin; y <= yMax; y++) {
-            for (int x = xMin; x <= xMax; x++) {
-                assert(x < mask->numCols && y < mask->numRows);
-                mask->data.PS_TYPE_MASK_DATA[y][x] = 0xff;
+    // Mask values:
+    // 0x0f: we think this is bad
+    // 0xf0: this is within the convolution kernel of a bad pixel
+    mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0x0f);
+    for (int i = 0; i < subRegions->n; i++) {
+        psRegion *region = subRegions->data[i]; // Subtraction region
+        pmSubtractionKernels *kernels = subKernels->data[i]; // Subtraction kernel
+
+        int size = kernels->size;           // Half-size of kernel
+        int fullSize = 2 * size + 1;        // Full size of kernel
+
+        // Get region for convolution: [xMin:xMax,yMin:yMax]
+        int xMin = PS_MAX(region->x0, size), xMax = PS_MIN(region->x1, numCols - size);
+        int yMin = PS_MAX(region->y0, size), yMax = PS_MIN(region->y1, numRows - size);
+
+        psImage *polyValues = NULL;     // Pre-calculated polynomial values
+        for (int j = yMin; j < yMax; j += fullSize) {
+            int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
+            for (int i = xMin; i < xMax; i += fullSize) {
+                int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
+
+                polyValues = p_pmSubtractionPolynomialFromCoords(polyValues, kernels, numCols, numRows,
+                                                                 i + size + 1, j + size + 1);
+                int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues,
+                                                   false, poorFrac); // Radius of bad box
+                if (box > 0) {
+                    // Convolve a subimage, then stick it in the original
+                    psImage *subMask = psImageSubset(mask, psRegionSet(i - box, xSubMax + box,
+                                                                       j - box, ySubMax + box)); // Subimage
+                    psImage *convolved = psImageConvolveMask(NULL, subMask, 0x0f, 0xf0,
+                                                             -box, box, -box, box); // Convolved mask
+                    psFree(subMask);
+
+                    int numBytes = (xSubMax - i) * PSELEMTYPE_SIZEOF(PS_TYPE_MASK); // Number of bytes to copy
+                    psAssert(convolved->numCols - 2 * box == xSubMax - i, "Bad number of columns");
+                    psAssert(convolved->numRows - 2 * box == ySubMax - j, "Bad number of rows");
+
+                    for (int yTarget = j, ySource = box; yTarget < ySubMax; yTarget++, ySource++) {
+                        memcpy(&mask->data.PS_TYPE_MASK_DATA[yTarget][i],
+                               &convolved->data.PS_TYPE_MASK_DATA[ySource][box], numBytes);
+                    }
+                    psFree(convolved);
+                }
             }
         }
+        psFree(polyValues);
     }
     bad = psPixelsFromMask(bad, mask, 0xff);
-    psFree(mask);
-
-    // Convert coordinates to frame of original image
-    for (int i = 0; i < bad->n; i++) {
-        int x = bad->data[i].x + x0;
-        int y = bad->data[i].y + y0;
-        if (x < 0 || x >= numCols || y < 0 || y >= numRows) {
-            psWarning("Bad pixel coordinate %d: %d,%d --- ignored.",
-                      i, x, y);
-            continue;
-        }
-        bad->data[i].x = x;
-        bad->data[i].y = y;
-    }
 
     return bad;