Changeset 24066

Timestamp:

May 4, 2009, 6:29:42 PM (17 years ago)

Author:

Paul Price

Message:

Search around provided stamp centre for a suitable stamp. If there's
a lot of masking going on, the centre of the stamp may not be
suitable, but there's no reason why we can't offset a bit.
Attempt to address ticket 1235.

File:

• trunk/psModules/src/imcombine/pmSubtractionStamps.c (modified) (5 diffs)

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -78 +78 @@
 }
 
-// Is this position unmasked?
-static bool checkStampMask(int x, int y, // Coordinates of stamp
-                           int numCols, int numRows, // Size of image
-                           const psImage *mask, // Mask
-                           pmSubtractionMode mode, // Mode for subtraction
-                           int footprint, // Size of stamp
-                           int border // Size of border
-                           )
-{
-    // Check the footprint bounds
-    if (x < border || x >= numCols - border || y < border || y >= numRows - border) {
-        return false;
-    }
-
-    if (!mask) {
-        return true;
-    }
-
-    // Check the central pixel
-    if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) {
-        mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
-        return false;
-    }
-
-    return true;
+
+// Search a region for a suitable stamp
+bool stampSearch(int *xStamp, int *yStamp, // Coordinates of stamp, to return
+                 float *fluxStamp, // Flux of stamp, to return
+                 const psImage *image1, const psImage *image2, // Images to search
+                 float thresh1, float thresh2, // Thresholds for images
+                 const psImage *subMask, // Subtraction mask
+                 int xMin, int xMax, int yMin, int yMax, // Bounds of search
+                 int numCols, int numRows, // Size of images
+                 int border             // Border around image
+    )
+{
+    bool found = false;                 // Found a suitable stamp?
+    *fluxStamp = -INFINITY;             // Flux of best stamp
+
+    // Ensure we're not going to go outside the bounds of the image
+    xMin = PS_MAX(border, xMin);
+    xMax = PS_MIN(numCols - border - 1, xMax);
+    yMin = PS_MAX(border, yMin);
+    yMax = PS_MIN(numRows - border - 1, yMax);
+
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            if ((image1 && image1->data.F32[y][x] < thresh1) ||
+                (image2 && image2->data.F32[y][x] < thresh2)) {
+                continue;
+            }
+
+            if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] &
+                (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) {
+                return false;
+            }
+
+            // We take the MIN to attempt to avoid transients in both images
+            float flux = (image1 && image2) ? PS_MIN(image1->data.F32[y][x], image2->data.F32[y][x]) :
+                ((image1) ? image1->data.F32[y][x] : image2->data.F32[y][x]); // Flux at pixel
+            if (flux > *fluxStamp) {
+                *xStamp = x;
+                *yStamp = y;
+                *fluxStamp = flux;
+                found = true;
+            }
+        }
+    }
+
+    return found;
 }
 
@@ -314 +335 @@
             numSearch++;
 
-            float xStamp = 0, yStamp = 0;   // Coordinates of stamp
-            float fluxStamp = NAN;          // Flux of stamp
-            bool goodStamp = false;         // Found a good stamp?
+            int xStamp = 0, yStamp = 0; // Coordinates of stamp
+            float fluxStamp = -INFINITY; // Flux of stamp
+            bool goodStamp = false;     // Found a good stamp?
 
             // A couple different ways of finding stamps:
@@ -324 +345 @@
                 psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
 
-                // Take stamp off the top of the (sorted) list
-                if (xList->n > 0) {
-                    int index = xList->n - 1; // Index of new stamp
-                    xStamp = xList->data.F32[index];
-                    yStamp = yList->data.F32[index];
-                    fluxStamp = fluxList->data.F32[index];
+                // Take stamps off the top of the (sorted) list
+                for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) {
+                    int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre
 
                     // Chop off the top of the list
-                    xList->n = index;
-                    yList->n = index;
-                    fluxList->n = index;
-
-                    goodStamp = true;
-                    int x = xStamp + 0.5, y = yStamp + 0.5; // Pixel location
-                    if (image1 && image1->data.F32[y][x] < thresh1) {
-                        goodStamp = false;
-                    } else if (image2 && image2->data.F32[y][x] < thresh2) {
-                        goodStamp = false;
-                    }
-                } else {
-                    psTrace("psModules.imcombine", 9, "No sources in subregion %d", i);
+                    xList->n = j;
+                    yList->n = j;
+                    fluxList->n = j;
+
+                    // Fish around a bit to see if we can find a pixel that isn't masked
+                    psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n",
+                            i, xCentre, yCentre);
+
+                    // Search bounds
+                    int search = footprint - size; // Search radius
+                    int xMin = PS_MAX(border, xCentre - search);
+                    int xMax = PS_MIN(numCols - border -1, xCentre + search);
+                    int yMin = PS_MAX(border, yCentre - search);
+                    int yMax = PS_MIN(numRows - border - 1, yCentre + search);
+
+                    goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
+                                            subMask, xMin, xMax, yMin, yMax, numCols, numRows, border);
                 }
             } else {
                 // Use a simple method of automatically finding stamps --- take the highest pixel in the
                 // subregion
-                fluxStamp = 0.0;
                 psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest
 
-                const psImage *image = NULL; // Image for flux of stamp
-                switch (mode) {
-                  case PM_SUBTRACTION_MODE_1:
-                    image = image2 ? image2 : image1;
-                    break;
-                  case PM_SUBTRACTION_MODE_2:
-                    image = image1 ? image1 : image2;
-                    break;
-                  case PM_SUBTRACTION_MODE_UNSURE:
-                  case PM_SUBTRACTION_MODE_DUAL:
-                    // If both are available, we'll take the MIN value below in the hope of rejecting
-                    // transients in both
-                    image = (image1 && image2) ? NULL : (image1 ? image1 : image2);
-                  default:
-                    psAbort("Unrecognised mode: %x", mode);
-                }
-
-                for (int y = subRegion->y0; y <= subRegion->y1; y++) {
-                    for (int x = subRegion->x0; x <= subRegion->x1; x++) {
-                        if ((image1 && image1->data.F32[y][x] < thresh1) ||
-                            (image2 && image2->data.F32[y][x] < thresh2)) {
-                            continue;
-                        }
-                        float value = image ? image->data.F32[y][x] :
-                            PS_MIN(image1->data.F32[y][x], image2->data.F32[y][x]); // Value of image
-                        if (value > fluxStamp &&
-                            checkStampMask(x, y, numCols, numRows, subMask, mode, footprint, border)) {
-                            fluxStamp = image->data.F32[y][x];
-                            xStamp = x;
-                            yStamp = y;
-                            goodStamp = true;
-                        }
-                    }
-                }
+                goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
+                                        subMask, subRegion->x0, subRegion->x1, subRegion->y0, subRegion->y1,
+                                        numCols, numRows, border);
             }
 
@@ -453 +443 @@
                                                                  region, footprint, spacing); // Stamp list
     int numStamps = stamps->num;        // Number of stamps
-    int numCols = image->numCols, numRows = image->numRows; // Size of image
-    int border = footprint + size;      // Border size
 
     psString ds9name = NULL;            // Filename for ds9 region file
@@ -480 +468 @@
             psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because outside region",
                     xPix, yPix);
-            pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "magenta");
-            continue;
-        }
-        if (!checkStampMask(xPix, yPix, numCols, numRows, subMask, mode, footprint, border)) {
-            // Not a good stamp
-            psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask",
-                    xPix, yPix);
             pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "red");
             continue;