Changeset 30631 for branches/czw_branch/20101203/psModules

branches/czw_branch/20101203

Property svn:mergeinfo changed

/branches/eam_branches/ipp-20101205 (added)	merged: 29959,30006-30013,30015,30019-30023,30027-30030,30076,30081-30098,30101-30110,30115-30117,30138-30142,30151-30177,30194-30197,30199-30202,30208-30213,30218-30226,30228-30233,30237,30252,30265-30271,30285-30289,30299-30305,30307-30312,30322-30323,30332-30335,30341-30344,30348-30351,30362,30383-30385,30399-30401,30421-30425,30440-30452,30474-30482,30500-30516,30523-30526,30530,30554,30576,30579-30581,30584-30585
/trunk	merged: 30588,30590,30592-30595,30598-30626,30629-30630

branches/czw_branch/20101203/psModules/src/extras/pmVisual.c

-              r29545
+              r30631
 #include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmFPAfile.h"
+#include "pmAstrometryObjects.h"
+#include "pmSubtractionStamps.h"
+#include "pmTrend2D.h"
+#include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmSpan.h"
+#include "pmFootprintSpans.h"
+#include "pmFootprint.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmModelFuncs.h"
+#include "pmModel.h"
+#include "pmSourceMasks.h"
+#include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
+#include "pmSource.h"
+#include "pmSourceFitModel.h"
+#include "pmPSF.h"
+#include "pmPSFtry.h"
+#include "pmFPAExtent.h"
+#include "pmAstrometryVisual.h"
+#include "pmSubtractionVisual.h"
+#include "pmStackVisual.h"
+#include "pmSourceVisual.h"
+bool pmSubtractionVisualClose(void);
+bool pmAstromVisualClose(void);
+bool pmSubtractionVisualClose(void);
+bool pmStackVisualClose(void);
+bool pmSourceVisualClose(void);
+// #include "pmHDU.h"
+// #include "pmFPA.h"
+// #include "pmFPAfile.h"
+// #include "pmAstrometryObjects.h"
+// #include "pmSubtractionStamps.h"
+// #include "pmTrend2D.h"
+// #include "pmResiduals.h"
+// #include "pmGrowthCurve.h"
+// #include "pmSpan.h"
+// #include "pmFootprintSpans.h"
+// #include "pmFootprint.h"
+// #include "pmPeaks.h"
+// #include "pmMoments.h"
+// #include "pmModelFuncs.h"
+// #include "pmModel.h"
+// #include "pmSourceMasks.h"
+// #include "pmSourceExtendedPars.h"
+// #include "pmSourceDiffStats.h"
+// #include "pmSource.h"
+// #include "pmSourceFitModel.h"
+// #include "pmPSF.h"
+// #include "pmPSFtry.h"
+// #include "pmFPAExtent.h"
+// #include "pmAstrometryVisual.h"
+// #include "pmSubtractionVisual.h"
+// #include "pmStackVisual.h"
+// #include "pmSourceVisual.h"
 # if (HAVE_KAPA)
 …
+}
+bool pmVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max) {
+    KiiImage image;
+    KapaImageData data;
+    Coords coords;
+    strcpy (coords.ctype, "RA---TAN");
+    image.data2d = inImage->data.F32;
+    image.Nx = inImage->numCols;
+    image.Ny = inImage->numRows;
+    strcpy (data.name, name);
+    strcpy (data.file, name);
+    data.zero = min;
+    data.range = max - min;
+    data.logflux = 0;
+    KiiSetChannel (kapaFD, channel);
+    KiiNewPicture2D (kapaFD, &image, &data, &coords);
+    return true;
+}
 psImage* pmVisualImageToFloat(psImage *image) {

branches/czw_branch/20101203/psModules/src/extras/pmVisual.h

r23242	r30631
64	64	const char *name, int channel, bool clip);
65	65
	66	bool pmVisualRangeImage (int kapaFD, psImage inImage, const char name, int channel, float min, float max);
