Changeset 19165 for trunk/ppSub/src/ppSubReadout.c

Timestamp:

Aug 22, 2008, 12:46:26 PM (18 years ago)

Author:

Paul Price

Message:

Introducing bad/poor distinction to reduce explosion of bad pixels in output. Attempt to interpolate over bad pixels before subtracting.

File:

• trunk/ppSub/src/ppSubReadout.c (modified) (10 diffs)

trunk/ppSub/src/ppSubReadout.c

@@ -98 +98 @@
     int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
     float penalty = psMetadataLookupF32(NULL, recipe, "PENALTY"); // Penalty for wideness
-    psMaskType maskIn = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask going in to pmSubtractionMatch
-    psMaskType maskOut = pmConfigMaskGet("BLANK", config); // Bits to mask after pmSubtractionMatch
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
+    psMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
+    psMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
+    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
+    psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
     float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
     const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps
@@ -113 +117 @@
     int renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width for renormalise
 
-    pmSubtractionMode mode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtraction mode
+    psString interpModeStr = psMetadataLookupStr(&mdok, recipe, "INTERPOLATION"); // Interpolation mode
+    psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp
+    if (interpMode == PS_INTERPOLATE_NONE) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unknown interpolation mode: %s", interpModeStr);
+        return false;
+    }
+    float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
+
+    pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode
 
     // Generate masks if they don't exist
@@ -134 +146 @@
     }
 
-    if (!pmReadoutMaskNonfinite(inRO, maskIn)) {
+    if (!pmReadoutMaskNonfinite(inRO, pmConfigMaskGet("SAT", config))) {
         psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in input.");
         return false;
     }
-    if (!pmReadoutMaskNonfinite(refRO, maskIn)) {
+    if (!pmReadoutMaskNonfinite(refRO, pmConfigMaskGet("SAT", config))) {
         psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference.");
         return false;
@@ -153 +165 @@
     }
 
+    // Interpolate over bad pixels, so the bad pixels don't explode
+    if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image.");
+        return false;
+    }
+    if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image.");
+        return false;
+    }
+    maskVal |= maskBad;
+
+    // Match the PSFs
     if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, regionSize, spacing, threshold, sources,
                             stampsName, type, size, order, widths, orders, inner, ringsOrder, binning,
-                            penalty, optimum, optWidths, optOrder, optThresh, iter, rej, maskIn,
-                            maskOut, badFrac, mode)) {
+                            penalty, optimum, optWidths, optOrder, optThresh, iter, rej, maskVal,
+                            maskBad, maskPoor, poorFrac, badFrac, subMode)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
         psFree(inConv);
@@ -221 +245 @@
 #ifdef TESTING
     {
-        pmReadoutMaskApply(minuend, maskOut);
+        pmReadoutMaskApply(minuend, maskVal);
         psFits *fits = psFitsOpen("minuend.fits", "w");
         psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
@@ -227 +251 @@
     }
     {
-        pmReadoutMaskApply(subtrahend, maskOut);
+        pmReadoutMaskApply(subtrahend, maskVal);
         psFits *fits = psFitsOpen("subtrahend.fits", "w");
         psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
@@ -267 +291 @@
 
         // set maskValue and markValue in the psphot recipe
-        psMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
+        psMaskType maskValue = maskVal;
         psMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking
         psMetadataAddU8(recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue);
@@ -306 +330 @@
     outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
 
-    pmReadoutMaskApply(outRO, maskOut);
+    pmReadoutMaskApply(outRO, maskBad);
 
     psFree(inConv);
@@ -314 +338 @@
     for (int y = 0; y < outRO->image->numRows; y++) {
         for (int x = 0; x < outRO->image->numCols; x++) {
-            if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.U8[y][x] & maskOut)) {
+            if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.U8[y][x] & maskVal)) {
                 printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.U8[y][x]);
             }
@@ -329 +353 @@
     if (renorm) {
         psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
-        if (!pmReadoutWeightRenorm(outRO, maskOut, PS_STAT_ROBUST_MEDIAN, PS_STAT_ROBUST_STDEV,
+        if (!pmReadoutWeightRenorm(outRO, maskBad, PS_STAT_ROBUST_MEDIAN, PS_STAT_ROBUST_STDEV,
                                    renormWidth, rng)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to renormalise weights for photometry.");