Changeset 23937

Timestamp:

Apr 20, 2009, 2:48:42 PM (17 years ago)

Author:

Paul Price

Message:

Implement stamp thresholds on both images. This is to avoid transients being selected as stamps.

Location:

trunk/psModules/src/imcombine

Files:

• pmSubtractionMatch.c (modified) (5 diffs)
• pmSubtractionStamps.c (modified) (5 diffs)
• pmSubtractionStamps.h (modified) (1 diff)

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -64 +64 @@
                       psImage *variance,  // Variance map
                       const psRegion *region, // Region of interest, or NULL
-                      float threshold,  // Threshold for stamp finding
+                      float thresh1,  // Threshold for stamp finding on readout 1
+                      float thresh2,  // Threshold for stamp finding on readout 2
                       float stampSpacing, // Spacing between stamps
                       int size,         // Kernel half-size
@@ -72 +73 @@
 {
     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
-    *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, size, footprint,
-                                      stampSpacing, mode);
+
+    psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest
+
+    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
+                                      size, footprint, stampSpacing, mode);
     if (!*stamps) {
         psError(psErrorCodeLast(), false, "Unable to find stamps.");
@@ -392 +396 @@
     }
 
+    float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
     {
         psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun
-        if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
-            psFree(bg);
-            psFree(rng);
-            goto MATCH_ERROR;
-        }
-        threshold = bg->robustMedian + threshold * bg->robustStdev;
+        if (ro1) {
+            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
+                psFree(bg);
+                psFree(rng);
+                goto MATCH_ERROR;
+            }
+            stampThresh1 = bg->robustMedian + threshold * bg->robustStdev;
+        }
+        if (ro2) {
+            if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
+                psFree(bg);
+                psFree(rng);
+                goto MATCH_ERROR;
+            }
+            stampThresh2 = bg->robustMedian + threshold * bg->robustStdev;
+        }
         psFree(bg);
     }
@@ -428 +444 @@
             // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
             // doesn't matter.
-            if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, threshold, stampSpacing,
-                           size, footprint, subMode)) {
+            if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
+                           stampSpacing, size, footprint, subMode)) {
                 goto MATCH_ERROR;
             }
@@ -490 +506 @@
                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
 
-                if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, threshold, stampSpacing,
-                               size, footprint, subMode)) {
+                if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
+                               stampSpacing, size, footprint, subMode)) {
                     goto MATCH_ERROR;
                 }

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -232 +232 @@
 
 
-pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image,
-                                                const psImage *subMask, const psRegion *region,
-                                                float threshold, int size, int footprint, float spacing,
+pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image1,
+                                                const psImage *image2, const psImage *subMask,
+                                                const psRegion *region, float thresh1, float thresh2,
+                                                int size, int footprint, float spacing,
                                                 pmSubtractionMode mode)
 {
-    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
+    if (!image1 && !image2) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Must specify an image");
+        return NULL;
+    }
+    int numCols = 0, numRows = 0;       // Size of images
+    if (image1) {
+        PS_ASSERT_IMAGE_NON_NULL(image1, NULL);
+        PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, NULL);
+        if (subMask) {
+            PS_ASSERT_IMAGES_SIZE_EQUAL(image1, subMask, NULL);
+        }
+        numCols = image1->numCols;
+        numRows = image1->numRows;
+    }
+    if (image2) {
+        PS_ASSERT_IMAGE_NON_NULL(image2, NULL);
+        PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, NULL);
+        if (subMask) {
+            PS_ASSERT_IMAGES_SIZE_EQUAL(image2, subMask, NULL);
+        }
+        numCols = image2->numCols;
+        numRows = image2->numRows;
+    }
+    if (image1 && image2) {
+        PS_ASSERT_IMAGES_SIZE_EQUAL(image1, image2, NULL);
+    }
     if (subMask) {
         PS_ASSERT_IMAGE_NON_NULL(subMask, NULL);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL);
         PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_IMAGE_MASK, NULL);
+        PS_ASSERT_IMAGE_SIZE(subMask, numCols, numRows, NULL);
     }
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -254 +279 @@
             return false;
         }
-        if (region->x0 < 0 || region->x1 > image->numCols ||
-            region->y0 < 0 || region->y1 > image->numRows) {
+        if (region->x0 < 0 || region->x1 > numCols ||
+            region->y0 < 0 || region->y1 > numRows) {
             psString string = psRegionToString(*region);
             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)",
-                    string, image->numCols, image->numRows);
+                    string, numCols, numRows);
             psFree(string);
             return false;
@@ -264 +289 @@
     }
 
-    int numRows = image->numRows, numCols = image->numCols; // Size of image
     int border = size + footprint;      // Border size
 
@@ -312 +336 @@
                     fluxList->n = index;
 
-                    goodStamp = (fluxStamp > threshold) ? true : false;
+                    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);
@@ -319 +349 @@
                 // Use a simple method of automatically finding stamps --- take the highest pixel in the
                 // subregion
-                fluxStamp = threshold;
+                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 (image->data.F32[y][x] > fluxStamp &&
+                        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];

trunk/psModules/src/imcombine/pmSubtractionStamps.h

@@ -74 +74 @@
 /// Find stamps on an image
 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL
-                                                const psImage *image, ///< Image for which to find stamps
+                                                const psImage *image1, ///< Image for which to find stamps
+                                                const psImage *image2, ///< Image for which to find stamps
                                                 const psImage *mask, ///< Mask, or NULL
                                                 const psRegion *region, ///< Region to search, or NULL
-                                                float threshold, ///< Threshold for stamps in the image
+                                                float thresh1, ///< Threshold for stamps in image 1
+                                                float thresh2, ///< Threshold for stamps in image 2
                                                 int size, ///< Kernel half-size
                                                 int footprint, ///< Half-size for stamps