66	67
67	68	/** Calculate statistics on an image.

branches/czw_branch/20101203/psModules/src/imcombine/Makefile.am

r26893	r30631
30	30	pmStackReject.h \
31	31	pmSubtraction.h \
	32	pmSubtractionTypes.h \
32	33	pmSubtractionAnalysis.h \
33	34	pmSubtractionEquation.h \

branches/czw_branch/20101203/psModules/src/imcombine/pmReadoutCombine.c

r24907	r30631
33	33	params->blank = 0;
34	34	params->nKeep = 0;
35		params->frac~~High~~ = 0.0;
	35	params->fracLow = 0.0;
36	36	params->fracHigh = 0.0;
37	37	params->iter = 1;

branches/czw_branch/20101203/psModules/src/imcombine/pmStackReject.c

r28405	r30631
7	7	#include <pslib.h>
8	8
	9	#include "pmFPA.h"
	10	#include "pmSubtractionTypes.h"
9	11	#include "pmSubtraction.h"
10	12	#include "pmSubtractionThreads.h"

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtraction.c

-              r29777
+              r30631
 #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"
 …
     if (convolutions) {
-        // Already done
         return convolutions;
+    }
 …
+}
+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)
 …
+    }
+#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;
 …
     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
 …
     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
 …
             // 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) {
 …
+    }
+    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);
+    }
 …
     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
 …
                 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);
+            }
 …
                 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);
+            }
 …
     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);
 …
       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);
 …
   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;
+}

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtraction.h

-              r29601
+              r30631
 #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
 …
 /// @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
 …
     );
+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
 …
 bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2);
+pmSubtractionQuality *pmSubtractionQualityAlloc();
 /// @}
 #endif

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionAnalysis.c

r29595	r30631
7	7
8	8	#include "pmFPA.h"
	9	#include "pmSubtractionTypes.h"
9	10	#include "pmSubtraction.h"
10	11	#include "pmSubtractionKernels.h"

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionDeconvolve.c

r26893	r30631
10	10
11	11	#include "pmFPA.h"
	12	#include "pmSubtractionTypes.h"
12	13	#include "pmSubtractionKernels.h"
13	14	#include "pmSubtractionDeconvolve.h"

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionEquation.c

-              r29598
+              r30631
 #include "pmErrorCodes.h"
+#include "pmVisual.h"
+#include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"
 …
 #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
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 …
                         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];
+                    }
 …
                     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;
 …
+}
+# define RENORM_BY_FLUX 0
 // Calculate the least-squares matrix and vector for dual convolution
 …
             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);
 …
             // 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);
 …
                         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;
 …
                         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;
 …
                         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];
+                    }
 …
             // 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
 …
                 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;
 …
                 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;
+                }
+            }
 …
         // 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
 …
     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;
 …
     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);
 …
     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);
 …
+}
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
+                                    const pmSubtractionEquationCalculationMode mode)
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
+{
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
 …
     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;
 …
             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);
+        }
+    }
 …
+    }
-    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;
 …
 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
 …
+        }
+        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);
 …
+        }
+        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
 #if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
 …
 #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);
 …
         // 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]);
+        }
 …
         psFree(sumVector);
         psFree(sumMatrix);
+        psFree(invMatrix);
     } else {
 …
+        }
+        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
 #if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
 …
+        }
+        // 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);
 …
         // 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:
 …
             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]);
+        }
 …
         kernels->solution1->data.F64[bgIndex] = 0.0;
+        psFree(solution);
+        psFree(sumVector);
         psFree(sumMatrix);
+        psFree(sumVector);
+        psFree(solution);
+        psFree(invMatrix);
+    }
 …
     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
 …
+}
+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
 …
     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);
 …
+            }
-            // 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++) {
 …
+            }
-            // 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++) {
 …
+        }
+        double flux = 0.0;
         double deviation = 0.0;         // Sum of differences
         for (int y = - footprint; y <= footprint; y++) {
 …
                 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);
 …
         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);
 …
     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;
+}
 …
 // 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) {

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionEquation.h

-              r29543
+              r30631
 #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
 …
 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
     );
 …
 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

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionEquation.v0.c

r29543	r30631
9	9	#include "pmErrorCodes.h"
10	10	#include "pmSubtraction.h"
	11	#include "pmSubtractionTypes.h"
11	12	#include "pmSubtractionKernels.h"
12	13	#include "pmSubtractionStamps.h"

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionHermitian.c

r26893	r30631
8	8	#include <pslib.h>
9	9
	10	#include "pmSubtractionTypes.h"
10	11	#include "pmSubtractionHermitian.h"
11	12

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionIO.c

r27321	r30631
12	12	#include "pmConceptsRead.h"
13	13
	14	#include "pmSubtractionTypes.h"
14	15	#include "pmSubtraction.h"
15	16	#include "pmSubtractionKernels.h"

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionKernels.c

-              r29597
+              r30631
 #include <pslib.h>
+#include "pmFPA.h"
+#include "pmSubtractionTypes.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"
 …
+{
     psFree(kernels->description);
+    psFree(kernels->fwhms);
+    psFree(kernels->orders);
     psFree(kernels->u);
     psFree(kernels->v);
 …
     psFree(kernels->solution1);
     psFree(kernels->solution2);
+    psFree(kernels->solution1err);
+    psFree(kernels->solution2err);
     psFree(kernels->sampleStamps);
+}
 …
     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)
 …
     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
 …
     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
 …
     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)
 …
 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)
+{
 …
     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);
 …
     kernels->size = size;
     kernels->inner = 0;
+    kernels->binning = 0;
+    kernels->ringsOrder = 0;
     kernels->spatialOrder = spatialOrder;
     kernels->bgOrder = 0;
 …
     kernels->solution1 = NULL;
     kernels->solution2 = NULL;
+    kernels->solution1err = NULL;
+    kernels->solution2err = NULL;
     kernels->mean = NAN;
     kernels->rms = NAN;
 …
     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)) {
 …
     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);
 …
     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);
 …
     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
 …
     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
 …
+}
+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)
+{
 …
     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);

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionKernels.h

-              r29601
+              r30631
 #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
 …
                                                 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
 …
     );
+bool pmSubtractionKernelsMakeDescription(pmSubtractionKernels *kernels);
 /// Copy kernels
 ///

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionMask.c

r27402	r30631
7	7
8	8	#include "pmErrorCodes.h"
	9	#include "pmFPA.h"
	10	#include "pmSubtractionTypes.h"
9	11	#include "pmSubtraction.h"
10	12	#include "pmSubtractionKernels.h"

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionMatch.c

-              r29596
+              r30631
 #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"
 …
 static bool useFFT = true;              // Do convolutions using FFT
-# define SUBMODE PM_SUBTRACTION_EQUATION_ALL
 //#define TESTING
 …
+}
-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
 …
                                   float badFrac,   // Maximum fraction of bad input pixels to accept
                                   pmSubtractionMode subMode // Mode of subtraction
+                                  )
+    )
+{
     if (subMode != PM_SUBTRACTION_MODE_2) {
 …
+}
+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,
 …
                         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,
 …
     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");
 …
             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) {
 …
+            }
+            // 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);
 …
     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;
+    }
 …
 #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
 …
     psFree(variance);
     psFree(rng);
+    psFree(bestMatch);
     return false;
+}
 …
 // 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
+    )
+{
 …
     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;
+    }
 …
     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);
 …
 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);
 …
     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
 …
     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]);
+    }
 …
     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;
+    }
 …
     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);
 …
     // 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;
+    }
 …
 // 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;
+    }
 …
     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;
+}
 …
 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);
 …
     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);
 …
     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);
 …
     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;
+    }
 …
+}
+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;

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionMatch.h

-              r29543
+              r30631
 #include <pmSubtractionKernels.h>
 #include <pmSubtractionStamps.h>
+#include <pmSubtraction.h>
 /// Match two images
 …
                         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
 …
 /// 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
     );

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionParams.c

-              r27086
+              r30631
 #include "pmErrorCodes.h"
+#include "pmFPA.h"
+#include "pmSubtractionTypes.h"
+#include "pmSubtraction.h"
 #include "pmSubtractionStamps.h"
-#include "pmSubtraction.h"
 #include "pmSubtractionParams.h"
 …
 #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

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionStamps.c

-              r29543
+              r30631
 #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
 …
             if ((image1 && image1->data.F32[y][x] < thresh1) ||
                 (image2 && image2->data.F32[y][x] < thresh2)) {
+                // fprintf (stderr, "%f,%f : thresh\n", xRaw, yRaw);
                 continue;
+            }
 …
+}
+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) {
 …
         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);
+    }
 …
+    }
+    int nTotal = 0;
     // Sort the list by flux, with the brightest last
     for (int i = 0; i < numStamps; i++) {
 …
         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);
 …
     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));
+        }
+    }
 …
     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;
 …
     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;
+        }
+    }
 …
+        }
+    }
-#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);
 …
             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];
 …
             y->data.F32[numOut] = source->peak->yf;
+        }
-        // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
         numOut++;
+    }

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionStamps.h

-              r29543
+              r30631
 #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
 …
     );
-/// 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
 …
 /// 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

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionThreads.c

-              r27365
+              r30631
 #include <pslib.h>
+#include "pmSubtractionTypes.h"
 #include "pmSubtractionMatch.h"
 #include "pmSubtractionEquation.h"
 …
+    {
         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);

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionVisual.c

-              r29600
+              r30631
 #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"
 …
         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;
+            }
+        }
+    }
 …
+}
+/** 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) {
 …
+    }
+    // XXX clear the overlay(s) (red at least!)
+    // clear the overlay (red at least!)
+    KiiEraseOverlay (kapa2, "red");
     // get the kernel sizes
 …
     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
 …
             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);
 …
             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);
 …
     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;
+}
 …
+}
+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;
 …
     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");
 …
     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;
 …
         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;
 …
     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);

branches/czw_branch/20101203/psModules/src/imcombine/pmSubtractionVisual.h

-              r29601
+              r30631
 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

branches/czw_branch/20101203/psModules/src/objects

Property svn:mergeinfo deleted

branches/czw_branch/20101203/psModules/src/objects/pmFootprint.c

-              r29004
+              r30631
+}
+// create a new footprint with the same span set as the input footprint
+pmFootprint *pmFootprintCopyData(pmFootprint *inFoot, psImage *image) {
+    pmFootprint *outFoot = pmFootprintAlloc(inFoot->nspans, image);
+    for (int i = 0; i < inFoot->nspans; i++) {
+        pmSpan *span = inFoot->spans->data[i];
+        pmFootprintAddSpan(outFoot, span->y, span->x0, span->x1);
+    }
+    return outFoot;
+}
 /************************************************************************************************************/

