Changeset 30622

Timestamp:

Feb 13, 2011, 12:19:53 PM (15 years ago)

Author:

eugene

Message:

some code reorganzation (move all types into pmSubtractionTypes.h); add stage to check on different mode and order options, choosing best based on chisq of subtraction; careful with the window and the kernel sizes: if the measured kron radius is too large, allow the window to grow as needed; use the 1st radial moment measured on the stamps to set the kernel scaling; remove some cruft from the code; calling program now passes in kernel scaling options to be used when 1st radial moment is measured; new function to update the kernel description string used for kernel I/O; update the kernel description after the spatial order and direction is chosen; need to carry the kernel fwhms and spatial order to allow for update; calculate the psf-match solution errors; psf solution now uses the weights to get sensible chisq values, window is used to calculate the moments; penalty scale is now adjusted to be consistent with sensible reduced chisq; improved visualizations; use range-limiter in SVD of 1e-10; calculate chisq and moments for the solution and use in the evaluation; split out the stamp selection / convolution steps; reject sources after fitting a model to the chisq assuming non-zero systematic error

Location:

trunk/psModules/src/imcombine

Files:

• Makefile.am (modified) (1 diff)
• pmReadoutCombine.c (modified) (1 diff)
• pmStackReject.c (modified) (1 diff)
• pmSubtraction.c (modified) (14 diffs)
• pmSubtraction.h (modified) (4 diffs)
• pmSubtractionAnalysis.c (modified) (1 diff)
• pmSubtractionDeconvolve.c (modified) (1 diff)
• pmSubtractionEquation.c (modified) (42 diffs)
• pmSubtractionEquation.h (modified) (3 diffs)
• pmSubtractionEquation.v0.c (modified) (1 diff)
• pmSubtractionHermitian.c (modified) (1 diff)
• pmSubtractionIO.c (modified) (1 diff)
• pmSubtractionKernels.c (modified) (18 diffs)
• pmSubtractionKernels.h (modified) (3 diffs)
• pmSubtractionMask.c (modified) (1 diff)
• pmSubtractionMatch.c (modified) (28 diffs)
• pmSubtractionMatch.h (modified) (3 diffs)
• pmSubtractionParams.c (modified) (2 diffs)
• pmSubtractionStamps.c (modified) (12 diffs)
• pmSubtractionStamps.h (modified) (3 diffs)
• pmSubtractionThreads.c (modified) (2 diffs)
• pmSubtractionTypes.h (copied) (copied from branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionTypes.h ) (2 diffs)
• pmSubtractionVisual.c (modified) (13 diffs)
• pmSubtractionVisual.h (modified) (1 diff)

trunk/psModules/src/imcombine/Makefile.am

@@ -30 +30 @@
         pmStackReject.h         \
         pmSubtraction.h         \
+        pmSubtractionTypes.h            \
         pmSubtractionAnalysis.h \
         pmSubtractionEquation.h \

trunk/psModules/src/imcombine/pmReadoutCombine.c

@@ -33 +33 @@
     params->blank = 0;
     params->nKeep = 0;
-    params->fracHigh = 0.0;
+    params->fracLow = 0.0;
     params->fracHigh = 0.0;
     params->iter = 1;

trunk/psModules/src/imcombine/pmStackReject.c

@@ -7 +7 @@
 #include <pslib.h>
 
+#include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionThreads.h"

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -20 +20 @@
 #include "pmHDU.h"                      // Required for pmFPA.h
 #include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionEquation.h"
 #include "pmSubtractionVisual.h"
 #include "pmSubtractionThreads.h"
+
+bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header);
 
 #include "pmSubtraction.h"
@@ -773 +776 @@
 
     if (convolutions) {
-        // Already done
         return convolutions;
     }
@@ -787 +789 @@
 }
 
+
+bool pmSubtractionConvolveStampThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    pmSubtractionStamp *stamp = job->args->data[0]; // List of stamps
+    pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
+    int footprint = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
+
+    return pmSubtractionConvolveStamp(stamp, kernels, footprint);
+}
 
 bool pmSubtractionConvolveStamp (pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int footprint)
@@ -818 +831 @@
     }
 
+#ifdef TESTING
+    for (int j = 0; j < kernels->num; j++) {
+        if (stamp->convolutions1) {
+            psString convName = NULL;
+            psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
+            psFits *fits = psFitsOpen(convName, "w");
+            psFree(convName);
+            psKernel *conv = stamp->convolutions1->data[j];
+            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
+            psFitsClose(fits);
+        }
+
+        if (stamp->convolutions2) {
+            psString convName = NULL;
+            psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
+            psFits *fits = psFitsOpen(convName, "w");
+            psFree(convName);
+            psKernel *conv = stamp->convolutions2->data[j];
+            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
+            psFitsClose(fits);
+        }
+    }
+#endif
+
     return true;
 }
 
+bool pmSubtractionConvolveStamps(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) 
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
+
+    psTimerStart("pmSubtractionConvolveStamps");
+
+    int footprint = stamps->footprint;  // Half-size of stamps
+
+    // We iterate over each stamp and generate the convolution if needed.  We do NOT need the
+    // convolution if (a) it has already been calculated or (b) the stamp is not available for
+    // use (available = USED or CALCULATE)
+    
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+
+        bool keep = false;
+        keep |= (stamp->status == PM_SUBTRACTION_STAMP_USED);
+        keep |= (stamp->status == PM_SUBTRACTION_STAMP_CALCULATE);
+        if (!keep) continue;
+
+        bool haveConvolutions = false;
+        if (kernels->mode == PM_SUBTRACTION_MODE_1) {
+            haveConvolutions = (stamp->convolutions1 != NULL);
+        }
+        if (kernels->mode == PM_SUBTRACTION_MODE_2) {
+            haveConvolutions = (stamp->convolutions2 != NULL);
+        }
+        if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+            haveConvolutions = (stamp->convolutions1 != NULL) && (stamp->convolutions2 != NULL);
+        }
+        if (haveConvolutions) {
+            continue;
+        }
+
+        if (pmSubtractionThreaded()) {
+            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CONVOLVE_STAMP");
+            psArrayAdd(job->args, 1, stamp);
+            psArrayAdd(job->args, 1, kernels);
+            PS_ARRAY_ADD_SCALAR(job->args, footprint, PS_TYPE_S32);
+            if (!psThreadJobAddPending(job)) {
+                return false;
+            }
+        } else {
+            pmSubtractionConvolveStamp(stamp, kernels, footprint);
+        }
+    }
+    if (!psThreadPoolWait(true)) {
+        psError(psErrorCodeLast(), false, "Error waiting for threads.");
+        return false;
+    }
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolve stamps: %f sec", psTimerClear("pmSubtractionConvolveStamps"));
+    return true;
+}
 
 int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps,
-                              const psVector *deviations, psImage *subMask, float sigmaRej)
+                              pmSubtractionQuality *match, psImage *subMask, float sigmaRej)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1);
-    PS_ASSERT_VECTOR_NON_NULL(deviations, -1);
-    PS_ASSERT_VECTOR_TYPE(deviations, PS_TYPE_F32, -1);
     PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1);
     PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_IMAGE_MASK, -1);
 
-    // I used to measure the rms deviation about zero, and use that as the sigma against which to clip, but
-    // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like
-    // with large N.  Therefore, let's just treat this as a Gaussian distribution.
+    // Comment from PAP (r18287): I used to measure the rms deviation about zero, and use that as the
+    // sigma against which to clip, but the distribution is actually something like a chi^2 or
+    // Student's t, both of which become Gaussian-like with large N.  Therefore, let's just
+    // treat this as a Gaussian distribution.
+
+    // Comment from EAM (r29777): The residual distribution is only chisq-like if the model is
+    // a good fit to the data.  In the (likely) case that there is a systematic difference
+    // between the model and the data, the squared-residual distribution grows quadratically
+    // with increasing flux: the systematic residual flux is a constant factor times the source
+    // flux; the squared-residual is then of the form (k0 + k1*flux)^2, where k0 comes from the
+    // Gaussian distributed residual and k1*flux is the systematic residual error.
+
+    // By rejecting sources with the largest squared-residuals, the rejection biases against
+    // the brighter sources; in severe cases, this pushes the measurement to the weakest
+    // sources with the most noise.  To account for this, let's fit a 2nd order polynomial to
+    // the distribution of flux vs squared-residual, subtract that fit, and reject sources
+    // which are significantly deviant from that distribution.
 
     kernels->mean = NAN;
@@ -840 +943 @@
     kernels->numStamps = -1;
 
-    int numStamps = 0;                  // Number of used stamps
-    psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_VECTOR_MASK); // Mask, for statistics
-    psVectorInit(mask, 0);
-    for (int i = 0; i < stamps->num; i++) {
-        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
-            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
-            continue;
-        }
-        numStamps++;
-    }
-    psTrace("psModules.imcombine", 1, "Number of good stamps: %d\n", numStamps);
-
-    if (numStamps == 0) {
-        psError(PM_ERR_STAMPS, true, "No good stamps found.");
-        psFree(mask);
-        return -1;
-    }
-
-    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV |
-                                  PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics for deviatns
-    if (!psVectorStats(stats, deviations, NULL, mask, 0xff)) {
-        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
+    psTrace("psModules.imcombine", 1, "Number of good stamps: %d\n", match->nGood);
+
+    // the chisq & flux vectors are calculated by pmSubtractionCalculateChisqAndMoments
+
+    // use 3hi/3lo sigma clipping on the chisq fit
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    stats->clipSigma = 5.0;
+    stats->clipIter = 2;
+    psPolynomial1D *model = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
+
+    bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, match->chisq, NULL, match->fluxes);
+    if (!result) {
+        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
+        psFree(model);
         psFree(stats);
-        psFree(mask);
-        return -1;
-    }
-    psFree(mask);
-
-    // XXX raise an error?
+        return -1;
+    }
     if (isnan(stats->sampleMean)) {
+        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
+        psFree(model);
         psFree(stats);
         return -1;
     }
 
-    double mean, rms;                 // Mean and RMS of deviations
-    if (numStamps < MIN_SAMPLE_STATS) {
-        mean = stats->sampleMean;
-        rms = stats->sampleStdev;
-    } else {
-        mean = stats->sampleMedian;
-        rms = 0.74 * (stats->sampleUQ - stats->sampleLQ);
-    }
-    psFree(stats);
-
-    psTrace("psModules.imcombine", 1, "Mean: %f\n", mean);
-    psTrace("psModules.imcombine", 1, "RMS deviation: %f\n", rms);
-
-    kernels->mean = mean;
-    kernels->rms = rms;
-    kernels->numStamps = numStamps;
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf",
-             numStamps, mean, rms);
-
-    if (!isfinite(sigmaRej) || sigmaRej <= 0.0) {
-        // User just wanted to calculate and record the deviation for posterity
-        return 0;
-    }
-
-    float limit = sigmaRej * rms; // Limit on maximum deviation
-    psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
-
+    kernels->mean = stats->sampleMean;
+    kernels->rms = stats->sampleStdev;
+    kernels->numStamps = stats->clippedNvalues;
+
+    psLogMsg ("pmPSFtry", 4, "chisq vs flux model: %e + %e flux + %e flux^2\n", model->coeff[0], model->coeff[1], model->coeff[2]);
+    psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", stats->sampleMean, stats->sampleStdev);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf",  kernels->numStamps, kernels->mean, kernels->rms);
 
     psString ds9name = NULL;            // Filename for ds9 region file
@@ -914 +985 @@
     int numRejected = 0;                // Number of stamps rejected
     int numGood = 0;                    // Number of good stamps
-    double newMean = 0.0;               // New mean
     psString log = NULL;                // Log message
-    psStringAppend(&log, "Rejecting stamps, mean = %f, threshold = %f\n", mean, limit);
+
+    // save DS9 region files for the stamps and mark for rejection and replacement
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
+        if (stamp->status  != PM_SUBTRACTION_STAMP_USED) { continue; }
+        if (match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
             // Should we reject stars with low deviation?  Well, if this is really a Gaussian-like
             // distribution and they're low, then we have the right to ask why.  Isn't it suspicious that
@@ -926 +998 @@
             // subtract well, in which case very few (if any) stars will be legitimately rejected for being
             // low.
-            if (fabsf(deviations->data.F32[i] - mean) > limit) {
-                // Mask out the stamp in the image so you it's not found again
-                psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
-                        (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
-                psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i,
-                               (int)(stamp->x - 0.5), (int)(stamp->y - 0.5),
-                               fabsf(deviations->data.F32[i] - mean));
-                numRejected++;
-                for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
-                    for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
-                        subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
-                    }
-                }
-                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
-
-                // Set stamp for replacement
-                stamp->x = 0;
-                stamp->y = 0;
-                stamp->xNorm = NAN;
-                stamp->yNorm = NAN;
-                stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-                // Recalculate convolutions
-                psFree(stamp->convolutions1);
-                psFree(stamp->convolutions2);
-                stamp->convolutions1 = stamp->convolutions2 = NULL;
-                psFree(stamp->image1);
-                psFree(stamp->image2);
-                psFree(stamp->weight);
-                stamp->image1 = stamp->image2 = stamp->weight = NULL;
-                psFree(stamp->matrix);
-                stamp->matrix = NULL;
-                psFree(stamp->vector);
-                stamp->vector = NULL;
-            } else {
-                numGood++;
-                newMean += deviations->data.F32[i];
-                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
-            }
-        }
-    }
-    newMean /= numGood;
+            psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
+                    (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
+            psStringAppend(&log, "Stamp %d (%d,%d): %f : %f : %f\n", 
+                           i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5),
+                           match->chisq->data.F32[i], match->fluxes->data.F32[i], match->chisq->data.F32[i] - psPolynomial1DEval(model, match->fluxes->data.F32[i])); 
+            numRejected++;
+            for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
+                for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
+                    subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
+                }
+            }
+            pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
+
+            // Set stamp for replacement
+            stamp->x = 0;
+            stamp->y = 0;
+            stamp->xNorm = NAN;
+            stamp->yNorm = NAN;
+            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+            // Recalculate convolutions
+            psFree(stamp->convolutions1);
+            psFree(stamp->convolutions2);
+            stamp->convolutions1 = stamp->convolutions2 = NULL;
+            psFree(stamp->image1);
+            psFree(stamp->image2);
+            psFree(stamp->weight);
+            stamp->image1 = stamp->image2 = stamp->weight = NULL;
+            psFree(stamp->matrix);
+            stamp->matrix = NULL;
+            psFree(stamp->vector);
+            stamp->vector = NULL;
+        } else {
+            numGood++;
+            pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
+        }
+    }
 
     if (numRejected == 0) {
@@ -978 +1045 @@
     }
 
+    psFree(model);
+    psFree(stats);
+
     if (numRejected > 0) {
-        psLogMsg("psModules.imcombine", PS_LOG_INFO,
-                 "%d good stamps; %d rejected.\nMean deviation: %lf --> %lf\n",
-                 numGood, numRejected, mean, newMean);
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; %d rejected.\n", numGood, numRejected);
     } else {
-        psLogMsg("psModules.imcombine", PS_LOG_INFO,
-                 "%d good stamps; 0 rejected.\nMean deviation: %lf\n",
-                 numGood, mean);
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; 0 rejected.\n", numGood);
     }
 
@@ -1373 +1439 @@
     psFree(kernelErr2);
 
+    static int nOut1 = 0;
+    static int nOut2 = 0;
+
     // Calculate covariances
     // This can be fairly involved, so we only do it for a small number of instances
@@ -1386 +1455 @@
                 psKernelTruncate(kernel, covarFrac);
                 covars->data[i] = psImageCovarianceCalculate(kernel, ro1->covariance);
+                if (0) {
+                    char name[128];
+                    snprintf (name, 128, "covar.sample1.%03d.fits", nOut1);
+                    psKernel *cov = covars->data[i];
+                    psFitsWriteImageSimple (name, cov->image, NULL);
+
+                    snprintf (name, 128, "incovar.sample1.%03d.fits", nOut1);
+                    psFitsWriteImageSimple (name, ro1->covariance->image, NULL);
+
+                    snprintf (name, 128, "conv.sample1.%03d.fits", nOut1);
+                    psFitsWriteImageSimple (name, kernel->image, NULL);
+
+                    fprintf (stderr, "incov: %d,%d; kern: %d,%d, outcov: %d,%d\n", 
+                             ro1->covariance->image->numCols, ro1->covariance->image->numRows, 
+                             kernel->image->numCols, kernel->image->numRows,
+                             cov->image->numCols, cov->image->numRows);
+
+                    nOut1 ++;
+                }
                 psFree(kernel);
             }
@@ -1407 +1495 @@
                 psKernelTruncate(kernel, covarFrac);
                 covars->data[i] = psImageCovarianceCalculate(kernel, ro2->covariance);
+                if (0) {
+                    char name[128];
+                    snprintf (name, 128, "covar.sample2.%03d.fits", nOut2);
+                    psKernel *cov = covars->data[i];
+                    psFitsWriteImageSimple (name, cov->image, NULL);
+
+                    snprintf (name, 128, "incovar.sample2.%03d.fits", nOut2);
+                    psFitsWriteImageSimple (name, ro2->covariance->image, NULL);
+
+                    snprintf (name, 128, "conv.sample2.%03d.fits", nOut2);
+                    psFitsWriteImageSimple (name, kernel->image, NULL);
+
+                    fprintf (stderr, "incov: %d,%d; kern: %d,%d, outcov: %d,%d\n", 
+                             ro2->covariance->image->numCols, ro2->covariance->image->numRows, 
+                             kernel->image->numCols, kernel->image->numRows,
+                             cov->image->numCols, cov->image->numRows);
+
+                    nOut2 ++;
+                }
                 psFree(kernel);
             }
@@ -1422 +1529 @@
     psImageCovarianceSetThreads(oldThreads);
 
