Changeset 26747 for branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.c

Timestamp:

Jan 31, 2010, 5:00:42 PM (16 years ago)

Author:

eugene

Message:

updates relative to 20091201, fixes for all psphot variants

Location:

branches/eam_branches/psModules.stack.20100120

Files:

• . (modified) (1 prop)
• src/imcombine/pmSubtractionMatch.c (modified) (22 diffs)

branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.c

@@ -33 +33 @@
 # else
 # define SUBMODE PM_SUBTRACTION_EQUATION_ALL
-// # define SUBMODE PM_SUBTRACTION_EQUATION_KERNELS
 # endif
 
@@ -71 +70 @@
                                  const psImage *subMask, // Mask for subtraction, or NULL
                                  psImage *variance,  // Variance map
-                                 const psRegion *region, // Region of interest, or NULL
+                                 const psRegion *region, // Region of interest
                                  float thresh1,  // Threshold for stamp finding on readout 1
                                  float thresh2,  // Threshold for stamp finding on readout 2
                                  float stampSpacing, // Spacing between stamps
+                                 float normFrac,     // Fraction of flux in window for normalisation window
                                  float sysError,     // Relative systematic error in images
                                  float skyError,     // Relative systematic error in images
@@ -82 +82 @@
     )
 {
+    PS_ASSERT_PTR_NON_NULL(stamps, false);
+    PM_ASSERT_READOUT_NON_NULL(ro1, false);
+    PM_ASSERT_READOUT_NON_NULL(ro2, false);
+    PS_ASSERT_IMAGE_NON_NULL(subMask, false);
+    PS_ASSERT_IMAGE_NON_NULL(variance, false);
+    PS_ASSERT_PTR_NON_NULL(region, false);
+
     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
 
@@ -87 +94 @@
 
     *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
-                                      size, footprint, stampSpacing, sysError, skyError, mode);
+                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
     if (!*stamps) {
         psError(psErrorCodeLast(), false, "Unable to find stamps.");
@@ -96 +103 @@
 
     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) {
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
         return false;
@@ -110 +117 @@
                                   const pmReadout *ro1, const pmReadout *ro2, // Input images
                                   int stride, // Size for convolution patches
+                                  float normFrac,           // Fraction of window for normalisation window
                                   float sysError,           // Systematic error in images
                                   float skyError,           // Systematic error in images
@@ -158 +166 @@
     PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
     PS_ASSERT_INT_NONNEGATIVE(stride, false);
+    if (isfinite(normFrac)) {
+        PS_ASSERT_FLOAT_LARGER_THAN(normFrac, 0.0, false);
+        PS_ASSERT_FLOAT_LESS_THAN(normFrac, 1.0, false);
+    }
     if (isfinite(sysError)) {
         PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false);
@@ -210 +222 @@
 }
 
+bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) {
+
+    if (!readout) return true;
+
+    psImage *image = readout->image;
+    psImage *mask  = readout->mask;
+    psImage *variance = readout->variance;
+    for (int y = 0; y < image->numRows; y++) {
+        for (int x = 0; x < image->numCols; x++) {
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue;
+            bool valid = false;
+            valid = isfinite(image->data.F32[y][x]);
+            if (variance) {
+                valid &= isfinite(variance->data.F32[y][x]);
+            }
+            if (valid) continue;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal;
+        }
+    }
+
+    return true;
+}
 
 bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
@@ -282 +316 @@
     }
 
-    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, kernelError,
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, NAN, kernelError,
                                maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) {
         psFree(kernels);
@@ -289 +323 @@
     }
 
-    psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
+    pmSubtractionMaskInvalid(ro1, maskVal);
+    pmSubtractionMaskInvalid(ro2, maskVal);
+
+    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
+
+    psImage *subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, 0,
                                          badFrac, mode); // Subtraction mask
     if (!subMask) {
@@ -348 +387 @@
                         int inner, int ringsOrder, int binning, float penalty,
                         bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
-                        int iter, float rej, float sysError, float skyError, float kernelError, psImageMaskType maskVal,
-                        psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac,
-                        float badFrac, pmSubtractionMode subMode)
+                        int iter, float rej, float normFrac, float sysError, float skyError,
+                        float kernelError, psImageMaskType maskVal, psImageMaskType maskBad,
+                        psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode)
 {
-    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, skyError, kernelError,
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, normFrac, sysError, skyError, kernelError,
                                maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
         return false;
@@ -411 +450 @@
     // Putting important variable declarations here, since they are freed after a "goto" if there is an error.
     psImage *subMask = NULL;            // Mask for subtraction
-    psRegion *region = NULL;            // Iso-kernel region
+    psRegion *region = psRegionAlloc(NAN, NAN, NAN, NAN); // Iso-kernel region
     psString regionString = NULL;       // String for region
     pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF
@@ -424 +463 @@
     memCheck("start");
 
-    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
-                                badFrac, subMode);
+    pmSubtractionMaskInvalid(ro1, maskVal);
+    pmSubtractionMaskInvalid(ro2, maskVal);
+
+    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
+
+    subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, footprint, badFrac, subMode);
     if (!subMask) {
         psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
@@ -435 +478 @@
     // Get region of interest
     int xRegions = 1, yRegions = 1;     // Number of iso-kernel regions
-    float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions
+    float xRegionSize = NAN, yRegionSize = NAN; // Size of iso-kernel regions
     if (isfinite(regionSize) && regionSize != 0.0) {
-        xRegions = numCols / regionSize + 1;
-        yRegions = numRows / regionSize + 1;
-        xRegionSize = (float)numCols / (float)xRegions;
-        yRegionSize = (float)numRows / (float)yRegions;
-        region = psRegionAlloc(NAN, NAN, NAN, NAN);
-    }
-
+        xRegions = (bounds.x1 - bounds.x0) / regionSize + 1;
+        yRegions = (bounds.y1 - bounds.y0) / regionSize + 1;
+        xRegionSize = (float)(bounds.x1 - bounds.x0) / (float)xRegions;
+        yRegionSize = (float)(bounds.y1 - bounds.y0) / (float)yRegions;
+    } else {
+        xRegionSize = bounds.x1 - bounds.x0;
+        yRegionSize = bounds.y1 - bounds.y0;
+    }
+
+    // General background subtraction and measurement of stamp threshold
     float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
     {
-        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun
+        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
         if (ro1) {
+            psStatsInit(bg);
             if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
@@ -453 +500 @@
                 goto MATCH_ERROR;
             }
-            stampThresh1 = bg->robustMedian + threshold * bg->robustStdev;
+            stampThresh1 = threshold * bg->robustStdev;
+            psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
         }
         if (ro2) {
+            psStatsInit(bg);
             if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
@@ -461 +510 @@
                 goto MATCH_ERROR;
             }
-            stampThresh2 = bg->robustMedian + threshold * bg->robustStdev;
-        }
+            stampThresh2 = threshold * bg->robustStdev;
+            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
+       }
         psFree(bg);
     }
+
+    // Just in case the iso-kernel region doesn't cover the entire image...
+    if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
+        subMode == PM_SUBTRACTION_MODE_DUAL) {
+        if (!conv1->image) {
+            conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        psImageInit(conv1->image, NAN);
+        if (ro1->variance) {
+            if (!conv1->variance) {
+                conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            }
+            psImageInit(conv1->variance, NAN);
+        }
+        if (subMask) {
+            if (!conv1->mask) {
+                conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+            }
+            psImageInit(conv1->mask, maskBad);
+        }
+    }
+    if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
+        subMode == PM_SUBTRACTION_MODE_DUAL) {
+        if (!conv2->image) {
+            conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        psImageInit(conv2->image, NAN);
+        if (ro2->variance) {
+            if (!conv2->variance) {
+                conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            }
+            psImageInit(conv2->variance, NAN);
+        }
+        if (subMask) {
+            if (!conv2->mask) {
+                conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+            }
+            psImageInit(conv2->mask, maskBad);
+        }
+    }
+
 
     // Iterate over iso-kernel regions
@@ -471 +562 @@
             psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n",
                     j * xRegions + i + 1, xRegions * yRegions);
-            if (region) {
-                *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize),
-                                      (int)(j * yRegionSize), (int)((j + 1) * yRegionSize));
-                psFree(regionString);
-                regionString = psRegionToString(*region);
-                psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n",
-                        regionString, numCols, numRows);
-            }
+            *region = psRegionSet(bounds.x0 + (int)(i * xRegionSize),
+                                  bounds.x0 + (int)((i + 1) * xRegionSize),
+                                  bounds.y0 + (int)(j * yRegionSize),
+                                  bounds.y0 + (int)((j + 1) * yRegionSize));
+            psFree(regionString);
+            regionString = psRegionToString(*region);
+            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n",
+                    regionString, numCols, numRows);
 
             if (stampsName && strlen(stampsName) > 0) {
                 stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
-                                                        footprint, stampSpacing, sysError, skyError, subMode);
+                                                        footprint, stampSpacing, normFrac,
+                                                        sysError, skyError, subMode);
             } else if (sources) {
                 stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
-                                                           footprint, stampSpacing, sysError, skyError, subMode);
+                                                           footprint, stampSpacing, normFrac,
+                                                           sysError, skyError, subMode);
             }
 
             // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
             // doesn't matter.
-            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
-                                      stampSpacing, sysError, skyError, size, footprint, subMode)) {
+            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
+                                      stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
                 goto MATCH_ERROR;
             }
@@ -498 +591 @@
             // generate the window function from the set of stamps
             if (!pmSubtractionStampsGetWindow(stamps, size)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to get stamp window.");
                 goto MATCH_ERROR;
             }
@@ -503 +597 @@
             // Define kernel basis functions
             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
-                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
-                                                          stamps, footprint, optThreshold, penalty, subMode);
+                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
+                                                          optFWHMs, optOrder, stamps, footprint,
+                                                          optThreshold, penalty, bounds, subMode);
                 if (!kernels) {
                     psErrorClear();
@@ -513 +608 @@
                 // Not an ISIS/GUNK kernel, or the optimum kernel search failed
                 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
-                                                       inner, binning, ringsOrder, penalty, subMode);
+                                                       inner, binning, ringsOrder, penalty, bounds, subMode);
                 // pmSubtractionVisualShowKernels(kernels);
             }
@@ -563 +658 @@
 
                 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
-                                          stampThresh1, stampThresh2, stampSpacing, sysError, skyError,
-                                          size, footprint, subMode)) {
+                                          stampThresh1, stampThresh2, stampSpacing, normFrac,
+                                          sysError, skyError, size, footprint, subMode)) {
                     goto MATCH_ERROR;
                 }
@@ -570 +665 @@
                 // generate the window function from the set of stamps
                 if (!pmSubtractionStampsGetWindow(stamps, size)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to get stamps window.");
                     goto MATCH_ERROR;
                 }