branches/czw_branch/20101203/psModules/src/objects/pmFootprint.h

r29004	r30631
27	27
28	28	bool pmFootprintAllocEmptySpans (pmFootprint *footprint, int nSpans);
	29
	30	pmFootprint pmFootprintCopyData(pmFootprint inFoot, psImage *image);
29	31
30	32	pmFootprint pmFootprintNormalize(pmFootprint fp);

branches/czw_branch/20101203/psModules/src/objects/pmFootprintAssignPeaks.c

-              r29004
+              r30631
 /*
  * Given a psArray of pmFootprints and another of pmPeaks, assign the peaks to the
  * footprints in which that fall; if they _don't_ fall in a footprint, add a suitable
+ * footprints in which they fall; if they _don't_ fall in a footprint, add a suitable
  * one to the list.
  */
 …
         // XXX are we allowed to have peak-less footprints??
+        if (!fp->peaks->n) continue;
+        if (fp->peaks->n == 0) continue;
+        if (fp->peaks->n == 1) continue;
         fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN);
+        // XXX check for an assert on duplicates (I don't think they can happen, but
+        // let's double check for now)
+        for (int j = 1; j < fp->peaks->n; j++) {
+            psAssert (fp->peaks->data[j] != fp->peaks->data[j-1], "duplicate peak!");
+        }
+        continue;
+        // XXX WHY am I culling duplicates?  how can there be duplicates?
         // XXX EAM : the algorithm below should be much faster than using psArrayRemove if
         // the number of peaks in the footprint is large, or if there are no duplicates.
 …
         // track the number of good peaks in the footprint
         int nGood = 1;