-    // Copy anything that wasn't convolved
+    // Copy anything that wasn't convolved (they may have been allocated though, so free them)
     switch (kernels->mode) {
       case PM_SUBTRACTION_MODE_1:
         if (out2) {
+            psFree(out2->image);
+            psFree(out2->variance);
+            psFree(out2->mask);
+            psFree(out2->covariance);
             out2->image = psMemIncrRefCounter(ro2->image);
             out2->variance = psMemIncrRefCounter(ro2->variance);
@@ -1434 +1545 @@
       case PM_SUBTRACTION_MODE_2:
         if (out1) {
+            psFree(out1->image);
+            psFree(out1->variance);
+            psFree(out1->mask);
+            psFree(out1->covariance);
             out1->image = psMemIncrRefCounter(ro1->image);
             out1->variance = psMemIncrRefCounter(ro1->variance);
@@ -1479 +1594 @@
   return true;
 }
+
+static void pmSubtractionQualityFree(pmSubtractionQuality *quality) {
+
+    psFree (quality->fluxes);
+    psFree (quality->chisq);
+    psFree (quality->moments);
+    psFree (quality->stampMask);
+}    
+
+pmSubtractionQuality *pmSubtractionQualityAlloc() {
+
+    pmSubtractionQuality *quality = psAlloc(sizeof(pmSubtractionQuality)); // Stamp list to return
+    psMemSetDeallocator(quality, (psFreeFunc)pmSubtractionQualityFree);
+
+    quality->fluxes = NULL;
+    quality->chisq = NULL;
+    quality->moments = NULL;
+    quality->stampMask = NULL;
+
+    quality->score = NAN;
+    quality->mode = PM_SUBTRACTION_MODE_ERR;
+    quality->spatialOrder = -1;
+    quality->nGood = 0;
+    
+    return quality;
+}

trunk/psModules/src/imcombine/pmSubtraction.h

@@ -14 +14 @@
 #define PM_SUBTRACTION_H
 
-#include <pslib.h>
-
-#include <pmHDU.h>
-#include <pmFPA.h>
-#include <pmSubtractionKernels.h>
-#include <pmSubtractionStamps.h>
+// #include <pslib.h>
+// #include <pmHDU.h>
+// #include <pmFPA.h>
+// #include <pmSubtractionKernels.h>
+// #include <pmSubtractionStamps.h>
 
 // if we use the original ppSub implementation, we subtract a central delta-function for all
@@ -30 +29 @@
 /// @addtogroup imcombine Image Combinations
 /// @{
-
-/// Mask values for the subtraction mask
-typedef enum {
-    PM_SUBTRACTION_MASK_CLEAR          = 0x00, // No masking
-    PM_SUBTRACTION_MASK_BAD_1          = 0x01, // Image 1 is bad
-    PM_SUBTRACTION_MASK_BAD_2          = 0x02, // Image 2 is bad
-    PM_SUBTRACTION_MASK_CONVOLVE_1     = 0x04, // If image 1 is convolved, would be poor or bad
-    PM_SUBTRACTION_MASK_CONVOLVE_2     = 0x08, // If image 2 is convolved, would be poor or bad
-    PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 = 0x10, // If image 1 is convolved, would be bad
-    PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 = 0x20, // If image 2 is convolved, would be bad
-    PM_SUBTRACTION_MASK_BORDER         = 0x40, // Image border
-    PM_SUBTRACTION_MASK_REJ            = 0x80, // Previously tried as a stamp, and rejected
-} pmSubtractionMasks;
-
 
 /// Number of terms in a polynomial
@@ -70 +55 @@
     );
 
+bool pmSubtractionConvolveStamps(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels);
+
+bool pmSubtractionConvolveStampThread(psThreadJob *job);
+
 /// Reject stamps
 int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, ///< Kernel parameters to update
                               pmSubtractionStampList *stamps, ///< Stamps
-                              const psVector *deviations, ///< Deviations for each stamp
+                              pmSubtractionQuality *match, ///< data on the subtraction quality
                               psImage *subMask, ///< Subtraction mask
                               float sigmaRej ///< Number of RMS deviations above zero at which to reject
@@ -167 +156 @@
 bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2);
 
+pmSubtractionQuality *pmSubtractionQualityAlloc();
+
 /// @}
 #endif

trunk/psModules/src/imcombine/pmSubtractionAnalysis.c

@@ -7 +7 @@
 
 #include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"

trunk/psModules/src/imcombine/pmSubtractionDeconvolve.c

@@ -10 +10 @@
 
 #include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtractionKernels.h"
 #include "pmSubtractionDeconvolve.h"

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -8 +8 @@
 
 #include "pmErrorCodes.h"
+#include "pmVisual.h"
+#include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"
@@ -16 +19 @@
 #include "pmSubtractionVisual.h"
 
-//#define TESTING                         // TESTING output for debugging; may not work with threads!
-
-//#define USE_WEIGHT                      // Include weight (1/variance) in equation?
-
-// XXX TEST:
-//# define USE_WINDOW                      // window to avoid neighbor contamination
+//# define TESTING                         // TESTING output for debugging; may not work with threads!
+# define USE_WEIGHT                      // Include weight (1/variance) in equation?
+# define USE_WINDOW                      // window to avoid neighbor contamination
+
+/* I believe we want to apply the WEIGHT to the chisq portions of the calculation (but not the WINDOW),
+ * and the WINDOW to the moments portiosn of the calculations (but not the WEIGHT)
+ *
+ */
 
 # define PENALTY false
 # define MOMENTS (!PENALTY)
-# define MOMENTS_PENALTY_SCALE 2e-5 // XXX this value is not completely arbitrary, but I don't understand why it needs to be this value...
+# define MOMENTS_PENALTY_SCALE 20 // up-weight the moments somewhat
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -107 +112 @@
                         cc *= weight->kernel[y][x];
                     }
-                    if (window) {
+                    // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
+                    if (false && window) {
                         cc *= window->kernel[y][x];
                     }
@@ -138 +144 @@
                     rc *= wtVal;
                 }
-                if (window) {
+                // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
+                if (false && window) {
                     float winVal = window->kernel[y][x];
                     ic *= winVal;
@@ -173 +180 @@
 }
 
+# define RENORM_BY_FLUX 0
 
 // Calculate the least-squares matrix and vector for dual convolution
@@ -268 +276 @@
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
+
+                    // XXX NOTE: clipping low S/N pixels does not seem to work very well
+                    if (false && weight) {
+                        float i1 = image1->kernel[y][x];
+                        float i2 = image2->kernel[y][x];
+                        float sn = (i1 + i2) / sqrt (weight->kernel[y][x]);
+                        if (sn < 0.5) continue;
+                    }
+
                     double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] * PS_SQR(normValue);
                     double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue;
-                    if (weight) {
-                        float wtVal = weight->kernel[y][x];
-                        aa *= wtVal;
-                        bb *= wtVal;
-                        ab *= wtVal;
-                    }
-                    if (window) {
-                        float wtVal = window->kernel[y][x];
-                        aa *= wtVal;
-                        bb *= wtVal;
-                        ab *= wtVal;
-                    }
-                    sumAA += aa;
-                    sumBB += bb;
-                    sumAB += ab;
+
+                    float wtVal = (weight) ? weight->kernel[y][x] : 1.0;
+                    sumAA += wtVal*aa;
+                    sumBB += wtVal*bb;
+                    sumAB += wtVal*ab;
 
                     if (MOMENTS) {
-                        MxxAA += x*x*aa;
-                        MyyAA += y*y*aa;
-                        MxxBB += x*x*bb;
-                        MyyBB += y*y*bb;
+                        float winVal = (window) ? window->kernel[y][x] : 1.0;
+                        MxxAA += winVal*x*x*aa;
+                        MyyAA += winVal*y*y*aa;
+                        MxxBB += winVal*x*x*bb;
+                        MyyBB += winVal*y*y*bb;
                     }
                 }
             }
 
+            // XXX does normSquare1,2 mess up the relative scaling?
+            // XXX no: normSquare1,2 is the sum of the flux^2 for the source
             if (MOMENTS) {
                 MxxAA /= stamp->normSquare1 * PS_SQR(normValue);
@@ -305 +315 @@
             // XXX this makes the Chisq portion independent of the normalization and star flux
             // but may be mis-scaling between stars of different fluxes
+# if (RENORM_BY_FLUX)       
             sumAA /= PS_SQR(stamp->normI2);
             sumAB /= PS_SQR(stamp->normI2);
             sumBB /= PS_SQR(stamp->normI2);
+# endif
 
             // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB);
@@ -328 +340 @@
                         matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
                         matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
-
                         matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
                         matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
@@ -334 +345 @@
                         matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
                         matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
-
                         matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
                         matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
@@ -354 +364 @@
                         ab *= weight->kernel[y][x];
                     }
-                    if (window) {
+                    // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
+                    if (false && window) {
                         ab *= window->kernel[y][x];
                     }
@@ -363 +374 @@
             // XXX this makes the Chisq portion independent of the normalization and star flux
             // but may be mis-scaling between stars of different fluxes
+# if (RENORM_BY_FLUX)
             sumAB /= PS_SQR(stamp->normI2);
+# endif
 
             // Spatial variation of kernel coeffs
@@ -391 +404 @@
                 float i2 = image2->kernel[y][x];
 
+                // XXX NOTE: clipping low S/N pixels does not seem to work very well
+                if (false && weight) {
+                    float sn = (i1 + i2) / sqrt (weight->kernel[y][x]);
+                    if (sn < 0.5) continue;
+                }
+
                 double ai2 = a * i2 * normValue;
                 double bi2 = b * i2;
@@ -396 +415 @@
                 double bi1 = b * i1 * normValue;
 
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    ai2 *= wtVal;
-                    bi2 *= wtVal;
-                    ai1 *= wtVal;
-                    bi1 *= wtVal;
-                }
-                if (window) {
-                    float wtVal = window->kernel[y][x];
-                    ai2 *= wtVal;
-                    bi2 *= wtVal;
-                    ai1 *= wtVal;
-                    bi1 *= wtVal;
-                }
-                sumAI2 += ai2;
-                sumBI2 += bi2;
-                sumAI1 += ai1;
-                sumBI1 += bi1;
+                float wtVal = (weight) ? weight->kernel[y][x] : 1.0;
+                sumAI2 += wtVal*ai2;
+                sumBI2 += wtVal*bi2;
+                sumAI1 += wtVal*ai1;
+                sumBI1 += wtVal*bi1;
 
                 if (MOMENTS) {
-                    MxxAI1 += x*x*ai1;
-                    MyyAI1 += y*y*ai1;
-                    MxxBI2 += x*x*bi2;
-                    MyyBI2 += y*y*bi2;
+                    float winVal = (window) ? window->kernel[y][x] : 1.0;
+                    MxxAI1 += winVal*x*x*ai1;
+                    MyyAI1 += winVal*y*y*ai1;
+                    MxxBI2 += winVal*x*x*bi2;
+                    MyyBI2 += winVal*y*y*bi2;
                 }
             }
@@ -435 +442 @@
         // XXX this makes the Chisq portion independent of the normalization and star flux
         // but may be mis-scaling between stars of different fluxes
+# if (RENORM_BY_FLUX)
         sumAI1 /= PS_SQR(stamp->normI2);
         sumBI1 /= PS_SQR(stamp->normI2);
         sumAI2 /= PS_SQR(stamp->normI2);
         sumBI2 /= PS_SQR(stamp->normI2);
+# endif
 
         // Spatial variation
@@ -788 +797 @@
     stamp->normSquare2 = normSquare2;
 
-    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f  (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2);
+    // psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f  (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2);
 
     return true;
@@ -800 +809 @@
     pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
     int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
-    pmSubtractionEquationCalculationMode mode  = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model
-
-    return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode);
-}
-
-bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
-                                         int index, const pmSubtractionEquationCalculationMode mode)
+
+    return pmSubtractionCalculateEquationStamp(stamps, kernels, index);
+}
+
+bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, int index)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -833 +840 @@
     psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state.");
 
-    // Generate convolutions: these are generated once and saved
-    if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
-        psError(psErrorCodeLast(), false, "Unable to convolve stamp %d.", index);
-        return NULL;
-    }
-
-#ifdef TESTING
-    for (int j = 0; j < numKernels; j++) {
-        if (stamp->convolutions1) {
-            psString convName = NULL;
-            psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
-            psFits *fits = psFitsOpen(convName, "w");
-            psFree(convName);
-            psKernel *conv = stamp->convolutions1->data[j];
-            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
-            psFitsClose(fits);
-        }
-
-        if (stamp->convolutions2) {
-            psString convName = NULL;
-            psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
-            psFits *fits = psFitsOpen(convName, "w");
-            psFree(convName);
-            psKernel *conv = stamp->convolutions2->data[j];
-            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
-            psFitsClose(fits);
-        }
-    }
-#endif
-
-    // XXX visualize the set of convolved stamps
-
-    psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder,
-                                                    stamp->xNorm, stamp->yNorm); // Polynomial terms
-
-    bool new = stamp->vector ? false : true; // Is this a new run?
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms
+
+    // Is this a new run? Have we allocated the correct sized vector/matrix?
+    bool new = stamp->vector ? false : true;
+    if (!new && (stamp->vector->n != numParams)) {
+        psFree (stamp->vector);
+        psFree (stamp->matrix);
+        new = true;
+    }
+
     if (new) {
         stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
@@ -941 +922 @@
 }
 
-bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
-                                    const pmSubtractionEquationCalculationMode mode)
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -953 +933 @@
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+
         if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
             continue;
@@ -969 +950 @@
             psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array
             PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32);
             if (!psThreadJobAddPending(job)) {
                 return false;
             }
         } else {
-            pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode);
+            pmSubtractionCalculateEquationStamp(stamps, kernels, i);
         }
     }
@@ -983 +963 @@
     }
 
-    pmSubtractionVisualPlotLeastSquares(stamps);
     pmSubtractionVisualShowKernels((pmSubtractionKernels  *)kernels);
     pmSubtractionVisualShowBasis(stamps);
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec",
-             psTimerClear("pmSubtractionCalculateEquation"));
-
+    pmSubtractionVisualPlotLeastSquares(stamps);
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", psTimerClear("pmSubtractionCalculateEquation"));
 
     return true;
@@ -998 +976 @@
 bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header);
 
-psImage *p_pmSubSolve_wUt (psVector *w, psImage *U);
-psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt);
-
-bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask);
-
-bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B);
-bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB);
-bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB);
-
-bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x);
-bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y);
-bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
-
-psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w);
-
-double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
-
-bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
-                                const pmSubtractionStampList *stamps,
-                                const pmSubtractionEquationCalculationMode mode)
+bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+
+    psTimerStart("pmSubtractionSolveEquation");
 
     // Check inputs
@@ -1042 +1003 @@
         }
 
+        if (stamp->vector->n != numParams) {
+            fprintf (stderr, "mismatch length\n");
+        }
         PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false);
         PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false);
@@ -1080 +1044 @@
         }
 
+        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
+
 #if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
@@ -1087 +1053 @@
 #endif
 
+        // XXX TEST : print the matrix & vector
+        if (0) {
+            for (int iy = 0; iy < sumMatrix->numRows; iy++) {
+                for (int ix = 0; ix < sumMatrix->numCols; ix++) {
+                    fprintf (stderr, "%e  ", sumMatrix->data.F64[iy][ix]);
+                }
+                fprintf (stderr, " : %e\n", sumVector->data.F64[iy]);
+            }
+        }
+
+        psImage *invMatrix = NULL;
         psVector *solution = NULL;                       // Solution to equation!
         solution = psVectorAlloc(numParams, PS_TYPE_F64);
@@ -1094 +1071 @@
         // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
         // SINGLE solution
-        if (1) {
-            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
-        } else {
-            psVector *PERM = NULL;
-            psImage *LU = psMatrixLUDecomposition(NULL, &PERM, sumMatrix);
-            solution = psMatrixLUSolution(solution, LU, sumVector, PERM);
-            psFree (LU);
-            psFree (PERM);
-        }
+# if (1)
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10);
+        invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
+# endif
+# if (0)
+        psMatrixLUSolve(sumMatrixLU, sumVector);
+        solution = psMemIncrRefCounter(sumVector);
+        invMatrix = psMemIncrRefCounter(sumMatrix);
+# endif
+# if (0)
+        psMatrixGJSolve(sumMatrix, sumVector);
+        invMatrix = psMemIncrRefCounter(sumMatrix);
+        solution = psMemIncrRefCounter(sumVector);
+# endif
 
 # if (0)
         for (int i = 0; i < solution->n; i++) {
-            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]);
+            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(fabs(invMatrix->data.F64[i][i])));
         }
 # endif
 
-        if (!kernels->solution1) {
-            kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
-            psVectorInit(kernels->solution1, 0.0);
-        }
+        // ensure we have a solution vector of the right size
+        kernels->solution1    = psVectorRecycle(kernels->solution1,    sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+        kernels->solution1err = psVectorRecycle(kernels->solution1err, sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+        psVectorInit(kernels->solution1, 0.0);
+        psVectorInit(kernels->solution1err, 0.0);
 
         int numKernels = kernels->num;
         int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
         int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
+
         for (int i = 0; i < numKernels * numPoly; i++) {
             kernels->solution1->data.F64[i] = solution->data.F64[i];
+            kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]);
         }
 
@@ -1131 +1116 @@
         psFree(sumVector);
         psFree(sumMatrix);
+        psFree(invMatrix);
 
     } else {
@@ -1160 +1146 @@
         }
 
+        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
+
 #if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
@@ -1171 +1159 @@
         }
 
+        // XXX TEST : print the matrix & vector
+        if (0) {
+            for (int iy = 0; iy < sumMatrix->numRows; iy++) {
+                for (int ix = 0; ix < sumMatrix->numCols; ix++) {
+                    fprintf (stderr, "%e  ", sumMatrix->data.F64[iy][ix]);
+                }
+                fprintf (stderr, " : %e\n", sumVector->data.F64[iy]);
+            }
+        }
+
+        psImage *invMatrix = NULL;
         psVector *solution = NULL;                       // Solution to equation!
         solution = psVectorAlloc(numParams, PS_TYPE_F64);
