Changeset 26035 for trunk/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Nov 4, 2009, 4:52:05 PM (17 years ago)

Author:

Paul Price

Message:

Soften errors, both in the least-squares matrix/vector construction and in the rejection stage. Softening in the rejection stage seems to be very important so that bright stamps aren't rejected (after which the solution can go crazy). Now carrying around a weight image = 1/variance instead of the variance image. This can be softened, saving many divisions in the matrix/vector construction.

File:

• trunk/psModules/src/imcombine/pmSubtractionMatch.c (modified) (12 diffs)

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -68 +68 @@
                                  float thresh2,  // Threshold for stamp finding on readout 2
                                  float stampSpacing, // Spacing between stamps
+                                 float sysError,     // Relative systematic error in images
                                  int size,         // Kernel half-size
                                  int footprint,     // Convolution footprint for stamps
@@ -78 +79 @@
 
     *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
-                                      size, footprint, stampSpacing, mode);
+                                      size, footprint, stampSpacing, sysError, mode);
     if (!*stamps) {
         psError(psErrorCodeLast(), false, "Unable to find stamps.");
@@ -101 +102 @@
                                   const pmReadout *ro1, const pmReadout *ro2, // Input images
                                   int stride, // Size for convolution patches
-                                  float sysError, // Relative systematic error
+                                  float sysError,           // Systematic error in images
+                                  float kernelError, // Systematic error in kernel
                                   psImageMaskType maskVal, // Value to mask for input
                                   psImageMaskType maskBad, // Mask for output bad pixels
@@ -151 +153 @@
         PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false);
     }
+    if (isfinite(kernelError)) {
+        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false);
+        PS_ASSERT_FLOAT_LESS_THAN(kernelError, 1.0, false);
+    }
     // Don't care about maskVal
     // Don't care about maskBad
@@ -194 +200 @@
 
 bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
-                               psMetadata *analysis, int stride, float sysError,
+                               psMetadata *analysis, int stride, float kernelError,
                                psImageMaskType maskVal, psImageMaskType maskBad, psImageMaskType maskPoor,
                                float poorFrac, float badFrac)
@@ -264 +270 @@
     }
 
-    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor,
-                               poorFrac, badFrac, mode)) {
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, kernelError,
+                               maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) {
         psFree(kernels);
         psFree(regions);
@@ -300 +306 @@
 
         if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
-                                   sysError, region, kernel, true, useFFT)) {
+                                   kernelError, region, kernel, true, useFFT)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
             psFree(outAnalysis);
@@ -330 +336 @@
                         int inner, int ringsOrder, int binning, float penalty,
                         bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
-                        int iter, float rej, float sysError, psImageMaskType maskVal, psImageMaskType maskBad,
-                        psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode)
-{
-    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor,
-                               poorFrac, badFrac, subMode)) {
+                        int iter, float rej, float sysError, float kernelError, psImageMaskType maskVal,
+                        psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac,
+                        float badFrac, pmSubtractionMode subMode)
+{
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, kernelError,
+                               maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
         return false;
     }
@@ -463 +470 @@
             if (stampsName && strlen(stampsName) > 0) {
                 stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
-                                                        footprint, stampSpacing, subMode);
+                                                        footprint, stampSpacing, sysError, subMode);
             } else if (sources) {
                 stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
-                                                           footprint, stampSpacing, subMode);
+                                                           footprint, stampSpacing, sysError, subMode);
             }
 
@@ -472 +479 @@
             // doesn't matter.
             if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
-                           stampSpacing, size, footprint, subMode)) {
+                                      stampSpacing, sysError, size, footprint, subMode)) {
                 goto MATCH_ERROR;
             }
@@ -538 +545 @@
 
                 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
-                                          stampThresh1, stampThresh2, stampSpacing,
+                                          stampThresh1, stampThresh2, stampSpacing, sysError,
                                           size, footprint, subMode)) {
                     goto MATCH_ERROR;
@@ -608 +615 @@
             psTrace("psModules.imcombine", 2, "Convolving...\n");
             if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
-                                       sysError, region, kernels, true, useFFT)) {
+                                       kernelError, region, kernels, true, useFFT)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
                 goto MATCH_ERROR;