+        int lastGood = 0;
         // check for duplicates
 …
         // (if sorted list has A, B, A, B ...)
         for (int j = 1; j < fp->peaks->n; j++) {
             if (fp->peaks->data[j] == fp->peaks->data[nGood]) {
+            if (fp->peaks->data[j] == fp->peaks->data[lastGood]) {
                 // everything on the array has its own mem reference; free and drop this one
                 psFree (fp->peaks->data[j]);
                 fp->peaks->data[j] = NULL;
             } else {
                 nGood ++;
+                lastGood ++;
+            }
+        }
+        int nGood = lastGood + 1;
         // no deleted peaks, go to next footprint

branches/czw_branch/20101203/psModules/src/objects/pmFootprintCullPeaks.c

-              r29004
+              r30631
   * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
   * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
   * to reach a still higher peak --- and if that coll's more than nsigma DN below your
+  * to reach a still higher peak --- and if that coll's more (less?) than nsigma DN below your
   * starting point, discard the peak.
   */
 …
         float threshold = flux - nsigma_delta*stdev_pad;
         if (isnan(threshold) || threshold < min_threshold) {
+        if (isnan(threshold)) {
             // min_threshold is assumed to be below the detection threshold,
             // so all the peaks are pmFootprint, and this isn't the brightest
 …
             threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
+        }
+        if (threshold < min_threshold) {
+            threshold = min_threshold;
+        }
         // init peakFootprint here?

branches/czw_branch/20101203/psModules/src/objects/pmModel.c

-              r29004
+              r30631
     tmp->fitRadius = 0;
     tmp->flags = PM_MODEL_STATUS_NONE;
+    tmp->residuals = NULL;              // XXX should the model own this memory?
+    tmp->residuals = NULL;              // do not free: the model does not own this memory
+    tmp->isPCM = false;
     psS32 Nparams = pmModelClassParameterCount(type);
 …
 pmModel *pmModelCopy (pmModel *model)
+{
+    PS_ASSERT_PTR_NON_NULL(model, NULL);
+    if (model == NULL) {
+        return NULL;
+    }
     pmModel *new = pmModelAlloc (model->type);

branches/czw_branch/20101203/psModules/src/objects/pmModel.h

r29004	r30631
43	43	float fitRadius; ///< fit radius actually used
44	44	pmResiduals *residuals; ///< normalized PSF residuals
	45	bool isPCM; ///< is this model fitted with PSF-convolution?
45	46
46	47	// functions for this model which depend on the model class

branches/czw_branch/20101203/psModules/src/objects/pmPCMdata.h

r29004	r30631
83	83	pmPCMdata *pcm);
84	84
	85	psKernel pmPCMkernelFromPSF (pmSource source, int nPix);
	86
85	87	bool pmSourceModelGuessPCM (pmPCMdata pcm, pmSource source, psImageMaskType maskVal, psImageMaskType markVal);
86	88

branches/czw_branch/20101203/psModules/src/objects/pmPSFtry.c

-              r29004
+              r30631
     for (int i = 0; i < sources->n; i++) {
+        test->sources->data[i] = pmSourceCopy (sources->data[i]);
+        pmSource *sourceOld = sources->data[i];
+        pmSource *sourceNew = pmSourceCopy (sourceOld);
+        // save a reference so we can get back to the original
+        // this is specifically used in psphotChooosePSF to unflag the candidate PSF sources
+        // which were not actually used to generate a PSF model
+        sourceNew->parent = sourceOld;
+        test->sources->data[i] = sourceNew;
+    }

branches/czw_branch/20101203/psModules/src/objects/pmPSFtryFitEXT.c

-              r29004
+              r30631
         if (!source->moments) {
             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y);
             continue;
+        }
         if (!source->moments->nPixels) {
+            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y);
             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
             continue;
 …
         source->modelEXT = pmSourceModelGuess (source, options->type);
         if (source->modelEXT == NULL) {
+            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
             continue;
+        }
 …
         // exclude the poor fits
         if (!status) {
+            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
             continue;
+        }

branches/czw_branch/20101203/psModules/src/objects/pmPSFtryFitPSF.c

r29004	r30631
81	81	if (source->modelPSF == NULL) {
82	82	psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
	83	psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
83	84	continue;
84	85	}

branches/czw_branch/20101203/psModules/src/objects/pmPSFtryModel.c

-              r30118
+              r30631
 // mask values indicate the reason the source was rejected:
+// XXX some test code (delete eventually)
+// bool countMaskedSources(pmPSFtry *try, char *msg) {
+//
+//     int N1 = 0;
+//     int N2 = 0;
+//     for (int i = 0; i < try->sources->n; i++) {
+//         pmSource *source = try->sources->data[i];
+//      // fprintf (stderr, "%llx : %d\n", (long long int) source, (source->mode & PM_SOURCE_MODE_PSFSTAR));
+//      if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
+//          N1 ++;
+//         }
+//      if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+//          N2 ++;
+//         }
+//     }
+//     fprintf (stderr, "%s : masked: %d or %d of %ld\n", msg, N1, N2, try->sources->n);
+//     return true;
+// }
 // generate a pmPSFtry with a copy of the test PSF sources
 pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal)
 …
             chisq->data.F32[i] = 0.0;
             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
