Changeset 26035 for trunk/psModules/src/imcombine/pmSubtraction.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/pmSubtraction.c (modified) (17 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -29 +29 @@
 #define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
 #define MIN_SAMPLE_STATS    7           // Minimum number to use sample statistics; otherwise use quartiles
-#define USE_SYS_ERR                   // Use systematic error image?
+#define USE_KERNEL_ERR                  // Use kernel error image?
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -266 +266 @@
 static void convolveVarianceFFT(psImage *target,// Place the result in here
                               psImage *variance, // Variance map to convolve
-                              psImage *sys, // Systematic error image
+                              psImage *kernelErr, // Kernel error image
                               psImage *mask, // Mask image
                               psImageMaskType maskVal, // Value to mask
@@ -278 +278 @@
 
     psImage *subVariance = variance ? psImageSubset(variance, border) : NULL; // Variance map
-    psImage *subSys = sys ? psImageSubset(sys, border) : NULL; // Systematic error image
+    psImage *subKE = kernelErr ? psImageSubset(kernelErr, border) : NULL; // Kernel error image
     psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Mask
 
     // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
     psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance
-    psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys
+    psImage *convKE = subKE ? psImageConvolveFFT(NULL, subKE, subMask, maskVal, kernel) : NULL; // Conv KE
 
     psFree(subVariance);
-    psFree(subSys);
+    psFree(subKE);
     psFree(subMask);
 
     // Now, we have to stick it in where it belongs
     int xMin = region.x0, xMax = region.x1, yMin = region.y0, yMax = region.y1; // Bounds of region
-    if (convSys) {
+    if (convKE) {
         for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
             for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) {
                 target->data.F32[yTarget][xTarget] = convVariance->data.F32[ySource][xSource] +
-                    convSys->data.F32[ySource][xSource];
+                    convKE->data.F32[ySource][xSource];
             }
         }
@@ -306 +306 @@
 
     psFree(convVariance);
-    psFree(convSys);
+    psFree(convKE);
 
     return;
@@ -342 +342 @@
                                   psImage *image, // Image to convolve
                                   psImage *variance, // Variance map to convolve, or NULL
-                                  psImage *sys, // Systematic error image, or NULL
+                                  psImage *kernelErr, // Kernel error image, or NULL
                                   psImage *subMask, // Subtraction mask
                                   const pmSubtractionKernels *kernels, // Kernels
@@ -379 +379 @@
         convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size);
         if (variance) {
-            convolveVarianceFFT(convVariance, variance, sys, subMask, subBad, *kernelVariance, region, kernels->size);
+            convolveVarianceFFT(convVariance, variance, kernelErr, subMask, subBad, *kernelVariance,
+                                region, kernels->size);
         }
     } else {
@@ -429 +430 @@
 }
 
-#ifdef USE_SYS_ERR
-// Generate an image that can be used to track systematic errors
-static psImage *subtractionSysErrImage(const psImage *image, // Image from which to make sys err image
-                                       float sysError // Relative systematic error
-                                       )
-{
-    if (!isfinite(sysError) || sysError == 0.0) {
+#ifdef USE_KERNEL_ERR
+// Generate an image that can be used to track systematic errors in the kernel
+static psImage *subtractionKernelErrImage(const psImage *image, // Image from which to make kernel error image
+                                          float kernelError // Relative systematic error in kernel
+    )
+{
+    if (!isfinite(kernelError) || kernelError == 0.0) {
         return NULL;
     }
 
     int numCols = image->numCols, numRows = image->numRows; // Size of image
-    psImage *sys = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Systematic error image
-
-    float sysError2 = PS_SQR(sysError); // Square of the systematic error
+    psImage *kernelErr = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Kernel error image
+
+    float kernelError2 = PS_SQR(kernelError); // Square of the kernel error
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            sys->data.F32[y][x] = PS_SQR(image->data.F32[y][x]) * sysError2;
-        }
-    }
-
-    return sys;
+            kernelErr->data.F32[y][x] = PS_SQR(image->data.F32[y][x]) * kernelError2;
+        }
+    }
+
+    return kernelErr;
 }
 #endif
@@ -892 +893 @@
                 psFree(stamp->image1);
                 psFree(stamp->image2);
-                psFree(stamp->variance);
-                stamp->image1 = stamp->image2 = stamp->variance = NULL;
+                psFree(stamp->weight);
+                stamp->image1 = stamp->image2 = stamp->weight = NULL;
                 psFree(stamp->matrix1);
                 psFree(stamp->matrix2);
@@ -1040 +1041 @@
                                      psImage *convMask, // Output convolved mask
                                      const pmReadout *ro1, const pmReadout *ro2, // Input readouts
-                                     psImage *sys1, psImage *sys2, // Systematic error images
+                                     psImage *kernelErr1, psImage *kernelErr2, // Kernel error images
                                      psImage *subMask, // Input subtraction mask
                                      psImageMaskType maskBad, // Mask value to give bad pixels
@@ -1066 +1067 @@
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         convolveRegion(out1->image, out1->variance, convMask, &kernelImage, &kernelVariance,
-                       ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region,
-                       maskBad, maskPoor, poorFrac, useFFT, false);
+                       ro1->image, ro1->variance, kernelErr1, subMask, kernels, polyValues, background,
+                       *region, maskBad, maskPoor, poorFrac, useFFT, false);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         convolveRegion(out2->image, out2->variance, convMask, &kernelImage, &kernelVariance,
-                       ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region,
-                       maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL);
+                       ro2->image, ro2->variance, kernelErr2, subMask, kernels, polyValues, background,
+                       *region, maskBad, maskPoor, poorFrac, useFFT,
+                       kernels->mode == PM_SUBTRACTION_MODE_DUAL);
     }
 
@@ -1117 +1119 @@
     const pmReadout *ro1 = args->data[7]; // Input readout 1
     const pmReadout *ro2 = args->data[8]; // Input readout 2
-    psImage *sys1 = args->data[9]; // Systematic error image 1
-    psImage *sys2 = args->data[10]; // Systematic error image 2
+    psImage *kernelErr1 = args->data[9]; // Kernel error image 1
+    psImage *kernelErr2 = args->data[10]; // Kernel error image 2
     psImage *subMask = args->data[11]; // Subtraction mask
     psImageMaskType maskBad = PS_SCALAR_VALUE(args->data[12], PS_TYPE_IMAGE_MASK_DATA); // Output mask value for bad pixels
@@ -1128 +1130 @@
     bool useFFT = PS_SCALAR_VALUE(args->data[18], U8); // Use FFT for convolution?
 
-    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
-                                    subMask, maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT);
+    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, kernelErr1,
+                                    kernelErr2, subMask, maskBad, maskPoor, poorFrac, region, kernels,
+                                    doBG, useFFT);
 }
 
 bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2,
                            psImage *subMask, int stride, psImageMaskType maskBad, psImageMaskType maskPoor,
-                           float poorFrac, float sysError, const psRegion *region,
+                           float poorFrac, float kernelError, const psRegion *region,
                            const pmSubtractionKernels *kernels, bool doBG, bool useFFT)
 {
@@ -1172 +1175 @@
     PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(poorFrac, 0.0, false);
     PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(poorFrac, 1.0, false);
-    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false);
-    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(sysError, 1.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false);
+    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(kernelError, 1.0, false);
     if (region && psRegionIsNaN(*region)) {
         psString string = psRegionToString(*region);
@@ -1232 +1235 @@
     }
 
-    psImage *sys1 = NULL, *sys2 = NULL; // Systematic error images
-#ifdef USE_SYS_ERR
+    psImage *kernelErr1 = NULL, *kernelErr2 = NULL; // Kernel error images
+#ifdef USE_KERNEL_ERR
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        sys1 = subtractionSysErrImage(ro1->image, sysError);
+        kernelErr1 = subtractionKernelErrImage(ro1->image, kernelError);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        sys2 = subtractionSysErrImage(ro2->image, sysError);
+        kernelErr2 = subtractionKernelErrImage(ro2->image, kernelError);
     }
 #endif
@@ -1289 +1292 @@
                 psArrayAdd(args, 1, (pmReadout*)ro1); // Casting away const
                 psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const
-                psArrayAdd(args, 1, sys1);
-                psArrayAdd(args, 1, sys2);
+                psArrayAdd(args, 1, kernelErr1);
+                psArrayAdd(args, 1, kernelErr2);
                 psArrayAdd(args, 1, subMask);
                 PS_ARRAY_ADD_SCALAR(args, maskBad, PS_TYPE_IMAGE_MASK);
@@ -1306 +1309 @@
                 psFree(job);
             } else {
-                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
-                                         subMask, maskBad, maskPoor, poorFrac, subRegion, kernels, doBG,
-                                         useFFT);
+                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2,
+                                         kernelErr1, kernelErr2, subMask, maskBad, maskPoor, poorFrac,
+                                         subRegion, kernels, doBG, useFFT);
             }
             psFree(subRegion);
@@ -1333 +1336 @@
     psImageConvolveSetThreads(oldThreads);
 
-    psFree(sys1);
-    psFree(sys2);
+    psFree(kernelErr1);
+    psFree(kernelErr2);
 
     // Calculate covariances