@@ -1176 +1175 @@
 
         // DUAL solution
-        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+# if (1)
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10);
+        invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
+# endif
+# if (0)
+        psMatrixLUSolve(sumMatrix, sumVector);
+        solution = psMemIncrRefCounter(sumVector);
+        invMatrix = psMemIncrRefCounter(sumMatrix);
+# endif
 
 #if (0)
         for (int i = 0; i < solution->n; i++) {
-            fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
+            fprintf(stderr, "Dual solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i]));
         }
 #endif
 
-        if (!kernels->solution1) {
-            kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64);
-            psVectorInit (kernels->solution1, 0.0);
-        }
-        if (!kernels->solution2) {
-            kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
-            psVectorInit (kernels->solution2, 0.0);
-        }
+        // XXX TEST: manually set the coeffs to a desired solution
+        // solution->data.F64[0] = +1.826;
+        // solution->data.F64[1] = -0.115;
+        // solution->data.F64[2] =  0.0;
+        // solution->data.F64[3] =  0.0;
+
+        // ensure we have solution vectors of the right size
+        kernels->solution1    = psVectorRecycle(kernels->solution1,    numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+        kernels->solution1err = psVectorRecycle(kernels->solution1err, numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+        kernels->solution2    = psVectorRecycle(kernels->solution2,    numSolution2,     PS_TYPE_F64); // 1 for norm, 1 for bg
+        kernels->solution2err = psVectorRecycle(kernels->solution2err, numSolution2,     PS_TYPE_F64); // 1 for norm, 1 for bg
+
+        psVectorInit(kernels->solution1, 0.0);
+        psVectorInit(kernels->solution1err, 0.0);
+        psVectorInit(kernels->solution2, 0.0);
+        psVectorInit(kernels->solution2err, 0.0);
 
         // for DUAL convolution analysis, we apply the normalization to I1 as follows:
@@ -1205 +1220 @@
             kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue;
             kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
+
+            kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]) * stamps->normValue;
+            int i2 = i + numSolution1;
+            kernels->solution2err->data.F64[i] = sqrt(invMatrix->data.F64[i2][i2]);
         }
 
@@ -1213 +1232 @@
         kernels->solution1->data.F64[bgIndex] = 0.0;
 
+        psFree(solution);
+        psFree(sumVector);
         psFree(sumMatrix);
-        psFree(sumVector);
-        psFree(solution);
+        psFree(invMatrix);
     }
 
@@ -1224 +1244 @@
     if (psTraceGetLevel("psModules.imcombine") >= 7) {
         for (int i = 0; i < kernels->solution1->n; i++) {
-            psTrace("psModules.imcombine", 7, "Solution 1 %d: %f\n", i, kernels->solution1->data.F64[i]);
+            psTrace("psModules.imcombine", 7, "Solution 1 %d: %f +/- %f\n", i, kernels->solution1->data.F64[i], kernels->solution1err->data.F64[i]);
         }
         if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
             for (int i = 0; i < kernels->solution2->n; i++) {
-                psTrace("psModules.imcombine", 7, "Solution 2 %d: %f\n", i, kernels->solution2->data.F64[i]);
+                psTrace("psModules.imcombine", 7, "Solution 2 %d: %f +/- %f\n", i, kernels->solution2->data.F64[i], kernels->solution2err->data.F64[i]);
             }
         }
      }
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Solve equation: %f sec", psTimerClear("pmSubtractionSolveEquation"));
 
     // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
@@ -1283 +1305 @@
 }
 
-psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps,
-                                           pmSubtractionKernels *kernels)
+// given the convolved image(s) and the residual image, calculate the second moment(s) and the chisq
+bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window) {
+
+# ifndef USE_WEIGHT
+    psAssert(weight == NULL, "impossible!");
+# endif
+# ifndef USE_WINDOW
+    psAssert(window == NULL, "impossible!");
+# endif
+
+    int npix = 0;
+    float chisqR = 0;
+    float chisqD = 0;
+
+    // get the chisq
+    for (int y = residual->yMin; y <= residual->yMax; y++) {
+        for (int x = residual->xMin; x <= residual->xMax; x++) {
+            float valueR = PS_SQR(residual->kernel[y][x]);
+            if (weight) {
+                valueR *= weight->kernel[y][x];
+            }
+            // XXX NOTE: do NOT apply the window to the chisq portions of the calculation (that would bias the chisq)
+            chisqR += valueR;
+
+            float valueD = PS_SQR(difference->kernel[y][x]);
+            if (weight) {
+                valueD *= weight->kernel[y][x];
+            }
+            chisqD += valueD;
+            npix ++;
+        }
+    }
+    psVectorAppend(chisqRVector, chisqR / npix);
+    psVectorAppend(chisqDVector, chisqD / npix);
+
+    float value1 = 0;
+    float value2 = 0;
+    float flux2 = 0;
+    float fluxX = 0;
+    float fluxY = 0;
+    float fluxX2 = 0;
+    float fluxY2 = 0;
+
+    float fluxC1 = 0;
+    float fluxC2 = 0;
+
+    float moment = 0;
+
+    // get the moments from convolved1
+    if (convolved1) {
+        for (int y = residual->yMin; y <= residual->yMax; y++) {
+            for (int x = residual->xMin; x <= residual->xMax; x++) {
+                value1  = convolved1->kernel[y][x];
+                value2  = PS_SQR(value1);
+
+                if (window) {
+                    value1 *= window->kernel[y][x];
+                    value2 *= window->kernel[y][x];
+                }
+
+                fluxC1 += value1;
+                flux2  += value2;
+                fluxX  += x * value2;
+                fluxY  += y * value2;
+                // fluxX2 += PS_SQR(x) * value2;
+                // fluxY2 += PS_SQR(y) * value2;
+                fluxX2 += PS_SQR(x) * value1;
+                fluxY2 += PS_SQR(y) * value1;
+            }
+        }
+        // float Mx = fluxX / flux2;
+        // float My = fluxY / flux2;
+        // float Mxx = fluxX2 / flux2;
+        // float Myy = fluxY2 / flux2;
+        float Mxx = fluxX2 / fluxC1;
+        float Myy = fluxY2 / fluxC1;
+
+        // fprintf (stderr, "conv1, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix);
+        moment += Mxx + Myy;
+    }
+
+    // get the moments from convolved1
+    if (convolved2) {
+        for (int y = residual->yMin; y <= residual->yMax; y++) {
+            for (int x = residual->xMin; x <= residual->xMax; x++) {
+                value1  = convolved2->kernel[y][x];
+                value2  = PS_SQR(value1);
+
+                // XXX NOTE: do NOT apply the weight to the moments calculation
+                if (false && weight) {
+                    value2 *= weight->kernel[y][x];
+                }
+                if (window) {
+                    value1 *= window->kernel[y][x];
+                    value2 *= window->kernel[y][x];
+                }
+
+                fluxC2 += value1;
+                flux2  += value2;
+                fluxX  += x * value2;
+                fluxY  += y * value2;
+                // fluxX2 += PS_SQR(x) * value2;
+                // fluxY2 += PS_SQR(y) * value2;
+                fluxX2 += PS_SQR(x) * value1;
+                fluxY2 += PS_SQR(y) * value1;
+            }
+        }
+        // float Mx = fluxX / flux2;
+        // float My = fluxY / flux2;
+        // float Mxx = fluxX2 / flux2;
+        // float Myy = fluxY2 / flux2;
+        float Mxx = fluxX2 / fluxC2;
+        float Myy = fluxY2 / fluxC2;
+
+        // fprintf (stderr, "conv2, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix);
+        moment += Mxx + Myy;
+    }
+
+    float flux = fluxC1 + fluxC2;
+    
+    if (convolved1 && convolved2) {
+        moment *= 0.5;
+        flux *= 0.5;
+    }
+    psVectorAppend(momentVector, moment);
+    psVectorAppend(fluxesVector, flux);
+
+    // check that the last appended values are ok:
+    int Nelem = fluxesVector->n - 1;
+    bool valid = true;
+    valid &= isfinite(chisqRVector->data.F32[Nelem]);
+    valid &= isfinite(fluxesVector->data.F32[Nelem]);
+    valid &= isfinite(momentVector->data.F32[Nelem]);
+    if (valid) {
+      psVectorAppend(stampMask, 0);
+    } else {
+      psVectorAppend(stampMask, 0x02);
+    }
+    return true;
+}
+
+bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, 
+                                           pmSubtractionStampList *stamps,
+                                           pmSubtractionKernels *kernels)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
     PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
+
+    psTimerStart("pmSubtractionCalculateChisqAndMoments");
+
+    // XXX need to save these somewhere
+    psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *chisqD = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *chisqR = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *stampMask = psVectorAllocEmpty(stamps->num, PS_TYPE_VECTOR_MASK);
+
+    int footprint = stamps->footprint; // Half-size of stamps
+    int numKernels = kernels->num;      // Number of kernels
+
+    psImage *polyValues = NULL;         // Polynomial values
+
+    // storage for the image (convolved2 is not used in SINGLE mode)
+    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+    psKernel *difference = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+    psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+    psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+
+    int nGood = 0;
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            // mark this stamp as unused (note that we have to append NANs to the other vectors to keep the lengths in sync)
+            psVectorAppend(moments, NAN);
+            psVectorAppend(fluxes, NAN);
+            psVectorAppend(chisqD, NAN);
+            psVectorAppend(chisqR, NAN);
+            psVectorAppend(stampMask, 0x01);
+            continue;
+        }
+        nGood ++;
+
+        // Calculate coefficients of the kernel basis functions
+        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
+        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
+
+        // Calculate residuals
+        psImageInit(residual->image, 0.0);
+        psImageInit(difference->image, 0.0);
+
+        psKernel *weight = NULL;
+        psKernel *window = NULL;
+    
+#ifdef USE_WEIGHT
+    weight = stamp->weight;
+#endif
+#ifdef USE_WINDOW
+    window = stamps->window;
+#endif
+
+        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
+
+            // the single-direction psf match code attempts to find the kernel such that:
+            // source * kernel = target.  we need to assign 'source' and 'target' correctly
+            // depending on which of image1 or image2 we asked to be convolved.
+
+            psKernel *target;           // Target postage stamp (convolve source to match the target)
+            psKernel *source;           // Source postage stamp (convolve source to match the target)
+            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
+
+            // init the accumulation image
+            psImageInit(convolved1->image, 0.0);
+
+            switch (kernels->mode) {
+              case PM_SUBTRACTION_MODE_1:
+                target = stamp->image2;
+                source = stamp->image1;
+                convolutions = stamp->convolutions1;
+                break;
+              case PM_SUBTRACTION_MODE_2:
+                target = stamp->image1;
+                source = stamp->image2;
+                convolutions = stamp->convolutions2;
+                break;
+              default:
+                psAbort("Unsupported subtraction mode: %x", kernels->mode);
+            }
+
+            // generate the convolved source image (sum over kernels)
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *convolution = convolutions->data[j]; // Convolution
+                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        convolved1->kernel[y][x] += convolution->kernel[y][x] * coefficient;
+                    }
+                }
+            }
+
+            // Generate the difference, residual, and convolved source images.  Note the we
+            // accumulate the convolution of (A-B), so we need to replace it to generate the
+            // images of the convolved source image.
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    difference->kernel[y][x] = target->kernel[y][x] - source->kernel[y][x] * norm - background;
+                    residual->kernel[y][x] = difference->kernel[y][x] - convolved1->kernel[y][x];
+                    convolved1->kernel[y][x] += source->kernel[y][x] * norm;
+                }
+            }
+
+            // XXX if we want to have a weight and window, we'll need to pass through to here
+            pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, NULL, difference, residual, weight, window);
+
+        } else {
+
+            // Dual convolution
+            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
+            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
+            psKernel *image1 = stamp->image1; // The first image
+            psKernel *image2 = stamp->image2; // The second image
+
+            // init the accumulation images
+            psImageInit(convolved1->image, 0.0);
+            psImageInit(convolved2->image, 0.0);
+
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
+                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
+                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
+                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
+
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        // NOTE sign for coeff2
+                        convolved1->kernel[y][x] += +conv1->kernel[y][x] * coeff1;
+                        convolved2->kernel[y][x] += -conv2->kernel[y][x] * coeff2;
+                    }
+                }
+            }
+
+            // Generate the difference, residual, and convolved source images.  Note the we
+            // accumulate the convolutions of (A-B), so we need to replace (A or B) to generate
+            // the images of the convolved source images.
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    difference->kernel[y][x] = image2->kernel[y][x] - image1->kernel[y][x] * norm - background;
+                    residual->kernel[y][x] = difference->kernel[y][x] + convolved2->kernel[y][x] - convolved1->kernel[y][x];
+                    convolved1->kernel[y][x] += image1->kernel[y][x] * norm;
+                    convolved2->kernel[y][x] += image2->kernel[y][x];
+                }
+            }
+
+            if (0) {
+                psFitsWriteImageSimple("conv1.fits", convolved1->image, NULL);
+                psFitsWriteImageSimple("conv2.fits", convolved2->image, NULL);
+                psFitsWriteImageSimple("resid.fits", residual->image,   NULL);
+                pmVisualAskUser(NULL);
+            } 
+
+            pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, convolved2, difference, residual, weight, window);
+        }
+    }
+
+    // find the mean chisq and mean moment
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
+    psVectorStats (stats, chisqD, NULL, stampMask, 0xff);
+    float chisqDValue = stats->sampleMean;
+
+    psStatsInit(stats);
+    psVectorStats (stats, chisqR, NULL, stampMask, 0xff);
+    float chisqRValue = stats->sampleMean;
+
+    psStatsInit(stats);
+    psVectorStats (stats, moments, NULL, stampMask, 0xff);
+    float momentValue = stats->sampleMean;
+
+    double sumKernel1 = 0.0, sumKernel2 = 0.0; // Sum of the kernel
+
+    // calculate the variance contribution from this smoothing kernel
+    psKernel *modelKernel = pmSubtractionKernel(kernels, 0.0, 0.0, false);
+    for (int y = modelKernel->yMin; y <= modelKernel->yMax; y++) {
+        for (int x = modelKernel->xMin; x <= modelKernel->xMax; x++) {
+            if (!isfinite(modelKernel->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Non-finite covariance matrix element in kernel at %d,%d", x, y);
+                return NULL;
+            }
+            sumKernel1 += PS_SQR(modelKernel->kernel[y][x]);
+        }
+    }
+    psFree (modelKernel);
+
+    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        psKernel *modelKernel = pmSubtractionKernel(kernels, 0.0, 0.0, true);
+        for (int y = modelKernel->yMin; y <= modelKernel->yMax; y++) {
+            for (int x = modelKernel->xMin; x <= modelKernel->xMax; x++) {
+                if (!isfinite(modelKernel->kernel[y][x])) {
+                    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Non-finite covariance matrix element in kernel at %d,%d", x, y);
+                    return NULL;
+                }
+                sumKernel2 += PS_SQR(modelKernel->kernel[y][x]);
+            }
+        }
+        psFree (modelKernel);
+    } else {
+        sumKernel2 = 1.0;
+    }
+
+    // if we modify the chisq value by the (sumKernel1 + sumKernel2), we account for the
+    // smoothing coming from larger kernels adding additional spatial fit terms should be
+    // penalized by increasing the score somewhat.  the 0.01 value is not well-chosen.
+    float orderFactor = 0.01 * kernels->spatialOrder;
+    float score = 2.0 * chisqRValue / (sumKernel1 + sumKernel2) + orderFactor;
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, chisqD: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqRValue, chisqDValue, momentValue, sumKernel1, sumKernel2, score);
+
+    // save this result if it is the first or the best (skip if bestMatch is NULL)
+    if (bestMatch) {
+        pmSubtractionQuality *match = *bestMatch;
+        bool keep = false;
+        if (match == NULL) {
+            *bestMatch = match = pmSubtractionQualityAlloc();
+            keep = true;
+        } else {
+            if (score < match->score) {
+                psFree(match->fluxes);
+                psFree(match->chisq);
+                psFree(match->moments);
+                psFree(match->stampMask);
+                keep = true;
+            }
+        }
+        if (keep) {
+            psLogMsg("psModules.imcombine", PS_LOG_INFO, "keeping order: %d, mode: %d, score: %f\n", kernels->spatialOrder, kernels->mode, score);
+            match->score        = score;
+            match->spatialOrder = kernels->spatialOrder;
+            match->mode         = kernels->mode;
+            match->nGood        = nGood;
+            match->fluxes       = psMemIncrRefCounter(fluxes);
+            match->chisq        = psMemIncrRefCounter(chisqR);
+            match->moments      = psMemIncrRefCounter(moments);
+            match->stampMask    = psMemIncrRefCounter(stampMask);
+        }           
+    }
+
+    pmSubtractionVisualPlotChisqAndMoments(fluxes, chisqR, moments);
+
+    psFree(stats);
+    psFree(chisqR);
+    psFree(chisqD);
+    psFree(fluxes);
+    psFree(moments);
+    psFree(stampMask);
+
+    psFree(residual);
+    psFree(difference);
+    psFree(convolved1);
+    psFree(convolved2);
+    psFree(polyValues);
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Chisq and Moments: %f sec", psTimerClear("pmSubtractionCalculateChisqAndMoments"));
+
+    return true;
+}
+
+// XXX for now, let's not use this, and let's instead just use values from pmSubtractionCalculateChisqAndMoments
+psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
+    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
+
+    psTimerStart("pmSubtractionCalculateDeviations");
 
     psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps
@@ -1298 +1727 @@
     psImage *polyValues = NULL;         // Polynomial values
     psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
-
-    // set up holding images for the visualization
-    pmSubtractionVisualShowFitInit (stamps);
 
     psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