+            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : no PSF model\n", i, source->peak->x, source->peak->y);
         } else {
             flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0];

branches/czw_branch/20101203/psModules/src/objects/pmSource.c

-              r30118
+              r30631
     psFree(tmp->moments);
     psFree(tmp->diffStats);
     psFree(tmp->radial);
+    psFree(tmp->radialAper);
     psFree(tmp->blends);
     psTrace("psModules.objects", 10, "---- end ----\n");
 …
     source->extpars = NULL;
     source->diffStats = NULL;
+    source->radial = NULL;
+    source->radialAper = NULL;
+    source->parent = NULL;
     source->region = psRegionSet(NAN, NAN, NAN, NAN);
 …
 affecting the original.  This Copy can be used to allow multiple fit attempts on the same
 object.  The pixels, variance, and mask arrays all point to the same original subarrays.  The
 peak and moments point at the original values.
+peak and moments point at the original values.  The models, blends, and XXX are NOT copied
 *****************************************************************************/
 pmSource *pmSourceCopy(pmSource *in)
 …
+    }
     pmSource *source = pmSourceAlloc ();
-    // keep the original ID so we can find map back to the original
-    P_PM_SOURCE_SET_ID(source, in->id);
     // peak has the same values as the original
 …
+}
-/******************************************************************************
-pmSourceCopyData(): this creates a new, duplicate source with the same parameters as the
-original (but is actually a new source at the same location)
-*****************************************************************************/
-pmSource *pmSourceCopyData(pmSource *in)
+{
-    if (in == NULL) {
-        return(NULL);
+    }
-    // this copy is used to allow multiple fit attempts on the
-    // same object.  the pixels, variance, and mask arrays all point to
-    // the same original subarrays.  the peak and moments point at
-    // the original values.
-    pmSource *source = pmSourceAlloc ();
-    // this is actually the same peak as the original, is this the correct way to handle this?
-    if (in->peak != NULL) {
-        source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->value, in->peak->type);
-        source->peak->xf = in->peak->xf;
-        source->peak->yf = in->peak->yf;
-        source->peak->flux = in->peak->flux;
-        source->peak->SN = in->peak->SN;
+    }
-    // copy the values in the moments structure
-    if (in->moments != NULL) {
-        source->moments  =  pmMomentsAlloc();
-        *source->moments = *in->moments;
+    }
-    // These images are all views to the parent.
-    // We want a new view, but pointing at the same pixels.
-    source->pixels   = psImageCopyView(NULL, in->pixels);
-    source->variance   = psImageCopyView(NULL, in->variance);
-    source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;
-    // the maskObj is a unique mask array; create a new mask image
-    source->maskObj = in->maskObj ? psImageCopy (NULL, in->maskObj, PS_TYPE_IMAGE_MASK) : NULL;
-    source->type = in->type;
-    source->mode = in->mode;
-    source->imageID = in->imageID;
-    return(source);
+}
 // x,y are defined in the parent image coords of readout->image
 bool pmSourceDefinePixels(pmSource *mySource,
 …
     newRegion = psRegionForImage (readout->image, newRegion);
+    // re-define if required by region or absence of pixels
     extend = false;
     extend |= (int)(newRegion.x0) < (int)(mySource->region.x0);
 …
+    }
     return extend;
+}
+bool pmSourceRedefinePixelsByRegion(pmSource *mySource,
+                                    const pmReadout *readout,
+                                    psRegion newRegion)
+{
+    PS_ASSERT_PTR_NON_NULL(mySource, false);
+    PS_ASSERT_PTR_NON_NULL(readout, false);
+    PS_ASSERT_PTR_NON_NULL(readout->image, false);
+    // re-create the subimage
+    psFree (mySource->pixels);
+    psFree (mySource->variance);
+    psFree (mySource->maskView);
+    mySource->pixels   = psImageSubset(readout->image,    newRegion);
+    mySource->variance = psImageSubset(readout->variance, newRegion);
+    mySource->maskView = psImageSubset(readout->mask,     newRegion);
+    mySource->region   = newRegion;
+    // re-copy the main mask pixels.  NOTE: the user will need to reset the object mask
+    // pixels (eg, with psImageKeepCircle)
+    mySource->maskObj = psImageCopy (mySource->maskObj, mySource->maskView, PS_TYPE_IMAGE_MASK);
+    // drop the old modelFlux pixels and force the user to re-create
+    psFree (mySource->modelFlux);
+    mySource->modelFlux = NULL;
+    // drop the old psfImage pixels and force the user to re-create
+    psFree (mySource->psfImage);
+    mySource->psfImage = NULL;
+    return true;
+}
 …
             if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
                 source->type = PM_SOURCE_TYPE_STAR;
                 source->mode |= PM_SOURCE_MODE_PSFSTAR;
