Changeset 26739 for branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Jan 29, 2010, 6:02:38 PM (16 years ago)

Author:

Paul Price

Message:

Fairly extensive change to make subtractions more stable on images with limited lit area in the presence of spatial variation. Took the polynomials which were scaled over the entire image and scale them only over the lit area. Removed some fprintf statements. Fixed the problem causing the background term to always come out 1.0.

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (modified) (15 diffs)

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c

@@ -70 +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
@@ -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");
 
@@ -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;
@@ -319 +326 @@
     pmSubtractionMaskInvalid(ro2, maskVal);
 
-    psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
+    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) {
@@ -441 +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
@@ -457 +466 @@
     pmSubtractionMaskInvalid(ro2, maskVal);
 
-    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
-                                badFrac, subMode);
+    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.");
@@ -468 +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;
     }
 
@@ -480 +492 @@
     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);
@@ -504 +516 @@
     }
 
+    // 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_1 || 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
     for (int j = 0; j < yRegions; j++) {
@@ -509 +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) {
@@ -530 +583 @@
             // 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,
+            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
                                       stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
                 goto MATCH_ERROR;
@@ -544 +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();
@@ -554 +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);
             }
@@ -603 +657 @@
                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
 
-#if 0
                 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
                                           stampThresh1, stampThresh2, stampSpacing, normFrac,
@@ -609 +662 @@
                     goto MATCH_ERROR;
                 }
-#endif
 
                 // generate the window function from the set of stamps