@@ -1401 +1827 @@
             }
 
-            // XXX visualize the target, source, convolution and residual
-            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
-
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
@@ -1438 +1861 @@
             }
 
-            // XXX visualize the target, source, convolution and residual
-            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
-
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
@@ -1455 +1875 @@
         }
 
+        double flux = 0.0;
         double deviation = 0.0;         // Sum of differences
         for (int y = - footprint; y <= footprint; y++) {
@@ -1460 +1881 @@
                 double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
                 deviation += dev;
-#ifdef TESTING
-                residual->kernel[y][x] = dev;
-#endif
+                flux += stamp->image1->kernel[y][x] + stamp->image2->kernel[y][x];
             }
         }
         deviations->data.F32[i] = devNorm * deviation;
-        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
-                i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]);
+        psTrace("psModules.imcombine", 5, "Deviation and Flux for stamp %d (%d,%d): %f %f\n",
+                i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i], flux);
         psStringAppend(&log, "Stamp %d (%d,%d): %f\n",
                        i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]);
         if (!isfinite(deviations->data.F32[i])) {
             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-            psTrace("psModules.imcombine", 5,
-                    "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
-                    i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
+            psTrace("psModules.imcombine", 5, "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
             continue;
         }
-
-#ifdef TESTING
-        {
-            psString filename = NULL;
-            psStringAppend(&filename, "resid_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (stamp->image1) {
-            psString filename = NULL;
-            psStringAppend(&filename, "stamp_image1_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (stamp->image2) {
-            psString filename = NULL;
-            psStringAppend(&filename, "stamp_image2_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (stamp->weight) {
-            psString filename = NULL;
-            psStringAppend(&filename, "stamp_weight_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
-            psFitsClose(fits);
-        }
-#endif
-
     }
 
     psFree(keepStamps);
 
-    psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log);
+    psLogMsg("psModules.imcombine", PS_LOG_MINUTIA, "%s", log);
     psFree(log);
 
@@ -1527 +1908 @@
         psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
 
-        pmSubtractionVisualShowFit(norm);
-        pmSubtractionVisualPlotFit(kernels);
-
         psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
         psVectorStats (stats, fResSigma, NULL, NULL, 0);
@@ -1560 +1938 @@
     psFree(polyValues);
 
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Deviations: %f sec", psTimerClear("pmSubtractionCalculateDeviations"));
+
     return deviations;
-}
-
-// we are supplied U, not Ut; w represents a diagonal matrix (also, we apply 1/w instead of w)
-psImage *p_pmSubSolve_wUt (psVector *w, psImage *U) {
-
-    psAssert (w->n == U->numCols, "w and U dimensions do not match");
-
-    // wUt has dimensions transposed relative to Ut.
-    psImage *wUt = psImageAlloc (U->numRows, U->numCols, PS_TYPE_F64);
-    psImageInit (wUt, 0.0);
-
-    for (int i = 0; i < wUt->numCols; i++) {
-        for (int j = 0; j < wUt->numRows; j++) {
-            if (!isfinite(w->data.F64[j])) continue;
-            if (w->data.F64[j] == 0.0) continue;
-            wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
-        }
-    }
-    return wUt;
-}
-
-// XXX this is just standard matrix multiplication: use psMatrixMultiply?
-psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt) {
-
-    psAssert (V->numCols == wUt->numRows, "matrix dimensions do not match");
-
-    psImage *Ainv = psImageAlloc (wUt->numCols, V->numRows, PS_TYPE_F64);
-
-    for (int i = 0; i < Ainv->numCols; i++) {
-        for (int j = 0; j < Ainv->numRows; j++) {
-            double sum = 0.0;
-            for (int k = 0; k < V->numCols; k++) {
-                sum += V->data.F64[j][k] * wUt->data.F64[k][i];
-            }
-            Ainv->data.F64[j][i] = sum;
-        }
-    }
-    return Ainv;
-}
-
-// we are supplied U, not Ut
-bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B) {
-
-    psAssert (U->numRows == B->n, "U and B dimensions do not match");
-
-    UtB[0] = psVectorRecycle (UtB[0], U->numCols, PS_TYPE_F64);
-
-    for (int i = 0; i < U->numCols; i++) {
-        double sum = 0.0;
-        for (int j = 0; j < U->numRows; j++) {
-            sum += B->data.F64[j] * U->data.F64[j][i];
-        }
-        UtB[0]->data.F64[i] = sum;
-    }
-    return true;
-}
-
-// w is diagonal
-bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB) {
-
-    psAssert (w->n == UtB->n, "w and UtB dimensions do not match");
-
-    // wUt has dimensions transposed relative to Ut.
-    wUtB[0] = psVectorRecycle (wUtB[0], w->n, PS_TYPE_F64);
-    psVectorInit (wUtB[0], 0.0);
-
-    for (int i = 0; i < w->n; i++) {
-        if (!isfinite(w->data.F64[i])) continue;
-        if (w->data.F64[i] == 0.0) continue;
-        wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
-    }
-    return true;
-}
-
-// this is basically matrix * vector
-bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB) {
-
-    psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match");
-
-    VwUtB[0] = psVectorRecycle (*VwUtB, V->numRows, PS_TYPE_F64);
-
-    for (int j = 0; j < V->numRows; j++) {
-        double sum = 0.0;
-        for (int i = 0; i < V->numCols; i++) {
-            sum += V->data.F64[j][i] * wUtB->data.F64[i];
-        }
-        VwUtB[0]->data.F64[j] = sum;
-    }
-    return true;
-}
-
-// this is basically matrix * vector
-bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x) {
-
-    psAssert (A->numCols == x->n, "A and x dimensions do not match");
-
-    B[0] = psVectorRecycle (*B, A->numRows, PS_TYPE_F64);
-
-    for (int j = 0; j < A->numRows; j++) {
-        double sum = 0.0;
-        for (int i = 0; i < A->numCols; i++) {
-            sum += A->data.F64[j][i] * x->data.F64[i];
-        }
-        B[0]->data.F64[j] = sum;
-    }
-    return true;
-}
-
-// this is basically Vector * vector
-bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y) {
-
-    psAssert (x->n == y->n, "x and y dimensions do not match");
-
-    double sum = 0.0;
-    for (int i = 0; i < x->n; i++) {
-        sum += x->data.F64[i] * y->data.F64[i];
-    }
-    *value = sum;
-    return true;
-}
-
-bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
-
-    int footprint = stamps->footprint; // Half-size of stamps
-
-    double sum = 0.0;
-    for (int i = 0; i < stamps->num; i++) {
-
-        pmSubtractionStamp *stamp = stamps->stamps->data[i];
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
-
-        psKernel *weight = NULL;
-        psKernel *window = NULL;
-        psKernel *input = NULL;
-
-#ifdef USE_WEIGHT
-        weight = stamp->weight;
-#endif
-#ifdef USE_WINDOW
-        window = stamps->window;
-#endif
-
-        switch (kernels->mode) {
-            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
-          case PM_SUBTRACTION_MODE_1:
-            input = stamp->image2;
-            break;
-          case PM_SUBTRACTION_MODE_2:
-            input = stamp->image1;
-            break;
-          default:
-            psAbort ("programming error");
-        }
-
-        for (int y = - footprint; y <= footprint; y++) {
-            for (int x = - footprint; x <= footprint; x++) {
-                double in = input->kernel[y][x];
-                double value = in*in;
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    value *= wtVal;
-                }
-                if (window) {
-                    float  winVal = window->kernel[y][x];
-                    value *= winVal;
-                }
-                sum += value;
-            }
-        }
-    }
-    *y2 = sum;
-    return true;
-}
-
-double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
-
-    int footprint = stamps->footprint; // Half-size of stamps
-    int numKernels = kernels->num;      // Number of kernels
-
-    double sum = 0.0;
-
-    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
-    psImageInit(residual->image, 0.0);
-
-    psImage *polyValues = NULL;         // Polynomial values
-
-    for (int i = 0; i < stamps->num; i++) {
-
-        pmSubtractionStamp *stamp = stamps->stamps->data[i];
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
-
-        psKernel *weight = NULL;
-        psKernel *window = NULL;
-        psKernel *target = NULL;
-        psKernel *source = NULL;
-
-        psArray *convolutions = NULL;
-
-#ifdef USE_WEIGHT
-        weight = stamp->weight;
-#endif
-#ifdef USE_WINDOW
-        window = stamps->window;
-#endif
-
-        switch (kernels->mode) {
-            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
-          case PM_SUBTRACTION_MODE_1:
-            target = stamp->image2;
-            source = stamp->image1;
-            convolutions = stamp->convolutions1;
-            break;
-          case PM_SUBTRACTION_MODE_2:
-            target = stamp->image1;
-            source = stamp->image2;
-            convolutions = stamp->convolutions2;
-            break;
-          default:
-            psAbort ("programming error");
-        }
-
-        // Calculate coefficients of the kernel basis functions
-        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
-        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
-        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
-
-        psImageInit(residual->image, 0.0);
-        for (int j = 0; j < numKernels; j++) {
-            psKernel *convolution = convolutions->data[j]; // Convolution
-            double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
-                }
-            }
-        }
-
-        for (int y = - footprint; y <= footprint; y++) {
-            for (int x = - footprint; x <= footprint; x++) {
-                double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
-                double value = PS_SQR(resid);
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    value *= wtVal;
-                }
-                if (window) {
-                    float  winVal = window->kernel[y][x];
-                    value *= winVal;
-                }
-                sum += value;
-            }
-        }
-    }
-    psFree (polyValues);
-    psFree (residual);
-
-    return sum;
-}
-
-bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask) {
-
-    for (int i = 0; i < w->n; i++) {
-        wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
-    }
-    return true;
-}
-
-// we are supplied V and w; w represents a diagonal matrix (also, we apply 1/w instead of w)
-psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w) {
-
-    psAssert (w->n == V->numCols, "w and U dimensions do not match");
-
-    psImage *Vn = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
-    psImageInit (Vn, 0.0);
-
-    // generate Vn = V * w^{-1}
-    for (int j = 0; j < Vn->numRows; j++) {
-        for (int i = 0; i < Vn->numCols; i++) {
-            if (!isfinite(w->data.F64[i])) continue;
-            if (w->data.F64[i] == 0.0) continue;
-            Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
-        }
-    }
-
-    psImage *Xvar = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
-    psImageInit (Xvar, 0.0);
-
-    // generate Xvar = Vn * Vn^T
-    for (int j = 0; j < Vn->numRows; j++) {
-        for (int i = 0; i < Vn->numCols; i++) {
-            double sum = 0.0;
-            for (int k = 0; k < Vn->numCols; k++) {
-                sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
-            }
-            Xvar->data.F64[j][i] = sum;
-        }
-    }
-    return Xvar;
 }
 
@@ -1864 +1946 @@
 // of the elements of an image A(x,y) = A->data.F64[y][x] = A_x,y, a matrix
 // multiplication is: A_k,j * B_i,k = C_i,j
-
 
 bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header) {

trunk/psModules/src/imcombine/pmSubtractionEquation.h

@@ -4 +4 @@
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionKernels.h"
-
-typedef enum {
-    PM_SUBTRACTION_EQUATION_NONE    = 0x00,
-    PM_SUBTRACTION_EQUATION_NORM    = 0x01,
-    PM_SUBTRACTION_EQUATION_BG      = 0x02,
-    PM_SUBTRACTION_EQUATION_KERNELS = 0x04,
-    PM_SUBTRACTION_EQUATION_ALL     = 0x07, // value should be NORM | BG | KERNELS
-} pmSubtractionEquationCalculationMode;
+#include "pmSubtraction.h"
 
 /// Execute a thread job to calculate the least-squares equation for a stamp
@@ -20 +13 @@
 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps
                                          pmSubtractionKernels *kernels, ///< Kernel parameters
-                                         int index, ///< Index of stamp
-                                         const pmSubtractionEquationCalculationMode mode
+                                         int index ///< Index of stamp
     );
 
 /// Calculate the least-squares equation to match the image quality
 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
-                                    pmSubtractionKernels *kernels, ///< Kernel parameters
-                                    const pmSubtractionEquationCalculationMode mode
+                                    pmSubtractionKernels *kernels ///< Kernel parameters
     );
 
 /// Solve the least-squares equation to match the image quality
 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters
-                                const pmSubtractionStampList *stamps, ///< Stamps
-                                const pmSubtractionEquationCalculationMode mode
+                                const pmSubtractionStampList *stamps ///< Stamps
     );
 
@@ -92 +82 @@
 bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window);
 
+bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window);
+
+bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, pmSubtractionStampList *stamps, pmSubtractionKernels *kernels);
 #endif

trunk/psModules/src/imcombine/pmSubtractionEquation.v0.c

@@ -9 +9 @@
 #include "pmErrorCodes.h"
 #include "pmSubtraction.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtractionKernels.h"
 #include "pmSubtractionStamps.h"

trunk/psModules/src/imcombine/pmSubtractionHermitian.c

@@ -8 +8 @@
 #include <pslib.h>
 
+#include "pmSubtractionTypes.h"
 #include "pmSubtractionHermitian.h"
 

trunk/psModules/src/imcombine/pmSubtractionIO.c

@@ -12 +12 @@
 #include "pmConceptsRead.h"
 
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -8 +8 @@
 #include <pslib.h>
 
+#include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"
@@ -20 +22 @@
 {
     psFree(kernels->description);
+    psFree(kernels->fwhms);
+    psFree(kernels->orders);
     psFree(kernels->u);
     psFree(kernels->v);
@@ -30 +34 @@
     psFree(kernels->solution1);
     psFree(kernels->solution2);
+    psFree(kernels->solution1err);
+    psFree(kernels->solution2err);
     psFree(kernels->sampleStamps);
 }
@@ -419 +425 @@
 
     int num = 0;                        // Number of basis functions
-    psString params = NULL;             // List of parameters
     for (int i = 0; i < numGaussians; i++) {
         int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
-        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
         num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
     }
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size,
-                                                              spatialOrder, penalty, bounds, mode); // Kernels
-    psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements",
-             params, spatialOrder, num);
-    psFree(params);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels
+    pmSubtractionKernelsMakeDescription(kernels);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
 
 # if (!CENTRAL_DELTA && !ZERO_KERNEL_ZERO_FLUX)
@@ -503 +503 @@
 
     int num = 0;                        // Number of basis functions
-    psString params = NULL;             // List of parameters
     for (int i = 0; i < numGaussians; i++) {
         int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
-        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
         num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
         num += (11 - gaussOrder - 1);   // include all higher order radial terms
     }
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size,
-                                                              spatialOrder, penalty, bounds, mode); // Kernels
-    psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS_RADIAL kernel: %s,%d --> %d elements", params, spatialOrder, num);
-    psFree(params);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels
+    pmSubtractionKernelsMakeDescription(kernels);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
 
     // Set the kernel parameters
@@ -569 +564 @@
 
     int num = 0;                        // Number of basis functions
-    psString params = NULL;             // List of parameters
     for (int i = 0; i < numGaussians; i++) {
         int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
-        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
         num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
     }
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size,
-                                                              spatialOrder, penalty, bounds, mode); // Kernels
-    psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements",
-             params, spatialOrder, num);
-    psFree(params);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels
+    pmSubtractionKernelsMakeDescription(kernels);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
 
     // Set the kernel parameters
@@ -627 +616 @@
 
     int num = 0;                        // Number of basis functions
-    psString params = NULL;             // List of parameters
     for (int i = 0; i < numGaussians; i++) {
         int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
-        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
         num += PS_SQR(gaussOrder + 1);
     }
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size,
-                                                              spatialOrder, penalty, bounds, mode); // Kernels
-    psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "DECONVOLVED HERM kernel: %s,%d --> %d elements", params, spatialOrder, num);
-    psFree(params);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels
+    pmSubtractionKernelsMakeDescription(kernels);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
 
     // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest)
@@ -713 +697 @@
 
 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
-                                                int size, int spatialOrder, float penalty, psRegion bounds,
+                                                int size, psVector *fwhms, psVector *orders, int spatialOrder, float penalty, psRegion bounds,
                                                 pmSubtractionMode mode)
 {
@@ -722 +706 @@
     kernels->description = NULL;
     kernels->num = numBasisFunctions;
+    kernels->fwhms = psMemIncrRefCounter(fwhms);
+    kernels->orders = psMemIncrRefCounter(orders);
     kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
     kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
@@ -740 +726 @@
     kernels->size = size;
     kernels->inner = 0;
+    kernels->binning = 0;
+    kernels->ringsOrder = 0;
     kernels->spatialOrder = spatialOrder;
     kernels->bgOrder = 0;
@@ -745 +733 @@
     kernels->solution1 = NULL;
     kernels->solution2 = NULL;
+    kernels->solution1err = NULL;
+    kernels->solution2err = NULL;
     kernels->mean = NAN;
     kernels->rms = NAN;
@@ -823 +813 @@
     int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size,
-                                                              spatialOrder, penalty, bounds, mode); // Kernels
-    psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty);
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
-             size, spatialOrder, num);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels
+    pmSubtractionKernelsMakeDescription(kernels);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
 
     if (!p_pmSubtractionKernelsAddGrid(kernels, 0, size)) {
@@ -871 +859 @@
     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size,
-                                                              spatialOrder, penalty, bounds, mode); // Kernels
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels
     kernels->inner = inner;
-    psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder,
-                   penalty);
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements",
-             size, inner, binning, spatialOrder, num);
+    kernels->binning = binning;
+    pmSubtractionKernelsMakeDescription(kernels);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
 
     kernels->uStop = psVectorAlloc(num, PS_TYPE_S32);
@@ -970 +955 @@
     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size,
-                                                              spatialOrder, penalty, bounds, mode); // Kernels
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels
     kernels->inner = inner;
-    psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty);
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements",
-             size, inner, spatialOrder, num);
+    pmSubtractionKernelsMakeDescription(kernels);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
 
     kernels->uStop = psVectorAlloc(num, PS_TYPE_S32);