+                source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR;
                 Npsf ++;
                 continue;
 …
+}
+// construct a realization of the source model
+// XXX this function should optionally save an existing psf image from modelFlux
+// construct a realization of the PSF model at the location of this source
 bool pmSourceCachePSF (pmSource *source, psImageMaskType maskVal) {
     PS_ASSERT_PTR_NON_NULL(source, false);
 …
+    }
+    if (!addNoise) {
+        if (add) {
+            status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+            source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+        } else {
+            status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+            source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+        }
+    if (add) {
+        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+        source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+    } else {
+        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+    }

branches/czw_branch/20101203/psModules/src/objects/pmSource.h

-              r29546
+              r30631
     PM_SOURCE_TMPF_SIZE_CR_CANDIDATE = 0x0008,
     PM_SOURCE_TMPF_MOMENTS_MEASURED  = 0x0010,
+    PM_SOURCE_TMPF_CANDIDATE_PSFSTAR = 0x0020,
 } pmSourceTmpF;
 …
  *  simplest measurement of a source is the location and flux of the peak pixel
  *  associated with the source:
+ *
+ * a pmSource is the information about a (possible) blob of flux in a specific image.  A source
+ * may represent an insignificant or undetected source.  There may be multiple representations
+ * of an image (eg, alternate smoothed copies); sources on alternate images may have a pointer
+ * to the version on the primary image (source->parent).  A set of sources on different, but
+ * related images (eg, multiple exposures or different filters) which (may) represent the same
+ * astronomical object are grouped together with the pmPhotObj type (set pmPhotObj.h).
+ *
+ * A single source may be fitted with multiple models (not at the same time!).  The PSF model
+ * fit is a fit of the position (optionally) and the flux to the PSF model at the location of
+ * the source.  Alternate model fits are extended source models. The best model fit is used to
+ * subtract the object from the image.
+ *
  *  XXX do I have to re-organize this (again!) to allow an arbitrary set of extended model fits??
 …
     pmPeak  *peak;                      ///< Description of peak pixel.
     psImage *pixels;                    ///< Rectangular region including object pixels.
     psImage *variance;                    ///< Image variance.
+    psImage *variance;                  ///< Image variance.
     psImage *maskObj;                   ///< unique mask for this object which marks included pixels associated with objects.
     psImage *maskView;                  ///< view into global image mask for this object region
     psImage *modelFlux;                 ///< cached copy of the best model for this source
     psImage *psfImage;                   ///< cached copy of the psf model for this source
+    psImage *psfImage;                  ///< cached copy of the psf model for this source
     pmMoments *moments;                 ///< Basic moments measured for the object.
     pmModel *modelPSF;                  ///< PSF Model fit (parameters and type)
 …
     pmSourceExtendedPars *extpars;      ///< extended source parameters
     pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
+    pmSourceRadialApertures *radial;    ///< radial flux in circular apertures
+    psArray *radialAper;                ///< radial flux in circular apertures
+    pmSource *parent;                   ///< reference to the master source from which this is derived
     int imageID;
 };
 …
 );
