Changeset 27840 for branches/simtest_nebulous_branches/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

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

branches/simtest_nebulous_branches/psModules/src/imcombine/pmSubtractionMatch.c

@@ -28 +28 @@
 static bool useFFT = true;              // Do convolutions using FFT
 
+# define SEPARATE 0
+# if (SEPARATE)
+# define SUBMODE PM_SUBTRACTION_EQUATION_NORM
+# else
+# define SUBMODE PM_SUBTRACTION_EQUATION_ALL
+# endif
 
 //#define TESTING
@@ -64 +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
                                  int size,         // Kernel half-size
                                  int footprint,     // Convolution footprint for stamps
@@ -73 +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");
 
@@ -78 +94 @@
 
     *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
-                                      size, footprint, stampSpacing, mode);
+                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
     if (!*stamps) {
         psError(psErrorCodeLast(), false, "Unable to find stamps.");
@@ -87 +103 @@
 
     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
+        psError(psErrorCodeLast(), false, "Unable to extract stamps.");
         return false;
     }
@@ -101 +117 @@
                                   const pmReadout *ro1, const pmReadout *ro2, // Input images
                                   int stride, // Size for convolution patches
-                                  float sysError, // Relative systematic error
+                                  float normFrac,           // Fraction of window for normalisation window
+                                  float sysError,           // Systematic error in images
+                                  float skyError,           // Systematic error in images
+                                  float kernelError, // Systematic error in kernel
+                                  float covarFrac,   // Fraction for kernel truncation before covariance
                                   psImageMaskType maskVal, // Value to mask for input
                                   psImageMaskType maskBad, // Mask for output bad pixels
@@ -112 +132 @@
     if (subMode != PM_SUBTRACTION_MODE_2) {
         PM_ASSERT_READOUT_NON_NULL(conv1, false);
+        PM_ASSERT_READOUT_NON_NULL(ro1, false);
+        PM_ASSERT_READOUT_IMAGE(ro1, false);
         if (conv1->image) {
             psFree(conv1->image);
@@ -127 +149 @@
     if (subMode != PM_SUBTRACTION_MODE_1) {
         PM_ASSERT_READOUT_NON_NULL(conv2, false);
+        PM_ASSERT_READOUT_NON_NULL(ro2, false);
+        PM_ASSERT_READOUT_IMAGE(ro2, false);
         if (conv2->image) {
             psFree(conv2->image);
@@ -141 +165 @@
     }
 
-    PM_ASSERT_READOUT_NON_NULL(ro1, false);
-    PM_ASSERT_READOUT_NON_NULL(ro2, false);
-    PM_ASSERT_READOUT_IMAGE(ro1, false);
-    PM_ASSERT_READOUT_IMAGE(ro2, false);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
+    if (ro1 && ro2) {
+        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);
         PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false);
     }
+    if (isfinite(sysError)) {
+        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(skyError, 0.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);
+    }
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(covarFrac, 0.0, false);
+    PS_ASSERT_FLOAT_LESS_THAN(covarFrac, 1.0, false);
     // Don't care about maskVal
     // Don't care about maskBad
@@ -164 +199 @@
 }
 
+
+/// Allocate images, as required
+static void subtractionMatchAlloc(pmReadout *conv1, pmReadout *conv2, // Output readouts
+                                  const pmReadout *ro1, const pmReadout *ro2, // Input readouts
+                                  const psImage *subMask,                     // Subtraction mask
+                                  psImageMaskType maskBad,                    // Mask value for bad pixels
+                                  pmSubtractionMode subMode,          // Subtraction mode
+                                  int numCols, int numRows            // Size of 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);
+        }
+    }
+
+    return;
+}
+
+
 static void subtractionAnalysisUpdate(pmReadout *conv1, pmReadout *conv2, // Convolved images
                                       const psMetadata *analysis, // Analysis metadata
@@ -192 +280 @@
 }
 
+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,
-                               psMetadata *analysis, int stride, float sysError,
+                               psMetadata *analysis, int stride, float kernelError, float covarFrac,
                                psImageMaskType maskVal, psImageMaskType maskBad, psImageMaskType maskPoor,
                                float poorFrac, float badFrac)
@@ -210 +320 @@
         while ((item = psMetadataGetAndIncrement(iter))) {
             if (item->type != PS_DATA_UNKNOWN) {
-                psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unexpected type for kernel.");
+                psError(PM_ERR_PROG, true, "Unexpected type for kernel.");
                 psFree(iter);
                 psFree(kernelList);
@@ -230 +340 @@
     }
     if (psListLength(kernelList) == 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find kernels");
+        psError(PM_ERR_PROG, true, "Unable to find kernels");
         psFree(kernelList);
         return false;
@@ -245 +355 @@
         while ((item = psMetadataGetAndIncrement(iter))) {
             if (item->type != PS_DATA_REGION) {
-                psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unexpected type for region.");
+                psError(PM_ERR_PROG, true, "Unexpected type for region.");
                 psFree(iter);
                 psFree(kernels);
@@ -257 +367 @@
     }
     if (regions->n != kernels->n) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Differing number of kernels (%ld) and regions (%ld)",
+        psError(PM_ERR_PROG, true, "Differing number of kernels (%ld) and regions (%ld)",
                 kernels->n, regions->n);
         psFree(regions);
@@ -264 +374 @@
     }
 
-    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor,
-                               poorFrac, badFrac, mode)) {
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, NAN, kernelError, covarFrac,
+                               maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) {
         psFree(kernels);
         psFree(regions);
@@ -271 +381 @@
     }
 
-    psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
+    int numCols, numRows;       // Size of image
+    if (ro1) {
+        numCols = ro1->image->numCols;
+        numRows = ro1->image->numRows;
+    } else if (ro2) {
+        numCols = ro2->image->numCols;
+        numRows = ro2->image->numRows;
+    } else {
+        psAbort("No input image provided.");
+    }
+
+    pmSubtractionMaskInvalid(ro1, maskVal);
+    pmSubtractionMaskInvalid(ro2, maskVal);
+
+    // General background subtraction, since this is done in pmSubtractionMatch
+    {
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for background
+        if (ro1) {
+            psStatsInit(bg);
+            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
+                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
+                psFree(bg);
+                psFree(rng);
+                return false;
+            }
+            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(PM_ERR_DATA, false, "Unable to measure background statistics.");
+                psFree(bg);
+                psFree(rng);
+                return false;
+            }
+            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
+        }
+        psFree(bg);
+        psFree(rng);
+    }
+
+    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) {
@@ -283 +437 @@
     psMetadata *outHeader = psMetadataAlloc(); // Output header values
 
+    subtractionMatchAlloc(conv1, conv2, ro1, ro2, subMask, maskBad, mode, numCols, numRows);
+
     psTrace("psModules.imcombine", 2, "Convolving...\n");
     for (int i = 0; i < kernels->n; i++) {
@@ -288 +444 @@
         psRegion *region = regions->data[i]; // Region of interest
 
-        if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region,
-                                   ro1->image->numCols, ro1->image->numRows)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
+        if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region, numCols, numRows)) {
+            psError(psErrorCodeLast(), false, "Unable to generate QA data");
             psFree(outAnalysis);
             psFree(outHeader);
@@ -300 +455 @@
 
         if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
-                                   sysError, region, kernel, true, useFFT)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
+                                   kernelError, covarFrac, region, kernel, true, useFFT)) {
+            psError(psErrorCodeLast(), false, "Unable to convolve image.");
             psFree(outAnalysis);
             psFree(outHeader);
@@ -330 +485 @@
                         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,
+                        int iter, float rej, float normFrac, float sysError, float skyError,
+                        float kernelError, float covarFrac, 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)) {
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, normFrac, sysError, skyError, kernelError,
+                               covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
         return false;
     }
+
+    // We need both inputs
+    PM_ASSERT_READOUT_NON_NULL(ro1, false);
+    PM_ASSERT_READOUT_NON_NULL(ro2, false);
 
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);
@@ -392 +552 @@
     // 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
@@ -405 +565 @@
     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.");
@@ -416 +580 @@
     // 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.");
+                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
                 psFree(bg);
                 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.");
+                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
                 psFree(bg);
                 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);
     }
+
+    subtractionMatchAlloc(conv1, conv2, ro1, ro2, subMask, maskBad, subMode, numCols, numRows);
 
     // Iterate over iso-kernel regions
@@ -452 +625 @@
             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, subMode);
+                                                        footprint, stampSpacing, normFrac,
+                                                        sysError, skyError, subMode);
             } else if (sources) {
                 stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
-                                                           footprint, stampSpacing, 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, size, footprint, subMode)) {
+            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
+                                      stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
                 goto MATCH_ERROR;
             }
 
+
+            // generate the window function from the set of stamps
+            if (!pmSubtractionStampsGetWindow(stamps, size)) {
+                psError(psErrorCodeLast(), false, "Unable to get stamp window.");
+                goto MATCH_ERROR;
+            }
+
+            // 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, bounds, subMode);
+                if (!kernels) {
+                    psErrorClear();
+                    psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
+                }
+            }
+            if (kernels == NULL) {
+                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
+                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
+                                                       inner, binning, ringsOrder, penalty, bounds, subMode);
+                // pmSubtractionVisualShowKernels(kernels);
+            }
+
+            memCheck("kernels");
+
             if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
+#if 0
                 // Get backgrounds
                 psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background
                 psVector *buffer = NULL;// Buffer for stats
                 if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskVal, rng)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 1.");
+                    psError(PM_ERR_DATA, false, "Unable to measure background of image 1.");
                     psFree(bgStats);
                     psFree(buffer);
@@ -488 +690 @@
                 float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1
                 if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskVal, rng)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 2.");
+                    psError(PM_ERR_DATA, false, "Unable to measure background of image 2.");
                     psFree(bgStats);
                     psFree(buffer);
@@ -498 +700 @@
 
                 pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use
+#endif
+                pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej);
                 switch (newMode) {
                   case PM_SUBTRACTION_MODE_1:
@@ -506 +710 @@
                     break;
                   default:
-                    psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order.");
+                    psError(psErrorCodeLast(), false, "Unable to determine subtraction order.");
                     goto MATCH_ERROR;
                 }
                 subMode = newMode;
             }
-
-            // 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);
-                if (!kernels) {
-                    psErrorClear();
-                    psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
-                }
-            }
-            if (kernels == NULL) {
-                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
-                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
-                                                       inner, binning, ringsOrder, penalty, subMode);
-            }
-
-            memCheck("kernels");
 
             int numRejected = -1;       // Number of rejected stamps in each iteration
@@ -534 +721 @@
 
                 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
-                                          stampThresh1, stampThresh2, stampSpacing,
-                                          size, footprint, subMode)) {
-                    goto MATCH_ERROR;
-                }
-
-                psTrace("psModules.imcombine", 3, "Calculating equation...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-
+                                          stampThresh1, stampThresh2, stampSpacing, normFrac,
+                                          sysError, skyError, size, footprint, subMode)) {
+                    goto MATCH_ERROR;
+                }
+
+                // generate the window function from the set of stamps
+                if (!pmSubtractionStampsGetWindow(stamps, size)) {
+                    psError(psErrorCodeLast(), false, "Unable to get stamps window.");
+                    goto MATCH_ERROR;
+                }
+
+                // XXX step 1: calculate normalization
+                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+
+                psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+                memCheck("  solve equation");
+
+# if (SEPARATE)
+                // set USED -> CALCULATE
+                pmSubtractionStampsResetStatus (stamps);
+
+                // XXX step 2: calculate kernel parameters
+                psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
                 memCheck("  calculate equation");
 
-                psTrace("psModules.imcombine", 3, "Solving equation...\n");
-
-                if (!pmSubtractionSolveEquation(kernels, stamps)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-
+                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
                 memCheck("  solve equation");
-
+# endif
                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+                    psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
                     goto MATCH_ERROR;
                 }
@@ -565 +774 @@
 
                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
-                numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint);
+                numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
                 if (numRejected < 0) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
+                    psError(psErrorCodeLast(), false, "Unable to reject stamps.");
                     psFree(deviations);
                     goto MATCH_ERROR;
@@ -576 +785 @@
             }
 
+            // if we hit the max number of iterations and we have rejected stamps, re-solve
             if (numRejected > 0) {
-                psTrace("psModules.imcombine", 3, "Solving equation...\n");
-                if (!pmSubtractionSolveEquation(kernels, stamps)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
+                // XXX step 1: calculate normalization
+                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+
+                // solve normalization
+                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+
+# if (SEPARATE)
+                // set USED -> CALCULATE
+                pmSubtractionStampsResetStatus (stamps);
+
+                // XXX step 2: calculate kernel parameters
+                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+
+                // solve kernel parameters
+                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+                memCheck("  solve equation");
+# endif
                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
-                    goto MATCH_ERROR;
-                }
-                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint);
+                    psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
+                    goto MATCH_ERROR;
+                }
+                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
                 psFree(deviations);
             }
@@ -596 +834 @@
 
             if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
+                psError(psErrorCodeLast(), false, "Unable to generate QA data");
                 goto MATCH_ERROR;
             }
@@ -604 +842 @@
             psTrace("psModules.imcombine", 2, "Convolving...\n");
             if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
-                                       sysError, region, kernels, true, useFFT)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
+                                       kernelError, covarFrac, region, kernels, true, useFFT)) {
+                psError(psErrorCodeLast(), false, "Unable to convolve image.");
                 goto MATCH_ERROR;
             }
@@ -624 +862 @@
     variance = NULL;
 
-    if (!pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
+    if (conv1 && !pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) {
+        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
+        goto MATCH_ERROR;
+    }
+    if (conv2 && !pmSubtractionBorder(conv2->image, conv2->variance, conv2->mask, size, maskBad)) {
+        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
         goto MATCH_ERROR;
     }
@@ -824 +1066 @@
         } else {
             if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to measure PSF width for stamp %d", i);
+                psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i);
                 psFree(models);
                 psFree(modelSums);
@@ -835 +1077 @@
 
     if (!psThreadPoolWait(true)) {
-        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        psError(psErrorCodeLast(), false, "Error waiting for threads.");
         psFree(models);
         psFree(modelSums);
@@ -848 +1090 @@
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     if (!psVectorStats(stats, ratios, NULL, mask, 0xff)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio.");
+        psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio.");
         psFree(mask);
         psFree(ratios);
@@ -869 +1111 @@
     return mode;
 }
+
+
+// Test a subtraction mode by performing a single iteration
+static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode
+                                pmSubtractionKernels *kernels, // Kernel description
+                                const char *description, // Description for trace
+                                psImage *subMask,  // Subtraction mask
+                                float rej               // Rejection threshold
+                                )
+{
+    assert(stamps);
+    assert(kernels);
+
+    psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
+    if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
+        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description);
+    if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
+        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+
+# if (SEPARATE)
+    // set USED -> CALCULATE
+    pmSubtractionStampsResetStatus (stamps);
+
+    psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description);
+    if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
+        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description);
+    if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+# endif
+
+    psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description);
+    psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
+    if (!deviations) {
+        psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description);
+    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
+    if (numRejected < 0) {
+        psError(psErrorCodeLast(), false, "Unable to reject stamps.");
+        psFree(deviations);
+        return false;
+    }
+    psFree(deviations);
+
+    if (numRejected > 0) {
+        // Allow re-fit with reduced stamps set
+        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
+        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
+            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+            return false;
+        }
+
+        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
+        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
+            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+            return false;
+        }
+        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
+
+# if (SEPARATE)
+        // set USED -> CALCULATE
+        pmSubtractionStampsResetStatus (stamps);
+
+        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
+        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
+            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+            return false;
+        }
+
+        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
+        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+            return false;
+        }
+        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
+# endif
+
+        psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
+        if (!deviations) {
+            psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
+            return false;
+        }
+        psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description);
+        long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
+        if (numRejected < 0) {
+            psError(psErrorCodeLast(), false, "Unable to reject stamps.");
+            psFree(deviations);
+            return false;
+        }
+        psFree(deviations);
+    }
+
+    return true;
+}
+
+
+pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels,
+                                        const psImage *subMask, float rej)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(*kernels, PM_SUBTRACTION_MODE_ERR);
+
+    // Copies of the inputs
+    pmSubtractionStampList *stamps1 = pmSubtractionStampListCopy(*stamps);
+    pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(*kernels);
+    psImage *subMask1 = psImageCopy(NULL, subMask, subMask->type.type);
+    kernels1->mode = PM_SUBTRACTION_MODE_1;
+
+    if (!subtractionModeTest(stamps1, kernels1, "convolve 1", subMask1, rej)) {
+        psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1");
+        psFree(stamps1);
+        psFree(kernels1);
+        psFree(subMask1);
+        return PM_SUBTRACTION_MODE_ERR;
+    }
+    psFree(subMask1);
+
+    // Copies of the inputs
+    pmSubtractionStampList *stamps2 = pmSubtractionStampListCopy(*stamps);
+    pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(*kernels);
+    psImage *subMask2 = psImageCopy(NULL, subMask, subMask->type.type);
+    kernels2->mode = PM_SUBTRACTION_MODE_2;
+
+    if (!subtractionModeTest(stamps2, kernels2, "convolve 2", subMask2, rej)) {
+        psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2");
+        psFree(stamps2);
+        psFree(kernels2);
+        psFree(subMask2);
+        psFree(stamps1);
+        psFree(kernels1);
+        return PM_SUBTRACTION_MODE_ERR;
+    }
+    psFree(subMask2);
+
+
+    pmSubtractionStampList *bestStamps = NULL; // Best choice for stamps
+    pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels
+    psLogMsg("psModules.imcombine", PS_LOG_INFO,
+             "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n",
+             kernels1->mean, kernels1->rms, kernels1->numStamps,
+             kernels2->mean, kernels2->rms, kernels2->numStamps);
+
+    if (kernels1->mean < kernels2->mean) {
+        bestStamps = stamps1;
+        bestKernels = kernels1;
+    } else {
+        bestStamps = stamps2;
+        bestKernels = kernels2;
+    }
+
+    psFree(*stamps);
+    psFree(*kernels);
+    *stamps = psMemIncrRefCounter(bestStamps);
+    *kernels = psMemIncrRefCounter(bestKernels);
+
+    psFree(stamps1);
+    psFree(stamps2);
+    psFree(kernels1);
+    psFree(kernels2);
+
+    return bestKernels->mode;
+}
+
+
+bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
+                              float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax)
+{
+    PS_ASSERT_PTR_NON_NULL(kernelSize, false);
+    PS_ASSERT_PTR_NON_NULL(stampSize, false);
+    PS_ASSERT_VECTOR_NON_NULL(widths, false);
+    PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
+
+//    float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
+    float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
+
+    if (isfinite(scaleMin) && scale < scaleMin) {
+        scale = scaleMin;
+    }
+    if (isfinite(scaleMax) && scale > scaleMax) {
+        scale = scaleMax;
+    }
+
+    for (int i = 0; i < widths->n; i++) {
+        widths->data.F32[i] *= scale;
+    }
+    *kernelSize = *kernelSize * scale + 0.5;
+    *stampSize = *stampSize * scale + 0.5;
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO,
+             "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
+
+    return true;
+}