@@ -1053 +1035 @@
     PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
 
-    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
-                                                                  penalty, bounds, mode); // Kernels
+    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode); // Kernels
     kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
-    psStringPrepend(&kernels->description, "GUNK=");
-    psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);
+    kernels->inner = inner;
+    pmSubtractionKernelsMakeDescription(kernels);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, (int) kernels->num);
 
     int numISIS = kernels->num;         // Number of ISIS kernels
@@ -1100 +1082 @@
     int num = numRings * numPoly; // Total number of basis functions
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size,
-                                                              spatialOrder, penalty, bounds, mode); // Kernels
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels
     kernels->inner = inner;
-    psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder,
-                   penalty);
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements",
-             size, inner, ringsOrder, spatialOrder, num);
+    kernels->ringsOrder = ringsOrder;
+    pmSubtractionKernelsMakeDescription(kernels);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
 
     // Set the Gaussian kernel parameters
@@ -1405 +1384 @@
 }
 
+bool pmSubtractionKernelsMakeDescription(pmSubtractionKernels *kernels) {
+
+    // free if it exists
+    psFree (kernels->description);
+
+    // generate the description parameter string
+    psString params = NULL;
+    if (kernels->fwhms) {
+        for (int i = 0; i < kernels->fwhms->n; i++) {
+            psStringAppend(&params, "(%.1f,%d)", kernels->fwhms->data.F32[i], kernels->orders->data.S32[i]);
+        }
+    }
+
+    switch (kernels->type) {
+      case PM_SUBTRACTION_KERNEL_ISIS:
+        psStringAppend (&kernels->description, "ISIS(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
+        break;
+
+      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
+        psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
+        break;
+
+      case PM_SUBTRACTION_KERNEL_HERM:
+        psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
+        break;
+
+      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
+        psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
+        break;
+
+      case PM_SUBTRACTION_KERNEL_POIS:
+        psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", kernels->size, kernels->spatialOrder, kernels->penalty);
+        break;
+
+      case PM_SUBTRACTION_KERNEL_SPAM:
+        psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", kernels->size, kernels->inner, kernels->binning, kernels->spatialOrder, kernels->penalty);
+        break;
+
+      case PM_SUBTRACTION_KERNEL_FRIES:
+        psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", kernels->size, kernels->inner, kernels->spatialOrder, kernels->penalty);
+        break;
+
+        // Grid United with Normal Kernel [description: GUNK=ISIS(...)+POIS(...)]
+      case PM_SUBTRACTION_KERNEL_GUNK:
+        psStringAppend(&kernels->description, "GUNK=ISIS(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
+        psStringAppend(&kernels->description, "+POIS(%d,%d)", kernels->inner, kernels->spatialOrder);
+        break;
+        
+      case PM_SUBTRACTION_KERNEL_RINGS:
+        psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", kernels->size, kernels->inner, kernels->ringsOrder, kernels->spatialOrder, kernels->penalty);
+        break;
+
+      default:
+        psAbort("unknown kernel");
+    }
+    psFree (params);
+    return true;
+}
+
 pmSubtractionKernels *pmSubtractionKernelsCopy(const pmSubtractionKernels *in)
 {
@@ -1435 +1473 @@
     out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
     out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;
+    out->solution1err = in->solution1err ? psVectorCopy(NULL, in->solution1err, PS_TYPE_F64) : NULL;
+    out->solution2err = in->solution2err ? psVectorCopy(NULL, in->solution2err, PS_TYPE_F64) : NULL;
     out->sampleStamps = psMemIncrRefCounter(in->sampleStamps);
 

trunk/psModules/src/imcombine/pmSubtractionKernels.h

@@ -2 +2 @@
 #define PM_SUBTRACTION_KERNELS_H
 
-#include <string.h>
-#include <pslib.h>
-
-/// Type of subtraction kernel
-typedef enum {
-    PM_SUBTRACTION_KERNEL_NONE,         ///< Nothing --- an error
-    PM_SUBTRACTION_KERNEL_POIS,         ///< Pan-STARRS Optimal Image Subtraction --- delta functions
-    PM_SUBTRACTION_KERNEL_ISIS,         ///< Traditional kernel --- gaussians modified by polynomials
-    PM_SUBTRACTION_KERNEL_ISIS_RADIAL,  ///< ISIS + higher-order radial Hermitians
-    PM_SUBTRACTION_KERNEL_HERM,         ///< Hermitian polynomial kernels
-    PM_SUBTRACTION_KERNEL_DECONV_HERM,  ///< Deconvolved Hermitian polynomial kernels
-    PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
-    PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
-    PM_SUBTRACTION_KERNEL_GUNK,         ///< Grid United with Normal Kernel --- POIS and ISIS hybrid
-    PM_SUBTRACTION_KERNEL_RINGS,        ///< Rings Instead of the Normal Gaussian Subtraction
-} pmSubtractionKernelsType;
-
-/// Modes --- specifies which image to convolve
-typedef enum {
-    PM_SUBTRACTION_MODE_ERR,            // Error in the mode
-    PM_SUBTRACTION_MODE_1,              // Convolve image 1
-    PM_SUBTRACTION_MODE_2,              // Convolve image 2
-    PM_SUBTRACTION_MODE_UNSURE,         // Not sure yet which image to convolve so try to satisfy both
-    PM_SUBTRACTION_MODE_DUAL,           // Dual convolution
-} pmSubtractionMode;
-
-/// Kernels specification
-typedef struct {
-    pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
-    psString description;               ///< Description of the kernel parameters
-    int xMin, xMax, yMin, yMax;         ///< Bounds of image (for normalisation)
-    long num;                           ///< Number of kernel components (not including the spatial ones)
-    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM)
-    psVector *widths;                   ///< Gaussian FWHMs (ISIS, HERM or DECONV_HERM)
-    psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
-    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM)
-    float penalty;                      ///< Penalty for wideness
-    psVector *penalties1;               ///< Penalty for each kernel component
-    psVector *penalties2;               ///< Penalty for each kernel component
-    bool havePenalties;                 ///< flag to test if we have already calculated the penalties or not.
-    int size;                           ///< The half-size of the kernel
-    int inner;                          ///< The size of an inner region
-    int spatialOrder;                   ///< The spatial order of the kernels
-    int bgOrder;                        ///< The order for the background fitting
-    pmSubtractionMode mode;             ///< Mode for subtraction
-    psVector *solution1, *solution2;    ///< Solution for the PSF matching
-    // Quality information
-    float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
-    int numStamps;                      ///< Number of good stamps
-    float fResSigmaMean;                ///< mean fractional stdev of residuals
-    float fResSigmaStdev;               ///< stdev of fractional stdev of residuals
-    float fResOuterMean;                ///< mean fractional positive swing in residuals
-    float fResOuterStdev;               ///< stdev of fractional positive swing in residuals
-    float fResTotalMean;                ///< mean fractional negative swing in residuals
-    float fResTotalStdev;               ///< stdev of fractional negative swing in residuals
-    psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
-} pmSubtractionKernels;
-
-// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures
-typedef struct {
-    psVector *uCoords;                  // used by RINGS
-    psVector *vCoords;                  // used by RINGS
-    psVector *poly;                     // used by RINGS
-
-    psVector *xKernel;                  // used by ISIS, HERM, DECONV_HERM
-    psVector *yKernel;                  // used by ISIS, HERM, DECONV_HERM
-    psKernel *kernel;                   // used by ISIS, HERM, DECONV_HERM
-} pmSubtractionKernelPreCalc;
+// #include <string.h>
+// #include <pslib.h>
 
 // Assertion to check pmSubtractionKernels
@@ -162 +96 @@
                                                 pmSubtractionKernelsType type, ///< Kernel type
                                                 int size, ///< Half-size of kernel
+                                                psVector *fwhms, ///< requested kernel basis function
+                                                psVector *orders,
                                                 int spatialOrder, ///< Order of spatial variations
                                                 float penalty, ///< Penalty for wideness
@@ -303 +239 @@
     );
 
+bool pmSubtractionKernelsMakeDescription(pmSubtractionKernels *kernels);
+
+
 /// Copy kernels
 ///

trunk/psModules/src/imcombine/pmSubtractionMask.c

@@ -7 +7 @@
 
 #include "pmErrorCodes.h"
+#include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -11 +11 @@
 #include "pmFPA.h"
 #include "pmHDUUtils.h"
+#include "pmSubtractionTypes.h"
+#include "pmSubtraction.h"
 #include "pmSubtractionParams.h"
 #include "pmSubtractionKernels.h"
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionEquation.h"
-#include "pmSubtraction.h"
 #include "pmSubtractionAnalysis.h"
 #include "pmSubtractionMask.h"
@@ -27 +28 @@
 
 static bool useFFT = true;              // Do convolutions using FFT
-
-# define SUBMODE PM_SUBTRACTION_EQUATION_ALL
 
 //#define TESTING
@@ -59 +58 @@
 }
 
-
-static bool subtractionGetStamps(pmSubtractionStampList **stamps, // Stamps to read
-                                 const pmReadout *ro1, // Readout 1
-                                 const pmReadout *ro2, // Readout 2
-                                 const psImage *subMask, // Mask for subtraction, or NULL
-                                 psImage *variance,  // Variance map
-                                 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
-                                 pmSubtractionMode mode // Mode for subtraction
-    )
-{
-    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");
-
-    psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest
-
-    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
-                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
-    if (!*stamps) {
-        psError(psErrorCodeLast(), false, "Unable to find stamps.");
-        return false;
-    }
-
-    memCheck("  find stamps");
-
-    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
-        psError(psErrorCodeLast(), false, "Unable to extract stamps.");
-        return false;
-    }
-
-    memCheck("   extract stamps");
-    pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1);
-    return true;
-}
-
 // Check input arguments
 static bool subtractionMatchCheck(pmReadout *conv1, pmReadout *conv2, // Convolved images
@@ -123 +73 @@
                                   float badFrac,   // Maximum fraction of bad input pixels to accept
                                   pmSubtractionMode subMode // Mode of subtraction
-                                  )
+    )
 {
     if (subMode != PM_SUBTRACTION_MODE_2) {
@@ -473 +423 @@
 }
 
+bool pmSubtractionMatchAttempt(pmSubtractionQuality **bestMatch, pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, pmSubtractionMode mode, int spatialOrder, bool final) {
+
+    pmSubtractionMode nativeMode = kernels->mode;
+    pmSubtractionMode nativeOrder = kernels->spatialOrder;
+
+    kernels->mode = mode;
+    kernels->spatialOrder = spatialOrder;
+
+    // we always need to recalculate the matrix equation elements...
+    pmSubtractionStampsResetStatus(stamps);
+
+    psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n");
+    if (!pmSubtractionConvolveStamps(stamps, kernels)) {
+        psError(psErrorCodeLast(), false, "Unable to convolve stamps.");
+        return false;
+    }
+
+    // step 1: generate the elements of the matrix equation Ax = B
+    psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");
+    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
+        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+                    
+    // step 2: solve the matrix equation Ax = B
+    psTrace("psModules.imcombine", 3, "Solving kernel equations...\n");
+    if (!pmSubtractionSolveEquation(kernels, stamps)) {
+        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+    memCheck("  solve equation");
+
+    // calculate the score for this model fit attempt
+    // XXX store the chisq, flux and moments for stamp rejection
+    pmSubtractionCalculateChisqAndMoments(bestMatch, stamps, kernels); // Stamp deviations
+
+    // display the input and model stamps
+    pmSubtractionVisualShowFit(stamps, kernels);
+    pmSubtractionVisualPlotFit(kernels);
+    pmSubtractionVisualPlotConvKernels(kernels);
+
+    // reset the kernel if desired (on final pass, do not reset)
+    if (!final) {
+        kernels->mode = nativeMode;
+        kernels->spatialOrder = nativeOrder;
+    } else {
+      pmSubtractionKernelsMakeDescription(kernels);
+      psLogMsg("psModules.imcombine", PS_LOG_INFO, "final kernel: %s", kernels->description);
+    }
+    return true;
+}
 
 bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
@@ -478 +479 @@
                         const psArray *sources, const char *stampsName,
                         pmSubtractionKernelsType type, int size, int spatialOrder,
-                        const psVector *isisWidths, const psVector *isisOrders,
+                        psVector *isisWidths, const psVector *isisOrders,
                         int inner, int ringsOrder, int binning, float penalty,
                         bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
@@ -559 +560 @@
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
 
+    pmSubtractionQuality *bestMatch = NULL;
+
+    int N_TEST_MODES;
+    int N_TEST_ORDER = spatialOrder;
+
+    pmSubtractionMode TestModes[3];
+    switch (subMode) {
+      case PM_SUBTRACTION_MODE_1:
+        N_TEST_MODES = 1;
+        TestModes[0] = PM_SUBTRACTION_MODE_1;
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        N_TEST_MODES = 1;
+        TestModes[0] = PM_SUBTRACTION_MODE_2;
+        break;
+      case PM_SUBTRACTION_MODE_DUAL:
+        N_TEST_MODES = 3;
+        TestModes[0] = PM_SUBTRACTION_MODE_1;
+        TestModes[1] = PM_SUBTRACTION_MODE_2;
+        TestModes[2] = PM_SUBTRACTION_MODE_DUAL;
+        break;
+      default:
+        psError(psErrorCodeLast(), false, "For now, only modes 1, 2, and DUAL are supported.");
+        goto MATCH_ERROR;
+    }
+    
     memCheck("start");
 
@@ -628 +655 @@
             regionString = psRegionToString(*region);
             psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n",
-                    regionString, numCols, numRows);
+                     regionString, numCols, numRows);
 
             if (stampsName && strlen(stampsName) > 0) {
@@ -640 +667 @@
             }
 
-            // 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, 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);
-            }
-
-            memCheck("kernels");
-
-            if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
-                pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej);
-                switch (newMode) {
-                  case PM_SUBTRACTION_MODE_1:
-                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2.");
-                    break;
-                  case PM_SUBTRACTION_MODE_2:
-                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1.");
-                    break;
-                  default:
-                    psError(psErrorCodeLast(), false, "Unable to determine subtraction order.");
-                    goto MATCH_ERROR;
-                }
-                subMode = newMode;
-            }
-
-            int numRejected = -1;       // Number of rejected stamps in each iteration
-            for (int k = 0; k < iter && numRejected != 0; k++) {
-                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
-
-                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 stamps window.");
-                    goto MATCH_ERROR;
-                }
+            bool tryAgain = true;
+            while (tryAgain) {
+                // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
+                // doesn't matter.
+                if (!pmSubtractionStampsSelect(&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(&tryAgain, stamps, size)) {
+                    // if we failed, it might be due to the desired normWindow being larger than the current footprint.
+                    // in this case, just adjust the footprint and try again.
+                    if (tryAgain) {
+                        // keep the border constant
+                        int border = footprint - size;
+                        size = PS_MAX(stamps->normWindow1, stamps->normWindow2) + 2;
+                        footprint = size + border;
+
+                        // we need to reconstruct everything, so just free the stamps here and retry
+                        psFree(stamps);
+                    } else {
+                        // unrecoverable error
+                        psError(psErrorCodeLast(), false, "Unable to get stamp window.");
+                        goto MATCH_ERROR;
+                    }
+                }
+            }
+
+            // check on the kernel scaling -- if the kron-based radial moments are very different, adjust to match them
+            { 
+                // float fwhm1;
+                // float fwhm2;
+                // pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
+                // psAssert(isfinite(fwhm1), "fwhm 1 not set");
+                // psAssert(isfinite(fwhm2), "fwhm 2 not set");
+
+                // XXX this is BAD: depends on the relationship below:
+                // stamps->normWindow1 = 2.75*R1;
+                // stamps->normWindow2 = 2.75*R2;
+                float radMoment1 = stamps->normWindow1 / 2.75;
+                float radMoment2 = stamps->normWindow2 / 2.75;
+                pmSubtractionParamsScale(NULL, NULL, isisWidths, radMoment1, radMoment2);
+
+                // float maxFWHM = PS_MAX(fwhm1, fwhm2);
+                // float maxRadial = PS_MAX(radMoment1, radMoment2);
+                
+                // if (fabs(2.0*(maxFWHM - maxRadial)/(maxFWHM + maxRadial)) > 0.25) {
+                // if (1) {
+                // 
+                //     float scale = maxRadial / maxFWHM;
+                //     psLogMsg ("psModules.imcombine", PS_LOG_INFO, "Kron and FWHM scales are quite different, re-scale by %f to use Kron", scale);
+                //     
+                //     for (int i = 0; i < isisWidths->n; i++) {
+                //      isisWidths->data.F32[i] *= scale;
+                //     }
+                // }
+            }
+
+            // 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);
+            }
+
+            memCheck("kernels");
+
+            if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
+                pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej);
+                switch (newMode) {
+                  case PM_SUBTRACTION_MODE_1:
+                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2.");
+                    break;
+                  case PM_SUBTRACTION_MODE_2:
+                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1.");
+                    break;
+                  default:
+                    psError(psErrorCodeLast(), false, "Unable to determine subtraction order.");
+                    goto MATCH_ERROR;
+                }
+                subMode = newMode;
+            }
+
+            int numRejected = -1;       // Number of rejected stamps in each iteration
+            for (int k = 0; (k < iter) && (numRejected != 0); k++) {
+                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
+
+                bool tryAgain = true;
+                while (tryAgain) {
+                    // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
+                    // doesn't matter.
+                    if (!pmSubtractionStampsSelect(&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(&tryAgain, stamps, size)) {
+                        // if we failed, it might be due to the desired normWindow being larger than the current footprint.
+                        // in this case, just adjust the footprint and try again.
+                        if (tryAgain) {
+                            footprint = PS_MAX(stamps->normWindow1, stamps->normWindow2) + 2;
+
+                            // we need to reconstruct everything, so just free the stamps here and retry
+                            psFree(stamps);
+                        } else {
+                            // unrecoverable error
+                            psError(psErrorCodeLast(), false, "Unable to get stamp window.");
+                            goto MATCH_ERROR;
+                        }
+                    }
+                }
 
                 // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue
-                psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
-                if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-
-                // step 0b : calculate the moments, pass along to the next steps via stamps->normValue
-                // XXX this step is now skipped -- delete
-                psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
-                if (!pmSubtractionCalculateMoments(kernels, stamps)) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-
-                // step 1: generate the elements of the matrix equation Ax = B
-                psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-
-                // step 2: solve the matrix equation Ax = B
-                psTrace("psModules.imcombine", 3, "Solving kernel equations...\n");
-                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-                memCheck("  solve equation");
-
-                pmSubtractionVisualPlotConvKernels(kernels);
-
-                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
-                if (!deviations) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
-                    goto MATCH_ERROR;
-                }
-                memCheck("   calculate deviations");
-
-                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
-                numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
-                if (numRejected < 0) {
-                    psError(psErrorCodeLast(), false, "Unable to reject stamps.");
-                    psFree(deviations);
-                    goto MATCH_ERROR;
-                }
-                psFree(deviations);
-
-                memCheck("  reject stamps");
-            }
-
-            // if we hit the max number of iterations and we have rejected stamps, re-solve
-            if (numRejected > 0) {
-
-                // step 1: generate the elements of the matrix equation Ax = B
-                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;
-                }
-
-                // step 2: solve the matrix equation Ax = B
-                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;
-                }
-
-                pmSubtractionVisualPlotConvKernels(kernels);
-
-                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
-                if (!deviations) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
-                    goto MATCH_ERROR;
-                }
-                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
-                psFree(deviations);
-            }
-            psFree(stamps);
-            stamps = NULL;
-
-            memCheck("solution");
-
-            if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
-                psError(psErrorCodeLast(), false, "Unable to generate QA data");
-                goto MATCH_ERROR;
-            }
-            memCheck("diag outputs");
-
-            psTrace("psModules.imcombine", 2, "Convolving...\n");
-            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
-                                       kernelError, covarFrac, region, kernels, true, useFFT)) {
-                psError(psErrorCodeLast(), false, "Unable to convolve image.");
-                goto MATCH_ERROR;
-            }
-
-            psFree(kernels);
-            kernels = NULL;
-        }
+                psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
+                if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+
+                // on each iteration, we start from scratch
+                psFree(bestMatch);
+
+                // choose the spatial order and subtraction direction (1, 2, dual)
+                // XXX need to make these respect recipe somewhat
+                for (int order = 0; order <= N_TEST_ORDER; order++) {
+                    for (int j = 0; j < N_TEST_MODES; j++) {
+                        if (!pmSubtractionMatchAttempt(&bestMatch, kernels, stamps, TestModes[j], order, false)) {
+                            goto MATCH_ERROR;
+                        }
+                    }
+                }
+                
+                // reject the deviant stamps based on the stats of the best match
+                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
+                numRejected = pmSubtractionRejectStamps(kernels, stamps, bestMatch, subMask, rej);
+                if (numRejected < 0) {
+                    psError(psErrorCodeLast(), false, "Unable to reject stamps.");
+                    goto MATCH_ERROR;
+                }
+                memCheck("  reject stamps");
+            }
+
+            // apply the best fit so we are ready to roll
+            psLogMsg("psModules.imcombine", PS_LOG_INFO, "applying order: %d, mode: %d\n", bestMatch->spatialOrder, bestMatch->mode);
+            if (!pmSubtractionMatchAttempt(NULL, kernels, stamps, bestMatch->mode, bestMatch->spatialOrder, true)) {
+                goto MATCH_ERROR;
+            }
+            psFree(stamps);
+            psFree(bestMatch);
+            memCheck("solution");
+
+            if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
+                psError(psErrorCodeLast(), false, "Unable to generate QA data");
+                goto MATCH_ERROR;
+            }
+            memCheck("diag outputs");
+
+            psTrace("psModules.imcombine", 2, "Convolving...\n");
+            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
+                                       kernelError, covarFrac, region, kernels, true, useFFT)) {
+                psError(psErrorCodeLast(), false, "Unable to convolve image.");
+                goto MATCH_ERROR;
+            }
+
+            psFree(kernels);
+            kernels = NULL;
+        }
     }
     psFree(rng);
@@ -816 +858 @@
 
     if (conv1 && !pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) {
-        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
-        goto MATCH_ERROR;
+        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;
+        psError(psErrorCodeLast(), false, "Unable to set border of convolved image.");
+        goto MATCH_ERROR;
     }
 
@@ -832 +874 @@
 #ifdef TESTING
     {
-        if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
-            psFits *fits = psFitsOpen("convolved1.fits", "w");
-            psFitsWriteImage(fits, NULL, conv1->image, 0, NULL);
-            psFitsClose(fits);
-        }
-
-        if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
-            psFits *fits = psFitsOpen("convolved2.fits", "w");
-            psFitsWriteImage(fits, NULL, conv2->image, 0, NULL);
-            psFitsClose(fits);
-        }
+        if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
+            psFits *fits = psFitsOpen("convolved1.fits", "w");
+            psFitsWriteImage(fits, NULL, conv1->image, 0, NULL);
+            psFitsClose(fits);
+        }
+
+        if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
+            psFits *fits = psFitsOpen("convolved2.fits", "w");
+            psFitsWriteImage(fits, NULL, conv2->image, 0, NULL);
+            psFitsClose(fits);
+        }
     }
 #endif
@@ -858 +900 @@
     psFree(variance);
     psFree(rng);
+    psFree(bestMatch);
     return false;
 }
@@ -866 +909 @@
 // increment).
 static int subtractionOrderWidth(const psKernel *kernel, // Image
-                                 float bg, // Background in image
-                                 int size, // Maximum size
-                                 const psArray *models, // Buffer of models
-                                 const psVector *modelSums // Buffer of model sums
+                                 float bg, // Background in image
+                                 int size, // Maximum size
+                                 const psArray *models, // Buffer of models
+                                 const psVector *modelSums // Buffer of model sums
     )
 {
@@ -882 +925 @@
     psVector *chi2 = psVectorAlloc(size, PS_TYPE_F32); // chi^2 as a function of radius
     for (int sigma = 0; sigma < size; sigma++) {
-        double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian
-        psKernel *model = models->data[sigma]; // Model of interest
-        for (int y = yMin; y <= yMax; y++) {
-            for (int x = xMin; x <= xMax; x++) {
-                sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg);
-            }
-        }
-        float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian
-        double sumDev2 = 0.0;           // Sum of square deviations
-        for (int y = yMin; y <= yMax; y++) {
-            for (int x = xMin; x <= xMax; x++) {
-                float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation
-                sumDev2 += PS_SQR(dev);
-            }
-        }
-        chi2->data.F32[sigma] = sumDev2;
+        double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian
+        psKernel *model = models->data[sigma]; // Model of interest
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg);
+            }
+        }
+        float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian
+        double sumDev2 = 0.0;           // Sum of square deviations
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation
+                sumDev2 += PS_SQR(dev);
+            }
+        }
+        chi2->data.F32[sigma] = sumDev2;
     }
 
@@ -904 +947 @@
     float bestChi2 = INFINITY;          // Best chi^2
     for (int i = 0; i < size; i++) {
-        if (chi2->data.F32[i] < bestChi2) {
-            bestChi2 = chi2->data.F32[i];
-            bestIndex = i;
-        }
+        if (chi2->data.F32[i] < bestChi2) {
+            bestChi2 = chi2->data.F32[i];
+            bestIndex = i;
+        }
     }
     psFree(chi2);
@@ -916 +959 @@
 
 bool pmSubtractionOrderStamp(psVector *ratios, psVector *mask, const pmSubtractionStampList *stamps,
-                             const psArray *models, const psVector *modelSums,
-                             int index, float bg1, float bg2)
+                             const psArray *models, const psVector *modelSums,
+                             int index, float bg1, float bg2)
 {
     PS_ASSERT_VECTOR_NON_NULL(ratios, false);
@@ -931 +974 @@
     pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest
     psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE || stamp->status == PM_SUBTRACTION_STAMP_USED,
-             "We checked this earlier.");
+             "We checked this earlier.");
 
     // Widths of stars
@@ -938 +981 @@
 
     if (width1 == 0 || width2 == 0) {
-        ratios->data.F32[index] = NAN;
-        mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff;
+        ratios->data.F32[index] = NAN;
+        mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff;
     } else {
-        ratios->data.F32[index] = (float)width1 / (float)width2;
-        mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0;
-        psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n",
-                index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]);
+        ratios->data.F32[index] = (float)width1 / (float)width2;
+        mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0;
+        psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n",
+                index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]);
     }
 
@@ -978 +1021 @@
     psVector *modelSums = psVectorAlloc(size, PS_TYPE_F64); // Gaussian model sums
     for (int sigma = 0; sigma < size; sigma++) {
-        psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model
-        float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared
-        double sumGG = 0.0;         // Sum of square of Gaussian
-        for (int y = -size; y <= size; y++) {
-            int y2 = PS_SQR(y);     // y squared
-            for (int x = -size; x <= size; x++) {
-                float rad2 = PS_SQR(x) + y2; // Radius squared
-                float value = expf(-rad2 * invSigma2); // Model value
-                model->kernel[y][x] = value;
-                sumGG += PS_SQR(value);
-            }
-        }
-        models->data[sigma] = model;
-        modelSums->data.F64[sigma] = 1.0 / sumGG;
+        psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model
+        float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared
+        double sumGG = 0.0;         // Sum of square of Gaussian
+        for (int y = -size; y <= size; y++) {
+            int y2 = PS_SQR(y);     // y squared
+            for (int x = -size; x <= size; x++) {
+                float rad2 = PS_SQR(x) + y2; // Radius squared
+                float value = expf(-rad2 * invSigma2); // Model value
+                model->kernel[y][x] = value;
+                sumGG += PS_SQR(value);
+            }
+        }
+        models->data[sigma] = model;
+        modelSums->data.F64[sigma] = 1.0 / sumGG;
     }
 
     // Fit models to stamps
     for (int i = 0; i < stamps->num; i++) {
-        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
-            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
-            continue;
-        }
-
-        if (pmSubtractionThreaded()) {
-            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER");
-            psArrayAdd(job->args, 1, ratios);
-            psArrayAdd(job->args, 1, mask);
-            psArrayAdd(job->args, 1, stamps);
-            psArrayAdd(job->args, 1, models);
-            psArrayAdd(job->args, 1, modelSums);
-            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32);
-            PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32);
-            if (!psThreadJobAddPending(job)) {
-                return false;
-            }
-        } else {
-            if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {
-                psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i);
-                psFree(models);
-                psFree(modelSums);
-                psFree(ratios);
-                psFree(mask);
-                return false;
-            }
-        }
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
+            continue;
+        }
+
+        if (pmSubtractionThreaded()) {
+            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER");
+            psArrayAdd(job->args, 1, ratios);
+            psArrayAdd(job->args, 1, mask);
+            psArrayAdd(job->args, 1, stamps);
+            psArrayAdd(job->args, 1, models);
+            psArrayAdd(job->args, 1, modelSums);
+            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32);
+            if (!psThreadJobAddPending(job)) {
+                return false;
+            }
+        } else {
+            if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {
+                psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i);
+                psFree(models);
+                psFree(modelSums);
+                psFree(ratios);
+                psFree(mask);
+                return false;
+            }
+        }
     }
 
     if (!psThreadPoolWait(true)) {
-        psError(psErrorCodeLast(), false, "Error waiting for threads.");
-        psFree(models);
-        psFree(modelSums);
-        psFree(ratios);
-        psFree(mask);
-            return false;
+        psError(psErrorCodeLast(), false, "Error waiting for threads.");
+        psFree(models);
+        psFree(modelSums);
+        psFree(ratios);
+        psFree(mask);
+        return false;
     }
 
@@ -1041 +1084 @@
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     if (!psVectorStats(stats, ratios, NULL, mask, 0xff)) {
-        psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio.");
-        psFree(mask);
-        psFree(ratios);
-        psFree(stats);
-        return PM_SUBTRACTION_MODE_ERR;
+        psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio.");
+        psFree(mask);
+        psFree(ratios);
+        psFree(stats);
+        return PM_SUBTRACTION_MODE_ERR;
     }
     psFree(ratios);
@@ -1052 +1095 @@
     // XXX raise an error here or not?
     if (isnan(stats->robustMedian)) {
-        psFree(stats);
-        return PM_SUBTRACTION_MODE_ERR;
+        psFree(stats);
+        return PM_SUBTRACTION_MODE_ERR;
     }
 
@@ -1066 +1109 @@
 // 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
-                                )
+                                pmSubtractionKernels *kernels, // Kernel description
+                                const char *description, // Description for trace
+                                psImage *subMask,  // Subtraction mask
+                                float rej               // Rejection threshold
+    )
 {
     assert(stamps);
     assert(kernels);
 
+    psAbort("this function is not working");
+# if (0)
+    psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n");
+    if (!pmSubtractionConvolveStamps(stamps, kernels)) {
+        psError(psErrorCodeLast(), false, "Unable to convolve stamps.");
+        return false;
+    }
+
     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;
+    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
+        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 (!pmSubtractionSolveEquation(kernels, stamps)) {
+        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+        return false;
     }
 
@@ -1090 +1141 @@
     psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
     if (!deviations) {
-        psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
-        return false;
-    }
-
+        psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
+        return false;
+    }
+
+    // XXX this needs to be made consistent with the modified 'reject stamps' function
     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;
+        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);
-
-        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);
-    }
-
+        // Allow re-fit with reduced stamps set
+        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
+        if (!pmSubtractionCalculateEquation(stamps, 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)) {
+            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
+            return false;
+        }
+        psTrace("psModules.imcombine", 3, "Recalculate %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, "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);
+    }
+# endif
     return true;
 }
@@ -1138 +1190 @@
 
 pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels,
-                                        const psImage *subMask, float rej)
+                                        const psImage *subMask, float rej)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR);
@@ -1150 +1202 @@
 
     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;
+        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);
@@ -1165 +1217 @@
 
     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;
+        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);
@@ -1179 +1231 @@
     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);
+             "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;
+        bestStamps = stamps1;
+        bestKernels = kernels1;
     } else {
-        bestStamps = stamps2;
-        bestKernels = kernels2;
+        bestStamps = stamps2;
+        bestKernels = kernels2;
     }
 
@@ -1204 +1256 @@
 }
 
-
-bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
-                              float scaleRef, float scaleMin, float scaleMax)
+static float scaleRefOption = NAN;
+static float scaleMinOption = NAN;
+static float scaleMaxOption = NAN;
+static bool  scaleOption = false;
+
+bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax) { 
+
+    if (scale) {
+        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);
+    }
+
+    scaleRefOption = scaleRef;
+    scaleMinOption = scaleMin;
+    scaleMaxOption = scaleMax;
+    scaleOption = scale;
+    
+    return true;
+}
+
+bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2)
 {
-    PS_ASSERT_PTR_NON_NULL(kernelSize, false);
-    PS_ASSERT_PTR_NON_NULL(stampSize, false);
+    // 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(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 fwhm1;
-    float fwhm2;
-
-    pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
-    psAssert(isfinite(fwhm1), "fwhm 1 not set");
-    psAssert(isfinite(fwhm2), "fwhm 2 not set");
+
+    if (!scaleOption) return true;
+
+    // pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
+    // psAssert(isfinite(fwhm1), "fwhm 1 not set");
+    // psAssert(isfinite(fwhm2), "fwhm 2 not set");
     
     // 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;
+    float scale = PS_MAX(fwhm1, fwhm2) / scaleRefOption;      // Scaling factor
+
+    if (isfinite(scaleMinOption) && scale < scaleMinOption) {
+        scale = scaleMinOption;
+    }
+    if (isfinite(scaleMaxOption) && scale > scaleMaxOption) {
+        scale = scaleMaxOption;
     }
 
     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);
+        widths->data.F32[i] *= scale;
+    }
+    if (kernelSize) {
+        *kernelSize = *kernelSize * scale + 0.5;
+    }
+    if (stampSize) {
+        *stampSize = *stampSize * scale + 0.5;
+    }
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Scaling kernel parameters by %f", scale);
+    if (kernelSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified kernel size %d", *kernelSize);
+    if (stampSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified stamp size %d", *stampSize);
 
     return true;

trunk/psModules/src/imcombine/pmSubtractionMatch.h

@@ -8 +8 @@
 #include <pmSubtractionKernels.h>
 #include <pmSubtractionStamps.h>
+#include <pmSubtraction.h>
 
 /// Match two images
@@ -26 +27 @@
                         int size,       ///< Kernel half-size
                         int order,      ///< Spatial polynomial order
-                        const psVector *widths, ///< ISIS Gaussian widths
+                        psVector *widths, ///< ISIS Gaussian widths
                         const psVector *orders, ///< ISIS Polynomial orders
                         int inner,      ///< Inner radius for various kernel types
@@ -100 +101 @@
 
 /// Scale subtraction parameters according to the FWHMs of the inputs
-bool pmSubtractionParamsScale(
-    int *kernelSize,                    ///< Half-size of the kernel
-    int *stampSize,                     ///< Half-size of the stamp (footprint)
-    psVector *widths,                   ///< ISIS widths
-    float scaleRef,                     ///< Reference width for scaling
-    float scaleMin,                     ///< Minimum scaling ratio, or NAN
-    float scaleMax                      ///< Maximum scaling ratio, or NAN
+// bool pmSubtractionParamsScale(
+//     int *kernelSize,                    ///< Half-size of the kernel
+//     int *stampSize,                     ///< Half-size of the stamp (footprint)
+//     psVector *widths,                   ///< ISIS widths
+//     float scaleRef,                     ///< Reference width for scaling
+//     float scaleMin,                     ///< Minimum scaling ratio, or NAN
+//     float scaleMax                      ///< Maximum scaling ratio, or NAN
+//     );
+
+bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2);
+
+bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax);
+
+bool pmSubtractionMatchAttempt(
+    pmSubtractionQuality **bestMatch,
+    pmSubtractionKernels *kernels, 
+    pmSubtractionStampList *stamps, 
+    pmSubtractionMode mode, 
+    int spatialOrder, 
+    bool final
     );
 

trunk/psModules/src/imcombine/pmSubtractionParams.c

@@ -9 +9 @@
 
 #include "pmErrorCodes.h"
+#include "pmFPA.h"
+#include "pmSubtractionTypes.h"
+#include "pmSubtraction.h"
 #include "pmSubtractionStamps.h"
-#include "pmSubtraction.h"
 #include "pmSubtractionParams.h"
 
@@ -16 +18 @@
 
 #if 0
+// XXX this was moved to pmSubtraction.c in r15443 -- delete 
 // Convolve the reference stamp by the kernel
 static psKernel *convolveStamp(const pmSubtractionStamp *stamp, // Stamp to be convolved

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -27 +27 @@
 #include "pmSource.h"
 
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionVisual.h"
 
 #define STAMP_LIST_BUFFER 20            // Number of stamps to add to list at a time
@@ -122 +124 @@
             if ((image1 && image1->data.F32[y][x] < thresh1) ||
                 (image2 && image2->data.F32[y][x] < thresh2)) {
+                // fprintf (stderr, "%f,%f : thresh\n", xRaw, yRaw);
                 continue;
             }
@@ -366 +369 @@
 }
 
-
-pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image1,
-                                                const psImage *image2, const psImage *subMask,
-                                                const psRegion *region, float thresh1, float thresh2,
-                                                int size, int footprint, float spacing, float normFrac,
-                                                float sysErr, float skyErr, pmSubtractionMode mode)
+bool pmSubtractionStampsSelect(pmSubtractionStampList **stamps, // Stamps to read
+                               const pmReadout *ro1, // Readout 1
+                               const pmReadout *ro2, // Readout 2
+                               const psImage *subMask, // Mask for subtraction, or NULL
+                               psImage *variance,  // Variance map
+                               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
+                               pmSubtractionMode mode // Mode for subtraction
+    )
+{
+    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");
+
+    psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest
+
+    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
+                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
+    if (!*stamps) {
+        psError(psErrorCodeLast(), false, "Unable to find stamps.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
+        psError(psErrorCodeLast(), false, "Unable to extract stamps.");
+        return false;
+    }
+
+    pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1);
+    return true;
+}
+
+pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, 
+                                                const psImage *image1,
+                                                const psImage *image2, 
+                                                const psImage *subMask,
+                                                const psRegion *region, 
+                                                float thresh1, 
+                                                float thresh2,
+                                                int size, 
+                                                int footprint, 
+                                                float spacing, 
+                                                float normFrac,
+                                                float sysErr, 
+                                                float skyErr, 
+                                                pmSubtractionMode mode)
 {
     if (!image1 && !image2) {
@@ -429 +485 @@
         stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing,
                                              normFrac, sysErr, skyErr);
+    }
+
+    // XXX TEST : dump all stars in the stamps here
+    if (0) {
+        FILE *f = fopen ("stamp.dat", "w");
+        for (int i = 0; i < stamps->num; i++) {
+            psVector *xList = stamps->x->data[i];
+            psVector *yList = stamps->y->data[i]; // Coordinate lists
+            psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
+
+            for (int j = 0; j < xList->n; j++) {
+                fprintf (f, "%d %d  %f %f  %f\n", i, j, xList->data.F32[j], yList->data.F32[j], fluxList->data.F32[j]);
+            }
+        }
+        fclose (f);
     }
 
@@ -636 +707 @@
     }
 
+    int nTotal = 0;
+
     // Sort the list by flux, with the brightest last
     for (int i = 0; i < numStamps; i++) {
@@ -662 +735 @@
         stamps->y->data[i] = ySorted;
         stamps->flux->data[i] = fluxSorted;
-    }
-
+        nTotal += num;
+    }
+    // fprintf (stderr, "nTotal %d\n", nTotal);
+    
     return stamps;
 }
 
-
-bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize)
+// we are essentially using aperture photometry to determine the photometric match between the
+// images.  we need to choose an appropriate-sized aperture for this analysis.  If it is too
+// large, the measurement will be noisy (and possibly biased) due to the sky noise.  If it is
+// too small, or inconsistent, the measurement will be biased.  We use Kron-mag like aperture
+// scaled by the first radial moment.
+bool pmSubtractionStampsGetWindow(bool *tryAgain, pmSubtractionStampList *stamps, int kernelSize)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
 
+    // if we succeed, or fail with an unrecoverable error, do not try again
+    if (tryAgain) {
+        *tryAgain = false;
+    }
+
     int size = stamps->footprint; // Size of postage stamps
 
+    // window for moments calculations downstream
+    psFree (stamps->window);
+    stamps->window = psKernelAlloc(-size, size, -size, size);
+    psImageInit(stamps->window->image, 0.0);
+
+    // window1 and window2 are mean stamp images used here to measure the 
+    // first radial moment, and thus the normalization window
     psFree (stamps->window1);
     stamps->window1 = psKernelAlloc(-size, size, -size, size);
@@ -683 +774 @@
     psImageInit(stamps->window2->image, 0.0);
 
-    // Generate a weighting window based on the fwhms (20% larger than the largest)
-    { 
-        float fwhm1, fwhm2;
-
-        // XXX this is annoyingly hack-ish
-        pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
-
-        float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35;
-
-        psFree (stamps->window);
-        stamps->window = psKernelAlloc(-size, size, -size, size);
-        psImageInit(stamps->window->image, 0.0);
-
-        for (int y = -size; y <= size; y++) {
-            for (int x = -size; x <= size; x++) {
-                stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma));
-            }
-        }
+    // Generate an initial weighting window based on the fwhms (50% larger than the largest)
+    float fwhm1, fwhm2;
+
+    // XXX this is annoyingly hack-ish
+    pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
+    
+    float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35;
+    
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+            stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma));
+        }
     }
 
@@ -790 +875 @@
     psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2);
 
-# if (1)
-    // this block attempts to calculate the radius based on the first radial moment
+    // attempt to calculate the normalization window based on the first radial moment
     double Sr1 = 0.0;
     double Sr2 = 0.0;
@@ -809 +893 @@
     float R2 = Sr2 / Sf2;
 
-    stamps->normWindow1 = 2.0*R1;
-    stamps->normWindow2 = 2.0*R2;
-    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2);
-
-# else
-    // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent).
-    // It did not do very well (though a true curve-of-growth analysis might be better...)
-    bool done1 = false;
-    bool done2 = false;
-    double prior1 = 0.0;
-    double prior2 = 0.0;
-    double delta1o = 1.0;
-    double delta2o = 1.0;
-    for (int radius = 1; radius <= size && !(done1 && done2); radius++) {
-        double within1 = 0.0;
-        double within2 = 0.0;
-        for (int y = -radius; y <= radius; y++) {
-            for (int x = -radius; x <= radius; x++) {
-                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
-                    within1 += stamps->window1->kernel[y][x];
-                }
-                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
-                    within2 += stamps->window2->kernel[y][x];
-                }
-            }
-        }
-        double delta1 = (within1 - prior1) / within1;
-        if (!done1 && (fabs(delta1) < stamps->normFrac)) {
-            // interpolate to the radius at which delta2 is normFrac:
-            stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1);
-            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
-            done1 = true;
-        }
-        double delta2 = (within2 - prior2) / within2;
-        if (!done2 && (fabs(delta2) < stamps->normFrac)) {
-            // interpolate to the radius at which delta2 is normFrac:
-            stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2);
-            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
-            done2 = true;
-        }
-        psTrace("psModules.imcombine", 5, "Radius %d: %f (%f) and %f (%f)\n", radius, within1, delta1, within2, delta2);
-
-        prior1 = within1;
-        prior2 = within2;
-        delta1o = delta1;
-        delta2o = delta2;
-
-        // if (!done1 && (within1 > (1.0 - stamps->normFrac) * sum1)) {
-        //     stamps->normWindow1 = radius;
-        //     done1 = true;
-        // }
-        // if (!done2 && (within2 > (1.0 - stamps->normFrac) * sum2)) {
-        //     stamps->normWindow2 = radius;
-        //     done2 = true;
-        // }
-
-    }
-# endif
-
-    psTrace("psModules.imcombine", 3, "Normalisation window radii set to %f and %f\n", stamps->normWindow1, stamps->normWindow2);
-    if (stamps->normWindow1 == 0 || stamps->normWindow1 >= size) {
+    // Compare the Kron Radii (R1 & R2) to above to the FWHMs : if they are too discrepant, we will need to rescale
+    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii vs FWHMs 1: fwhm: %f, kron %f\n", fwhm1, R1);
+    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii vs FWHMs 2: fwhm: %f, kron %f\n", fwhm2, R2);
+
+    // XXX CAREFUL : in pmSubtractionMatch.c:703, we rely on this factor of 2.75..
+    stamps->normWindow1 = 2.75*R1;
+    stamps->normWindow2 = 2.75*R2;
+    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Windows from Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2);
+
+    // if the calculated normWindows are too large, we will fall off the stamps.  In this case, we need to try again.
+    if ((stamps->normWindow1 > size) || (stamps->normWindow2 > size)) { 
+        if (tryAgain) {
+            *tryAgain = true;
+        }
+        psFree (stats);
+        psFree (flux1);
+        psFree (flux2);
+        psFree (norm1);
+        psFree (norm2);
+        return false; 
+    }
+
+    // this is an unrecoverable error : something really bogus in the data
+    if (stamps->normWindow1 == 0) {
         psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1).");
+        psFree (stats);
+        psFree (flux1);
+        psFree (flux2);
+        psFree (norm1);
+        psFree (norm2);
         return false;
     }
-    if (stamps->normWindow2 == 0 || stamps->normWindow2 >= size) {
+    if (stamps->normWindow2 == 0) {
         psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2).");
+        psFree (stats);
+        psFree (flux1);
+        psFree (flux2);
+        psFree (norm1);
+        psFree (norm2);
         return false;
+    }
+
+    // Generate a weighting window based on the kron radii
+    float radius = 2.0 * PS_MAX(R1, R2);
+    psImageInit(stamps->window->image, 0.0);
+
+    // we use a top-hat window for the moments analysis
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+            if (hypot(x,y) > radius) continue;
+            stamps->window->kernel[y][x] = 1.0;
+        }
     }
 
@@ -890 +959 @@
         }
     }
-
-#if 0
-    {
-        psFits *fits = NULL;
-        fits = psFitsOpen ("window1.norm.fits", "w");
-        psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
-        psFitsClose (fits);
-        fits = psFitsOpen ("window2.norm.fits", "w");
-        psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
-        psFitsClose (fits);
-    }
-#endif
 
     psFree (stats);
@@ -1128 +1185 @@
             continue;
         }
+       
+        // XXX this is somewhat arbitrary...
+        if (source->errMag > 0.05) continue;
+        if (fabs(source->psfMag - source->apMag) > 0.5) continue;
+
         if (source->modelPSF) {
             x->data.F32[numOut] = source->modelPSF->params->data.F32[PM_PAR_XPOS];
@@ -1135 +1197 @@
             y->data.F32[numOut] = source->peak->yf;
         }
-        // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
         numOut++;
     }

trunk/psModules/src/imcombine/pmSubtractionStamps.h

@@ -1 +1 @@
 #ifndef PM_SUBTRACTION_STAMPS_H
 #define PM_SUBTRACTION_STAMPS_H
-
-#include <pslib.h>
-
-#include "pmSubtractionKernels.h"
-
-/// Status of stamp
-typedef enum {
-    PM_SUBTRACTION_STAMP_INIT,          ///< Initial state
-    PM_SUBTRACTION_STAMP_FOUND,         ///< Found a suitable source for this stamp
-    PM_SUBTRACTION_STAMP_CALCULATE,     ///< Calculate matrix and vector values for this stamp
-    PM_SUBTRACTION_STAMP_USED,          ///< Use this stamp
-    PM_SUBTRACTION_STAMP_REJECTED,      ///< This stamp has been rejected
-    PM_SUBTRACTION_STAMP_NONE           ///< No stamp in this region
-} pmSubtractionStampStatus;
-
-/// A list of stamps
-typedef struct {
-    long num;                           ///< Number of stamps
-    psArray *stamps;                    ///< The stamps
-    psArray *regions;                   ///< Regions for each stamp
-    psArray *x, *y;                     ///< Coordinates for possible stamps (or NULL)
-    psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
-    int footprint;                      ///< Half-size of stamps
-    float normFrac;                     ///< Fraction of flux in window for normalisation window
-    float normValue;                    ///< calculated normalization
-    float normValue2;                   ///< calculated normalization
-    psKernel *window1;                  ///< window function generated from ensemble of stamps (input 1)
-    psKernel *window2;                  ///< window function generated from ensemble of stamps (input 2)
-    psKernel *window;                   ///< weighting window function (sigma = 1.1 * MAX(fwhm))
-    float normWindow1;                  ///< Size of window for measuring normalisation
-    float normWindow2;                  ///< Size of window for measuring normalisation
-    float sysErr;                       ///< Systematic error
-    float skyErr;                       ///< increase effective readnoise
-} pmSubtractionStampList;
 
 /// Allocate a list of stamps
@@ -74 +40 @@
     );
 
-
-/// A stamp for image subtraction
-typedef struct {
-    float x, y;                         ///< Position
-    float flux;                         ///< Flux
-    float xNorm, yNorm;                 ///< Normalised position
-    psKernel *image1;                   ///< Reference image postage stamp
-    psKernel *image2;                   ///< Input image postage stamp
-    psKernel *weight;                   ///< Weight image (1/variance) postage stamp, or NULL
-    psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
-    psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
-    psImage *matrix;                    ///< Least-squares matrix, or NULL
-    psVector *vector;                   ///< Least-squares vector, or NULL
-    double norm;                        ///< Normalisation difference
-    double normI1;                       ///< Sum(flux) for image 1
-    double normI2;                       ///< Sum(flux) for image 2
-    double normSquare1;                 ///< Sum(flux^2) for image 1 (used for penalty)
-    double normSquare2;                 ///< Sum(flux^2) for image 2 (used for penalty)
-    pmSubtractionStampStatus status;    ///< Status of stamp
-    psVector *MxxI1;                    ///< second moments of convolution images
-    psVector *MyyI1;                    ///< second moments of convolution images
-    psVector *MxxI2;                    ///< second moments of convolution images
-    psVector *MyyI2;                    ///< second moments of convolution images
-    double MxxI1raw;
-    double MyyI1raw;
-    double MxxI2raw;
-    double MyyI2raw;
-} pmSubtractionStamp;
-
 /// Allocate a stamp
 pmSubtractionStamp *pmSubtractionStampAlloc(void);
+
+// find and extract the stamps
+bool pmSubtractionStampsSelect(pmSubtractionStampList **stamps, // Stamps to read
+                               const pmReadout *ro1, // Readout 1
+                               const pmReadout *ro2, // Readout 2
+                               const psImage *subMask, // Mask for subtraction, or NULL
+                               psImage *variance,  // Variance map
+                               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
+                               pmSubtractionMode mode // Mode for subtraction
+    );
+
 
 /// Find stamps on an image
@@ -172 +128 @@
 /// Calculate the window and normalisation window from the stamps
 bool pmSubtractionStampsGetWindow(
+    bool *tryAgain,                     ///< re-try with new stamp size?
     pmSubtractionStampList *stamps,     ///< List of stamps
     int kernelSize                      ///< Half-size of kernel

trunk/psModules/src/imcombine/pmSubtractionThreads.c

@@ -6 +6 @@
 #include <pslib.h>
 
+#include "pmSubtractionTypes.h"
 #include "pmSubtractionMatch.h"
 #include "pmSubtractionEquation.h"
@@ -33 +34 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 4);
+        psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 3);
         task->function = &pmSubtractionCalculateEquationThread;
+        psThreadTaskAdd(task);
+        psFree(task);
+    }
+
+    {
+        psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CONVOLVE_STAMP", 3);
+        task->function = &pmSubtractionConvolveStampThread;
         psThreadTaskAdd(task);
         psFree(task);

trunk/psModules/src/imcombine/pmSubtractionTypes.h

@@ -80 +80 @@
     int xMin, xMax, yMin, yMax;         ///< Bounds of image (for normalisation)
     long num;                           ///< Number of kernel components (not including the spatial ones)
+    psVector *fwhms;                    ///< requested fwhms of the kernel Gaussians (ISIS, HERM or DECONV_HERM)
+    psVector *orders;                   ///< polynomial orders for each Gaussian (ISIS, HERM or DECONV_HERM)
     psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM)
-    psVector *widths;                   ///< Gaussian FWHMs (ISIS, HERM or DECONV_HERM)
+    psVector *widths;                   ///< measured Gaussian FWHMs of Gauss*poly (ISIS, HERM or DECONV_HERM)
     psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
     psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM)
@@ -90 +92 @@
     int size;                           ///< The half-size of the kernel
     int inner;                          ///< The size of an inner region
+    int binning;                        ///< Binning used for the SPAM kernels
+    int ringsOrder;                     ///< 
     int spatialOrder;                   ///< The spatial order of the kernels
     int bgOrder;                        ///< The order for the background fitting

trunk/psModules/src/imcombine/pmSubtractionVisual.c

@@ -16 +16 @@
 
 #include "pmKapaPlots.h"
+#include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionEquation.h"
 #include "pmSubtractionKernels.h"
+#include "pmSubtractionVisual.h"
 
 #include "pmVisual.h"
@@ -210 +213 @@
         psImageOverlaySection(canvas, im, x0, y0, "=");
         stampNum++;
+
+        // renormalize the section
+        float maxValue = 0;
+        for (int iy = 0; iy < im->numRows; iy++) {
+            for (int ix = 0; ix < im->numCols; ix++) {
+                maxValue = PS_MAX(maxValue, im->data.F64[iy][ix]);
+            }
+        }
+        if (maxValue == 0.0) continue;
+        for (int iy = 0; iy < im->numRows; iy++) {
+            for (int ix = 0; ix < im->numCols; ix++) {
+                canvas->data.F64[y0 + iy][x0 + ix] /= maxValue;
+            }
+        }
     }
 
@@ -231 +248 @@
 }
 
+/** Plot the least-squares matrix of each stamp */
+bool pmSubtractionVisualPlotLeastSquaresResid (const pmSubtractionStampList *stamps, psImage *matrixIn, int nUsed) {
+
+    if (!pmVisualTestLevel("ppsub.chisq", 1)) return true;
+
+    if (!plotLeastSquares) return true;
+
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
+        return false;
+    }
+
+    psImage *matrixNorm = psImageCopy(NULL, matrixIn, PS_TYPE_F64);
+    {
+        // renormalize the matrix
+        float maxValue = 0;
+        for (int iy = 0; iy < matrixNorm->numRows; iy++) {
+            for (int ix = 0; ix < matrixNorm->numCols; ix++) {
+                maxValue = PS_MAX(maxValue, matrixNorm->data.F64[iy][ix]);
+            }
+        }
+        for (int iy = 0; iy < matrixNorm->numRows; iy++) {
+            for (int ix = 0; ix < matrixNorm->numCols; ix++) {
+                matrixNorm->data.F64[iy][ix] /= maxValue;
+            }
+        }
+    }
+
+    // Find the stamp size
+    int imageMax = -1;
+    int numStamps = 0;
+    psElemType type = PS_TYPE_F64;
+
+    for(int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i];
+        if (stamp == NULL) continue;
+
+        psImage *im = stamp->matrix;
+        if (im == NULL) continue;
+
+        imageMax = PS_MAX(imageMax, im->numCols);
+        imageMax = PS_MAX(imageMax, im->numRows);
+        numStamps++;
+        type = im->type.type;
+    }
+    if (imageMax == -1) return false;
+
+    int border = 15;
+    imageMax += border;
+    int tileRowCount = (int) ceil(sqrt(numStamps));
+    int canvasX = tileRowCount * (imageMax);
+    int canvasY = tileRowCount * (imageMax);
+    psImage *canvas = psImageAlloc (canvasX, canvasY, type);
+    psImageInit (canvas, NAN);
+
+    // overlay the images
+    int stampNum = 0;
+    int stampListNum = 0;
+    while (stampNum < numStamps) {
+        int x0 = (imageMax) * (stampNum % tileRowCount);
+        int y0 = (imageMax) * (stampNum / tileRowCount);
+
+        pmSubtractionStamp *stamp = stamps->stamps->data[stampListNum++];
+        if (stamp == NULL) continue;
+
+        psImage *im = stamp->matrix;
+        if (im == NULL) continue;
+
+        stampNum++;
+
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
+
+        // renormalize the section
+        float maxValue = 0;
+        for (int iy = 0; iy < im->numRows; iy++) {
+            for (int ix = 0; ix < im->numCols; ix++) {
+                maxValue = PS_MAX(maxValue, im->data.F64[iy][ix]);
+            }
+        }
+        if (maxValue == 0.0) continue;
+        for (int iy = 0; iy < im->numRows; iy++) {
+            for (int ix = 0; ix < im->numCols; ix++) {
+                canvas->data.F64[y0 + iy][x0 + ix] = (im->data.F64[iy][ix] / maxValue - matrixNorm->data.F64[iy][ix]) / matrixNorm->data.F64[iy][ix];
+            }
+        }
+    }
+
+    psImage *canvas32 = pmVisualImageToFloat(canvas);
+    pmVisualRangeImage(kapa2, canvas32, "Least_Squares", 0, -100.0, 100.0);
+
+    if (0) {
+        static int count = 0;
+        char filename[64];
+        sprintf (filename, "chisq.%02d.fits", count);
+        count ++;
+        psFits *fits = psFitsOpen (filename, "w");
+        psFitsWriteImage (fits, NULL, canvas32, 0, NULL);
+        psFitsClose (fits);
+    }
+
+    pmVisualAskUser(&plotLeastSquares);
+    psFree(canvas);
+    psFree(canvas32);
+    return true;
+}
+
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {
 
@@ -304 +426 @@
     }
 
-    // XXX clear the overlay(s) (red at least!)
+    // clear the overlay (red at least!)
+    KiiEraseOverlay (kapa2, "red");
 
     // get the kernel sizes
@@ -337 +460 @@
     int nKernels = 0;
 
+    // paste in the kernel images, scaled by sum2
     if (maxStamp->convolutions1) {
         // output image is a grid of NXsub by NYsub sub-images
@@ -359 +483 @@
             int yPix = ySub * (2*footprint + 1 + 3) + footprint;
             
-            double sum = 0.0;
             double sum2 = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
-                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
-                    sum += kernel->kernel[y][x];
                     sum2 += PS_SQR(kernel->kernel[y][x]);
                 }
             }
-            // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
+            float scale = sqrt(sum2) / PS_SQR(2*footprint + 1);
+            for (int y = -footprint; y <= footprint; y++) {
+                for (int x = -footprint; x <= footprint; x++) {
+                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x] / scale;
+                }
+            }
         }               
         pmVisualScaleImage(kapa2, output, "Image", 0, true);
@@ -401 +527 @@
             int yPix = ySub * (2*footprint + 1 + 3) + footprint;
             
-            double sum = 0.0;
             double sum2 = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
-                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
-                    sum += kernel->kernel[y][x];
                     sum2 += PS_SQR(kernel->kernel[y][x]);
                 }
             }
-            // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
+            float scale = sqrt(sum2) / PS_SQR(2*footprint + 1);
+            for (int y = -footprint; y <= footprint; y++) {
+                for (int x = -footprint; x <= footprint; x++) {
+                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x] / scale;
+                }
+            }
         }               
         pmVisualScaleImage(kapa2, output, "Image", 1, true);
@@ -468 +596 @@
     KiiLoadOverlay (kapa2, overlay, Noverlay, "red");
     FREE (overlay);
-    return true;
-}
-
-static int footprint = 0;
-static int NX = 0;
-static int NY = 0;
-static psImage *sourceImage      = NULL;
-static psImage *targetImage      = NULL;
-static psImage *residualImage    = NULL;
-static psImage *fresidualImage   = NULL;
-static psImage *differenceImage  = NULL;
-static psImage *convolutionImage = NULL;
-
-bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
-
-    if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
-
-    // generate 4 storage images large enough to hold the N stamps:
-
-    footprint = stamps->footprint;
-
-    float NXf = sqrt(stamps->num);
-    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
-    
-    float NYf = stamps->num / NX;
-    NY = (int) NYf == NY ? NYf : NYf + 1.0;
-
-    int NXpix = (2*footprint + 1) * NX;
-    NXpix += (NX > 1) ? 3 * NX : 0;
-
-    int NYpix = (2*footprint + 1) * NY;
-    NYpix += (NY > 1) ? 3 * NY : 0;
-
-    sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    
-    psImageInit (sourceImage,      0.0);
-    psImageInit (targetImage,      0.0);
-    psImageInit (residualImage,    0.0);
-    psImageInit (fresidualImage,   0.0);
-    psImageInit (differenceImage,  0.0);
-    psImageInit (convolutionImage, 0.0);
-
-    return true;
-}
-
-bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
-
-    if (!pmVisualTestLevel("ppsub.stamp", 1)) return true;
-
-    double sum;
-
-    int NXoff = index % NX;
-    int NYoff = index / NX;
-
-    int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
-    int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;
-
-    // insert the (target) kernel into the (target) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
-            sum += targetImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-
-    // insert the (source) kernel into the (source) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
-            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-
-    // insert the (convolution) kernel into the (convolution) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            convolutionImage->data.F32[y + NYpix][x + NXpix] = convolution->kernel[y][x];
-            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-    
-    // insert the (difference) kernel into the (difference) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
-            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-
-    // insert the (residual) kernel into the (residual) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
-            sum += residualImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-
-    // insert the (fresidual) kernel into the (fresidual) image:
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
-        }
-    }
     return true;
 }
@@ -611 +621 @@
 }
 
-bool pmSubtractionVisualShowFit(double norm) {
-
-    // for testing, dump the residual image and exit
-    if (0) {
-        psMetadata *header = psMetadataAlloc();
-        psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm);
-        psFits *fits = psFitsOpen("resid.fits", "w");
-        psFitsWriteImage(fits, header, residualImage, 0, NULL);
-        psFitsClose(fits);
-        // exit (0);
-    }
+static int footprint = 0;
+static int NX = 0;
+static int NY = 0;
+static psImage *sourceImage      = NULL;
+static psImage *targetImage      = NULL;
+static psImage *residualImage    = NULL;
+static psImage *fresidualImage   = NULL;
+static psImage *differenceImage  = NULL;
+static psImage *convolutionImage = NULL;
+
+bool pmSubtractionVisualShowFit(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) {
 
     if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
@@ -627 +637 @@
     if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
     if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
+
+    // set up holding images for the visualization
+    pmSubtractionVisualShowFitInit (stamps);
+
+    int numKernels = kernels->num;      // Number of kernels
+
+    psImage *polyValues = NULL;         // Polynomial values
+    psKernel *residual = psKernelAlloc(-stamps->footprint, stamps->footprint, -stamps->footprint, stamps->footprint); // Residual image
+
+    double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) { continue; }
+
+        // Calculate coefficients of the kernel basis functions
+        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
+        double background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Difference in background
+
+        psImageInit(residual->image, 0.0);
+
+        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
+            psKernel *target;           // Target postage stamp
+            psKernel *source;           // Source postage stamp
+            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
+            switch (kernels->mode) {
+              case PM_SUBTRACTION_MODE_1:
+                target = stamp->image2;
+                source = stamp->image1;
+                convolutions = stamp->convolutions1;
+                break;
+              case PM_SUBTRACTION_MODE_2:
+                target = stamp->image1;
+                source = stamp->image2;
+                convolutions = stamp->convolutions2;
+                break;
+              default:
+                psAbort("Unsupported subtraction mode: %x", kernels->mode);
+            }
+
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *convolution = convolutions->data[j]; // Convolution
+                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
+                    }
+                }
+            }
+            // visualize the target, source, convolution and residual
+            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
+        } else {
+            // Dual convolution
+            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
+            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
+            psKernel *image1 = stamp->image1; // The first image
+            psKernel *image2 = stamp->image2; // The second image
+
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
+                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
+                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
+                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
+
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
+                    }
+                }
+            }
+            // visualize the target, source, convolution and residual
+            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
+        }
+    }
+    pmSubtractionVisualShowFitImage(norm);
+
+    return true;
+}
+
+// generate 4 storage images large enough to hold the stamps:
+bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
+
+    footprint = stamps->footprint;
+
+    float NXf = sqrt(stamps->num);
+    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
+    
+    float NYf = stamps->num / NX;
+    NY = (int) NYf == NY ? NYf : NYf + 1.0;
+
+    int NXpix = (2*footprint + 1) * NX;
+    NXpix += (NX > 1) ? 3 * NX : 0;
+
+    int NYpix = (2*footprint + 1) * NY;
+    NYpix += (NY > 1) ? 3 * NY : 0;
+
+    sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    
+    psImageInit (sourceImage,      0.0);
+    psImageInit (targetImage,      0.0);
+    psImageInit (residualImage,    0.0);
+    psImageInit (fresidualImage,   0.0);
+    psImageInit (differenceImage,  0.0);
+    psImageInit (convolutionImage, 0.0);
+
+    return true;
+}
+
+bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
+
+    double sum;
+
+    int NXoff = index % NX;
+    int NYoff = index / NX;
+
+    int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
+    int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;
+
+    // insert the (target) kernel into the (target) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
+            sum += targetImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (source) kernel into the (source) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
+            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (convolution) kernel into the (convolution) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            convolutionImage->data.F32[y + NYpix][x + NXpix] = convolution->kernel[y][x];
+            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+    
+    // insert the (difference) kernel into the (difference) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
+            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (residual) kernel into the (residual) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
+            sum += residualImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (fresidual) kernel into the (fresidual) image:
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
+        }
+    }
+    return true;
+}
+
+bool pmSubtractionVisualShowFitImage(double norm) {
 
     KiiEraseOverlay (kapa1, "red");
@@ -673 +866 @@
     psVector *x = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
     psVector *y = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
+    psVector *dy = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
 
     graphdata.xmin = -1.0;
@@ -685 +879 @@
         x->data.F32[i] = i;
         y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
+        dy->data.F32[i] = kernels->solution1err->data.F64[i];
         graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]);
         graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);
     }
-    x->n = y->n = kernels->num;
+    x->n = y->n = dy->n = kernels->num;
 
     float range;
@@ -709 +904 @@
     graphdata.size = 0.5;
     graphdata.style = 2;
+    graphdata.etype |= 0x01;
 
     KapaPrepPlot   (kapa3, x->n, &graphdata);
     KapaPlotVector (kapa3, x->n, x->data.F32, "x");
     KapaPlotVector (kapa3, x->n, y->data.F32, "y");
+    KapaPlotVector (kapa3, x->n, dy->data.F32, "dym");
+    KapaPlotVector (kapa3, x->n, dy->data.F32, "dyp");
 
     psFree (x);
     psFree (y);
+    psFree (dy);
     psFree (polyValues);
+
+    pmVisualAskUser(NULL);
+    return true;
+}
+
+// plot log(flux) vs log(chisq), log(flux) vs log(moments), log(chisq) vs log(moments)
+bool pmSubtractionVisualPlotChisqAndMoments(psVector *fluxes, psVector *chisq, psVector *moments) {
+
+    Graphdata graphdata;
+    KapaSection section;
+
+    if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
+
+    if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false;
+
+    KapaClearSections (kapa3);
+    KapaInitGraph (&graphdata);
+    KiiResize(kapa3, 1500, 500);
+
+    psVector *lchi = psVectorAlloc (fluxes->n, PS_TYPE_F32);
+    psVector *lflx = psVectorAlloc (fluxes->n, PS_TYPE_F32);
+    psVector *lMxx = psVectorAlloc (fluxes->n, PS_TYPE_F32);
+
+    // construct the plot vectors
+    for (int i = 0; i < fluxes->n; i++) {
+        lchi->data.F32[i] = log10(chisq->data.F32[i]);
+        lflx->data.F32[i] = log10(fluxes->data.F32[i]);
+        lMxx->data.F32[i] = log10(moments->data.F32[i]);
+    }
+
+    section.bg = KapaColorByName ("none"); // XXX probably should be 'none'
+
+    // section 1: lflux vs lchi
+    section.dx = 0.33;
+    section.dy = 1.00;
+    section.x  = 0.00;
+    section.y  = 0.00;
+    section.name = psStringCopy ("flux.v.chi");
+    KapaSetSection (kapa3, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    pmVisualScaleGraphdata(&graphdata, lflx, lchi, false);
+    KapaSetLimits (kapa3, &graphdata);
+
+    KapaSetFont (kapa3, "helvetica", 14);
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "log(flux)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "log(chisq)", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    KapaPrepPlot   (kapa3, lflx->n, &graphdata);
+    KapaPlotVector (kapa3, lflx->n, lflx->data.F32, "x");
+    KapaPlotVector (kapa3, lflx->n, lchi->data.F32, "y");
+
+    // section 2: lflux vs lMxx
+    section.dx = 0.33;
+    section.dy = 1.00;
+    section.x  = 0.33;
+    section.y  = 0.00;
+    section.name = psStringCopy ("flux.v.mom");
+    KapaSetSection (kapa3, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    pmVisualScaleGraphdata(&graphdata, lflx, lMxx, false);
+    KapaSetLimits (kapa3, &graphdata);
+
+    KapaSetFont (kapa3, "helvetica", 14);
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "log(flux)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "log(moments)", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    KapaPrepPlot   (kapa3, lflx->n, &graphdata);
+    KapaPlotVector (kapa3, lflx->n, lflx->data.F32, "x");
+    KapaPlotVector (kapa3, lflx->n, lMxx->data.F32, "y");
+
+    // section 1: lflux vs lchi
+    section.dx = 0.33;
+    section.dy = 1.00;
+    section.x  = 0.66;
+    section.y  = 0.00;
+    section.name = psStringCopy ("chi.v.mom");
+    KapaSetSection (kapa3, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    pmVisualScaleGraphdata(&graphdata, lchi, lMxx, false);
+    KapaSetLimits (kapa3, &graphdata);
+
+    KapaSetFont (kapa3, "helvetica", 14);
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "log(chisq)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "log(moments)", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    KapaPrepPlot   (kapa3, lflx->n, &graphdata);
+    KapaPlotVector (kapa3, lflx->n, lchi->data.F32, "x");
+    KapaPlotVector (kapa3, lflx->n, lMxx->data.F32, "y");
 
     pmVisualAskUser(NULL);

trunk/psModules/src/imcombine/pmSubtractionVisual.h

@@ -6 +6 @@
 bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro);
 bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps);
+bool pmSubtractionVisualPlotLeastSquaresResid (const pmSubtractionStampList *stamps, psImage *matrixIn, int nUsed);
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub);
+
+bool pmSubtractionVisualShowFit(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels);
 bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps);
 bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index);
-bool pmSubtractionVisualShowFit(double norm);
+bool pmSubtractionVisualShowFitImage(double norm);
+
 bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels);
 bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels);
 bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps);
+bool pmSubtractionVisualPlotChisqAndMoments(psVector *fluxes, psVector *chisq, psVector *moments);
 
 #endif