+bool pmSourceRedefinePixelsByRegion (
+    pmSource *mySource,   ///< source to be re-defined
+    const pmReadout *readout,   ///< base the source on this readout
+    psRegion newRegion ///< region for source pixel definition
+);
 /** pmSourcePSFClump()
+ *

branches/czw_branch/20101203/psModules/src/objects/pmSourceExtendedPars.c

-              r28013
+              r30631
 #endif
-// #include <stdio.h>
-// #include <math.h>
-// #include <string.h>
 #include <pslib.h>
-// #include "pmHDU.h"
-// #include "pmFPA.h"
-// #include "pmFPAMaskWeight.h"
-// #include "pmSpan.h"
-// #include "pmFootprint.h"
-// #include "pmPeaks.h"
-// #include "pmMoments.h"
-// #include "pmResiduals.h"
-// #include "pmGrowthCurve.h"
-// #include "pmTrend2D.h"
-// #include "pmPSF.h"
-// #include "pmModel.h"
-// #include "pmSource.h"
 #include "pmSourceExtendedPars.h"

branches/czw_branch/20101203/psModules/src/objects/pmSourceFitPCM.c

-              r29004
+              r30631
     if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;
+    // once we have fitted a model, we need to record that this model is a PCM model:
+    pcm->modelConv->isPCM = true;
     // models can go insane: reject these
     bool onPic = true;

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO.c

-              r29546
+              r30631
 // lookup the EXTNAME values used for table data and image header segments
+static bool sourceExtensions(psString *headname, // Extension name for header
+                             psString *dataname, // Extension name for data
+                             psString *deteffname, // Extension name for detection efficiency
+                             psString *xsrcname, // Extension name for extended sources
+                             psString *xfitname, // Extension name for extended fits
+                             const pmFPAfile *file, // File of interest
+                             const pmFPAview *view // View to level of interest
+                             )
+bool pmSourceIOextnames(psString *headname,    // Extension name for image header
+                        psString *dataname,    // Extension name for PSF table data
+                        psString *deteffname,  // Extension name for detection efficiency
+                        psString *xsrcname,    // Extension name for extended non-parametric measurements
+                        psString *xfitname,    // Extension name for extended fitted measurements
+                        psString *xradname,    // Extension name for radial apertures
+                        const pmFPAfile *file, // File of interest
+                        const pmFPAview *view  // View to level of interest
+    )
+{
     bool status;                        // Status of MD lookup
 …
+    }
     // EXTNAME for table data
+    // EXTNAME for PSF table data
     if (dataname) {
         const char *rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
 …
+    }
     // EXTNAME for extended source data table
+    // EXTNAME for extended source non-parametric measurements
     if (xsrcname) {
         const char *rule = psMetadataLookupStr(&status, menu, "CMF.XSRC");
 …
+    }
+    // EXTNAME for extended source fitted measurements
     if (xfitname) {
-        // EXTNAME for extended source data table
         const char *rule = psMetadataLookupStr(&status, menu, "CMF.XFIT");
         if (!rule) {
 …
+        }
         *xfitname = pmFPAfileNameFromRule (rule, file, view);
+    }
+    // EXTNAME for radial apertures
+    if (xradname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XRAD");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XRAD in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *xradname = pmFPAfileNameFromRule (rule, file, view);
+    }
 …
             status &= pmSourcesWrite_##TYPE##_XFIT (file->fits, readout, sources, file->header, xfitname); \
         }                                                               \
+        if (xradname) {                                                 \
+            status &= pmSourcesWrite_##TYPE##_XRAD (file->fits, readout, sources, file->header, xradname, recipe); \
+        }                                                               \
+    }
 …
         bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS");
         bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS");
+        bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES");
         // define the EXTNAME values for the different data segments:
 …
         psString xsrcname = NULL;
         psString xfitname = NULL;
+        if (!sourceExtensions(&headname, &dataname, &deteffname,
+                              XSRC_OUTPUT ? &xsrcname : NULL,
+                              XFIT_OUTPUT ? &xfitname : NULL,
+                              file, view)) {
+        psString xradname = NULL;
+        if (!pmSourceIOextnames(&headname, &dataname, &deteffname,
+                                XSRC_OUTPUT ? &xsrcname : NULL,
+                                XFIT_OUTPUT ? &xfitname : NULL,
+                                XRAD_OUTPUT ? &xradname : NULL,
+                                file, view)) {
             return false;
+        }
 …
                 psMetadataAddStr (outhead, PS_LIST_TAIL, "XFITNAME", PS_META_REPLACE, "name of XFIT table extension", xfitname);
+            }
+            if (xradname) {
+                psMetadataAddStr (outhead, PS_LIST_TAIL, "XRADNAME", PS_META_REPLACE, "name of XRAD table extension", xradname);
+            }
             // these are case-sensitive since the EXTYPE is case-sensitive
 …
+        }
         // write out the detection efficiency TABLE segments
         if (deteffname) {
 …
         psFree (xsrcname);
         psFree (xfitname);
+        psFree (xradname);
         psFree (deteffname);
 …
         psFree (xsrcname);
         psFree (xfitname);
+        psFree (xradname);
         psFree (deteffname);
         return false;
 …
         psString dataname = NULL;
         psString deteffname = NULL;
         if (!sourceExtensions(&headname, &dataname, &deteffname, NULL, NULL, file, view)) {
+        if (!pmSourceIOextnames(&headname, &dataname, &deteffname, NULL, NULL, NULL, file, view)) {
             return false;
+        }

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO.h

-              r29546
+              r30631
 bool pmSourcesWrite_SMPDATA_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_SMPDATA_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_SMPDATA_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_DEV_0(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_DEV_0_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_DEV_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_PS1_DEV_0_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_DEV_1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_DEV_1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_DEV_1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_PS1_DEV_1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_CAL_0(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_CAL_0_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_CAL_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_PS1_CAL_0_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_V1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_V1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_V1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_V1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_V2(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_V2_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_V2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_V2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_V3(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_V3_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_V3_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_SV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_SV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_DV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_DV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_DV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_DV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_DV2(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_DV2_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_DV2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_DV2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 psArray *pmSourcesReadCMP (char *filename, psMetadata *header);

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c

-              r29004
+              r30631
     return true;
+}
+bool pmSourcesWrite_CMF_PS1_DV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    return true;
+}

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c

-              r29620
+              r30631
     return true;
+}
+bool pmSourcesWrite_CMF_PS1_DV2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    return true;
+}

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

-              r29019
+              r30631
 #include "pmConfig.h"
+#include "pmErrorCodes.h"
 #include "pmDetrendDB.h"
 …
     psF32 errMag, chisq, apRadius;
     psS32 nPix, nDOF;
-    char keyword1[80], keyword2[80];
     pmChip *chip = readout->parent->parent;
 …
     table = psArrayAllocEmpty (sources->n);
+# if (0)
     // we use this just to define the output vectors (which must be present for all objects)
     bool status = false;
 …
       psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
+    }
+# endif
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
 …
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
+# if (0)
+        // XXX if we have raw radial apertures, write them out here
         psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
         psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
 …
         psFree (radFluxErr);
         psFree (radFill);
+# endif
         // XXX not sure how to get this : need to load Nimages with weight?
 …
     return true;
+}
+// **** write out the radial flux values for the sources for a given matched-PSF image
+// **** how do we distinguish the matched-PSF images from the non-matched version
+bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    psArray *table;
+    psMetadata *row;
+    psF32 xPos, yPos;
+    char keyword1[80], keyword2[80];
+    // create a header to hold the output data
+    psMetadata *outhead = psMetadataAlloc ();
+    // write the links to the image header
+    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "radial flux table extension", extname);
+    // we use this just to define the output vectors (which must be present for all objects)
+    bool status = false;
+    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
+    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    psAssert (radMax, "this must have been defined and tested earlier!");
+    psAssert (radMax->n, "this must have been defined and tested earlier!");
+    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
+    // write the radial profile apertures to header
+    for (int i = 0; i < radMax->n; i++) {
+      sprintf (keyword1, "RMIN_%02d", i);
+      sprintf (keyword2, "RMAX_%02d", i);
+      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
+      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
+    }
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    if (!fwhmValues) {
+        psError (PM_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
+        return false;
+    }
+    // let's write these out in S/N order
+    sources = psArraySort (sources, pmSourceSortBySN);
+    table = psArrayAllocEmpty (sources->n);
+    // we write out all sources, regardless of quality.  the source flags tell us the state
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        // skip sources without radial aper measurements (or insufficient)
+        if (source->radialAper == NULL) continue;
+        psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
+        for (int entry = 0; entry < fwhmValues->n; entry++) {
+            // choose the convolved EXT model, if available, otherwise the simple one
+            pmSourceRadialApertures *radialAper = source->radialAper->data[entry];
+            assert (radialAper);
+            bool useMoments = true;
+            useMoments = (useMoments && source->moments);          // can't if there are no moments
+            useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
+            useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
+            if (useMoments) {
+                xPos = source->moments->Mx;
+                yPos = source->moments->My;
+            } else {
+                xPos = source->peak->xf;
+                yPos = source->peak->yf;
+            }
+            row = psMetadataAlloc ();
+            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
+            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                        fwhmValues->data.F32[entry]);
+            // XXX if we have raw radial apertures, write them out here
+            psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVectorInit (radFlux,    NAN);
+            psVectorInit (radFluxErr, NAN);
+            psVectorInit (radFill,    NAN);
+            if (!radialAper->flux) goto write_annuli;
+            if (!radialAper->fill) goto write_annuli;
+            psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths");
+            psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths");
+            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
+            for (int j = 0; j < radialAper->flux->n; j++) {
+                radFlux->data.F32[j]    = radialAper->flux->data.F32[j];
+                radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j];
+                radFill->data.F32[j]    = radialAper->fill->data.F32[j];
+            }
+        write_annuli:
+            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",     PS_DATA_VECTOR, "flux within annuli",    radFlux);
+            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli",  radFluxErr);
+            psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",     PS_DATA_VECTOR, "fill factor of annuli", radFill);
+            psFree (radFlux);
+            psFree (radFluxErr);
+            psFree (radFill);
+            psArrayAdd (table, 100, row);
+            psFree (row);
+        }
+    }
+    if (table->n == 0) {
+        if (!psFitsWriteBlank (fits, outhead, extname)) {
+            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
+            psFree(outhead);
+            psFree(table);
+            return false;
+        }
+        psFree (outhead);
+        psFree (table);
+        return true;
+    }
+    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
+    if (!psFitsWriteTable (fits, outhead, table, extname)) {
+        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
+    }
+    psFree (outhead);
+    psFree (table);
+    return true;
+}

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c

-              r29004
+              r30631
     return false;
+}
+bool pmSourcesWrite_CMF_PS1_V1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    return true;
+}

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

-              r29004
+              r30631
     return true;
+}
+bool pmSourcesWrite_CMF_PS1_V2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    return true;
+}

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c

-              r29620
+              r30631
     return true;
+}
+bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    return true;
+}

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_PS1_CAL_0.c

-              r29004
+              r30631
     return true;
+}
+bool pmSourcesWrite_PS1_CAL_0_XRAD (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    return true;
+}

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_PS1_DEV_0.c

-              r29004
+              r30631
     return true;
+}
+bool pmSourcesWrite_PS1_DEV_0_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    return true;
+}

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

-              r29004
+              r30631
     return true;
+}
+bool pmSourcesWrite_PS1_DEV_1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    return true;
+}

branches/czw_branch/20101203/psModules/src/objects/pmSourceIO_SMPDATA.c

-              r29004
+              r30631
     return true;
+}
+bool pmSourcesWrite_SMPDATA_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+    return true;
+}

branches/czw_branch/20101203/psModules/src/psmodules.h

r29388	r30631
98	98	// the following headers are from psModule:imcombine
99	99	#include <pmStack.h>
	100	#include <pmSubtractionTypes.h>
100	101	#include <pmStackReject.h>
101	102	#include <pmSubtraction.h>

Context Navigation

Legend:

Download in other formats: