Changeset 27876 for branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFsUtils.c

Timestamp:

May 6, 2010, 8:02:38 PM (16 years ago)

Author:

eugene

Message:

working on psphotStack

File:

• branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFsUtils.c (modified) (12 diffs)

branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFsUtils.c

@@ -1 @@
-/***** defines *****/
-
-#define ARRAY_BUFFER 16                 // Number to add to array at a time
-#define MAG_IGNORE 50                   // Ignore magnitudes fainter than this --- they're not real!
-#define FAKE_SIZE 1                     // Size of fake convolution kernel
-#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
-#define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
-#define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
+# include "psphotInternal.h"
+# define ARRAY_BUFFER 16                 // Number to add to array at a time
 
 // XXX better name
@@ -18 +12 @@
     psFree(resolved);
     if (!fits) {
-        psError(PPSTACK_ERR_IO, false, "Unable to open previously produced image: %s", name);
+        psError(PSPHOT_ERR_IO, false, "Unable to open previously produced image: %s", name);
         return false;
     }
     psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
     if (!image) {
-        psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name);
+        psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image: %s", name);
         psFitsClose(fits);
         return false;
@@ -52 +46 @@
 
 psArray *stackSourcesFilter(psArray *sources, // Source list to filter
-                                   int exclusion // Exclusion zone, pixels
+                            int exclusion // Exclusion zone, pixels
     )
 {
@@ -90 +84 @@
 
         long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
-        psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
+        psTrace("psphotStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
                 coords->data.F64[0], coords->data.F64[1], numWithin);
         if (numWithin == 1) {
@@ -104 +98 @@
     psFree(y);
 
-    psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
+    psLogMsg("psphotStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
 
     return filtered;
@@ -111 +105 @@
 // Add background into the fake image
 // Based on ppSubBackground()
-static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model
+psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model
                                      const pmConfig *config // Configuration
     )
@@ -121 +115 @@
     int numCols = image->numCols, numRows = image->numRows; // Size of image
 
-    psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE);
+    psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, "PPSTACK");
     psAssert(ppStackRecipe, "Need PPSTACK recipe");
-    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
+    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, "PSPHOT");
     psAssert(psphotRecipe, "Need PSPHOT recipe");
 
@@ -138 +132 @@
     psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model
     if (!psImageUnbin(unbinned, binned, binning)) {
-        psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model");
+        psError(PSPHOT_ERR_DATA, false, "Unable to unbin background model");
         psFree(binned);
         psFree(unbinned);
@@ -153 +147 @@
     )
 {
-#if 1
     bool mdok; // Status of metadata lookups
 
-    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack
+    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, "PPSTACK"); // Recipe for ppStack
     psAssert(recipe, "Need PPSTACK recipe");
 
@@ -163 +156 @@
     int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM");
     if (!mdok) {
-        psError(PPSTACK_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe");
+        psError(PSPHOT_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe");
         return false;
     }
     float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN");
     if (!mdok) {
-        psError(PPSTACK_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe");
+        psError(PSPHOT_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe");
         return false;
     }
     float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX");
     if (!mdok) {
-        psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe");
+        psError(PSPHOT_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe");
         return false;
     }
@@ -181 +174 @@
     psImageCovarianceTransfer(readout->variance, readout->covariance);
     return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
-#else
-    return true;
-#endif
 }
 
@@ -190 +180 @@
 // It implicitly assumes the output root name is the same between invocations.
 
-bool loadKernel () {
-            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
-            psAssert(file, "Require file");
-
-            pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
-            view->chip = view->cell = view->readout = 0;
-            psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest
-
-            // Read convolution kernel
-            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
-            psFree(filename);
-            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
-            psFree(resolved);
-            if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) {
-                psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel");
-                psFitsClose(fits);
-                return false;
-            }
-            psFitsClose(fits);
-
-            if (!readImage(&readout->image, options->convImages->data[index], config) ||
-                !readImage(&readout->mask, options->convMasks->data[index], config) ||
-                !readImage(&readout->variance, options->convVariances->data[index], config)) {
-                psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image.");
-                return false;
-            }
-
-            psRegion *region = psMetadataLookupPtr(NULL, conv->analysis,
-                                                   PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
-            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis,
-                                                                PM_SUBTRACTION_ANALYSIS_KERNEL);
-
-            pmSubtractionAnalysis(conv->analysis, NULL, kernels, region,
-                                  readout->image->numCols, readout->image->numRows);
-
-            psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
-            bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
-            psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
-            psImageCovarianceSetThreads(oldThreads);
-            psFree(readout->covariance);
-            readout->covariance = covar;
-            psFree(kernel);
-}
-
-bool dumpImage() {
-    // XXX should be optional
-            {
-                pmHDU *hdu = pmHDUFromCell(readout->parent);
-                psString name = NULL;
-                psStringAppend(&name, "fake_%03d.fits", index);
-                pmStackVisualPlotTestImage(fake->image, name);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
-                psFitsClose(fits);
-            }
-            {
-                pmHDU *hdu = pmHDUFromCell(readout->parent);
-                psString name = NULL;
-                psStringAppend(&name, "real_%03d.fits", index);
-                pmStackVisualPlotTestImage(readout->image, name);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
-                psFitsClose(fits);
-            }
-}
-
-bool dumpImage2() {
-    // XXX should be optional
-
-            {
-                pmHDU *hdu = pmHDUFromCell(readout->parent);
-                psString name = NULL;
-                psStringAppend(&name, "conv_%03d.fits", index);
-                pmStackVisualPlotTestImage(conv->image, name);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL);
-                psFitsClose(fits);
-            }
-            {
-                pmHDU *hdu = pmHDUFromCell(readout->parent);
-                psString name = NULL;
-                psStringAppend(&name, "diff_%03d.fits", index);
-                pmStackVisualPlotTestImage(fake->image, name);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psBinaryOp(fake->image, conv->image, "-", fake->image);
-                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
-                psFitsClose(fits);
-            }
-}
-
-bool dumpImage3() 
+bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index) {
+
+    // Read the convolution kernel from the saved file
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
+    psAssert(file, "Require file");
+
+    pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
+    view->chip = view->cell = view->readout = 0;
+    psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest
+
+    // Read convolution kernel data
+    psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
+    psFree(filename);
+    psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
+    psFree(resolved);
+    if (!fits || !pmReadoutReadSubtractionKernels(readoutCnv, fits)) {
+        psError(PSPHOT_ERR_IO, false, "Unable to read previously produced kernel");
+        psFitsClose(fits);
+        return false;
+    }
+    psFitsClose(fits);
+
+    // read the convolved pixels (image, mask, variance) -- names are pre-defined
+    if (!readImage(&readoutCnv->image,    options->convImages->data[index],    config) ||
+        !readImage(&readoutCnv->mask,     options->convMasks->data[index],     config) ||
+        !readImage(&readoutCnv->variance, options->convVariances->data[index], config)) {
+        psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image.");
+        return false;
+    }
+
+    // XXX ??? not sure what is happening here -- consult Paul
+    psRegion *region = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
+    pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
+
+    pmSubtractionAnalysis(readoutCnv->analysis, NULL, kernels, region, readoutCnv->image->numCols, readoutCnv->image->numRows);
+
+    psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
+
+    // update the covariance matrix 
+    // XXX why is this needed if we have correctly read the saved data?
+    bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
+    psKernel *covar = psImageCovarianceCalculate(kernel, readoutCnv->covariance); // Covariance matrix
+    psImageCovarianceSetThreads(oldThreads);
+    psFree(readoutCnv->covariance);
+    readoutCnv->covariance = covar;
+    psFree(kernel);
+    return true;
+}
+
+bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname) {
+
+    pmHDU *hdu = pmHDUFromCell(readoutRef->parent);
+    psString name = NULL;
+    psStringAppend(&name, "%s_%03d.fits", rootname, index);
+    pmStackVisualPlotTestImage(readoutOut->image, name);
+    psFits *fits = psFitsOpen(name, "w");
+    psFree(name);
+    psFitsWriteImage(fits, hdu->header, readoutOut->image, 0, NULL);
+    psFitsClose(fits);
+    return true;
+}
+
+bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname) {
+
+    pmHDU *hdu = pmHDUFromCell(readoutRef->parent);
+    psString name = NULL;
+    psStringAppend(&name, "%s_%03d.fits", rootname, index);
+    pmStackVisualPlotTestImage(readoutFake->image, name);
+    psFits *fits = psFitsOpen(name, "w");
+    psFree(name);
+    psBinaryOp(readoutFake->image, readoutConv->image, "-", readoutFake->image);
+    psFitsWriteImage(fits, hdu->header, readoutFake->image, 0, NULL);
+    psFitsClose(fits);
+    return true;
+}
+
+// perform the bulk of the PSF-matching
+bool matchKernel(pmConfig *config, pmReadout *readoutCnv, pmReadout *readoutRaw, psphotStackOptions *options, int index) {
+
+    bool mdok;
+
+    psAssert(options->psf, "Require target PSF");
+    psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
+
+    psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
+    psAssert(stackRecipe, "We've thrown an error on this before.");
+
+    // Look up appropriate values from the ppSub recipe
+    psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    psAssert(subRecipe, "recipe missing");
+
+    psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
+    psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor
+    psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
+
+    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
+    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+    float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness
+    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
+
+    int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order
+    float regionSize = psMetadataLookupF32(NULL, subRecipe, "REGION.SIZE"); // Size of iso-kernel regs
+    float spacing = psMetadataLookupF32(NULL, subRecipe, "STAMP.SPACING"); // Typical stamp spacing
+
+    int footprint = psMetadataLookupS32(NULL, subRecipe, "STAMP.FOOTPRINT"); // Stamp half-size
+    int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size
+
+    float threshold = psMetadataLookupF32(NULL, subRecipe, "STAMP.THRESHOLD"); // Threshold for stmps
+    int stride = psMetadataLookupS32(NULL, subRecipe, "STRIDE"); // Size of convolution patches
+    int iter = psMetadataLookupS32(NULL, subRecipe, "ITER"); // Rejection iterations
+    float rej = psMetadataLookupF32(NULL, subRecipe, "REJ"); // Rejection threshold
+    float kernelError = psMetadataLookupF32(NULL, subRecipe, "KERNEL.ERR"); // Relative systematic error in kernel
+    float normFrac = psMetadataLookupF32(NULL, subRecipe, "NORM.FRAC"); // Fraction of window for normalisn windw
+    float sysError = psMetadataLookupF32(NULL, subRecipe, "SYS.ERR"); // Relative systematic error in images
+    float skyErr = psMetadataLookupF32(NULL, subRecipe, "SKY.ERR"); // Additional error in sky
+    float covarFrac = psMetadataLookupF32(NULL, subRecipe, "COVAR.FRAC"); // Fraction for covariance calculation
+
+    const char *typeStr = psMetadataLookupStr(NULL, subRecipe, "KERNEL.TYPE"); // Kernel type
+    pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
+    psVector *widths = psMetadataLookupPtr(NULL, subRecipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
+    psVector *orders = psMetadataLookupPtr(NULL, subRecipe, "ISIS.ORDERS"); // ISIS Polynomial orders
+    int inner = psMetadataLookupS32(NULL, subRecipe, "INNER"); // Inner radius
+    int ringsOrder = psMetadataLookupS32(NULL, subRecipe, "RINGS.ORDER"); // RINGS polynomial order
+    int binning = psMetadataLookupS32(NULL, subRecipe, "SPAM.BINNING"); // Binning for SPAM kernel
+    float badFrac = psMetadataLookupF32(NULL, subRecipe, "BADFRAC"); // Maximum bad fraction
+    float optThresh = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.TOL"); // Tolerance for search
+    int optOrder = psMetadataLookupS32(&mdok, subRecipe, "OPTIMUM.ORDER"); // Order for search
+    float poorFrac = psMetadataLookupF32(&mdok, subRecipe, "POOR.FRACTION"); // Fraction for "poor"
+
+    bool scale = psMetadataLookupBool(NULL, subRecipe, "SCALE");        // Scale kernel parameters?
+    float scaleRef = psMetadataLookupF32(NULL, subRecipe, "SCALE.REF"); // Reference for scaling
+    float scaleMin = psMetadataLookupF32(NULL, subRecipe, "SCALE.MIN"); // Minimum for scaling
+    float scaleMax = psMetadataLookupF32(NULL, subRecipe, "SCALE.MAX"); // Maximum for scaling
+    if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
+        psError(PSPHOT_ERR_CONFIG, false,
+                "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.",
+                scaleRef, scaleMin, scaleMax);
+        return false;
+    }
+
+    // These values are specified specifically for stacking
+    const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename
+
+    psVector *widthsCopy = NULL;
+    psVector *optWidths = NULL;
+    pmReadout *fake = NULL;
+    psArray *stampSources = NULL;
+
+    bool optimum = false;
+    optWidths = SetOptWidths(&optimum, subRecipe); // Vector with FWHMs for optimum search
+
+    // For the sake of stamps, remove nearby sources
+    stampSources = stackSourcesFilter(options->sourceLists->data[index], footprint); // Filtered list of sources
+
+    fake = makeFakeReadout(config, readoutRaw, stampSources, options->psf, maskVal | maskBad, footprint + size);
+    if (!fake) goto escape;
+
+    dumpImage(fake, readoutRaw, index, "fake");
+    dumpImage(readoutRaw,  readoutRaw, index, "real");
+
+    if (threads) pmSubtractionThreadsInit();
+
+    // Do the image matching
+    pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readoutRaw->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
+    if (kernel) {
+        if (!pmSubtractionMatchPrecalc(NULL, readoutCnv, fake, readoutRaw, readoutRaw->analysis, stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac)) {
+            psError(psErrorCodeLast(), false, "Unable to convolve images.");
+            goto escape;
+        }
+    } else {
+        // Scale the input parameters
+        widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
+        if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, options->inputSeeing->data.F32[index], options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
+            psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
+            goto escape;
+        }
+
+        if (!pmSubtractionMatch(NULL, readoutCnv, fake, readoutRaw, footprint, stride, regionSize, spacing, threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
+            psError(psErrorCodeLast(), false, "Unable to match images.");
+            goto escape;
+        }
+    }
+
+    // Reject image completely if the maximum deconvolution fraction exceeds the limit
+    float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
+    float deconv = psMetadataLookupF32(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
+    if (deconv > deconvLimit) {
+        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index);
+        goto escape;
+    }
+
+    dumpImage(readoutCnv, readoutRaw, index, "conv");
+    dumpImageDiff(readoutCnv, fake, readoutRaw, index, "diff");
+
+    psFree(fake);
+    psFree(optWidths);
+    psFree(stampSources);
+    psFree(widthsCopy);
+    pmSubtractionThreadsFinalize();
+    return true;
+
+escape:
+    psFree(fake);
+    psFree(optWidths);
+    psFree(stampSources);
+    psFree(widthsCopy);
+    pmSubtractionThreadsFinalize();
+    return false;
+}
+
+// Extract the regions and solutions used in the image matching
+// This stops them from being freed when we iterate back up the FPA
+// Record the chi-square value
+// XXX this function may not be needed for psphotStack
+bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index) {
+
+    psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions
     {
-        pmHDU *hdu = pmHDUFromCell(readout->parent);
-        psString name = NULL;
-        psStringAppend(&name, "convolved_%03d.fits", index);
-        pmStackVisualPlotTestImage(readout->image, name);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
-        psFitsClose(fits);
-    }
-
-bool matchKernel() {
-            // Normal operations here
-            psAssert(options->psf, "Require target PSF");
-            psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
-
-            int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order
-            float regionSize = psMetadataLookupF32(NULL, subRecipe, "REGION.SIZE"); // Size of iso-kernel regs
-            float spacing = psMetadataLookupF32(NULL, subRecipe, "STAMP.SPACING"); // Typical stamp spacing
-            int footprint = psMetadataLookupS32(NULL, subRecipe, "STAMP.FOOTPRINT"); // Stamp half-size
-            float threshold = psMetadataLookupF32(NULL, subRecipe, "STAMP.THRESHOLD"); // Threshold for stmps
-            int stride = psMetadataLookupS32(NULL, subRecipe, "STRIDE"); // Size of convolution patches
-            int iter = psMetadataLookupS32(NULL, subRecipe, "ITER"); // Rejection iterations
-            float rej = psMetadataLookupF32(NULL, subRecipe, "REJ"); // Rejection threshold
-            float kernelError = psMetadataLookupF32(NULL, subRecipe, "KERNEL.ERR"); // Relative systematic error in kernel
-            float normFrac = psMetadataLookupF32(NULL, subRecipe, "NORM.FRAC"); // Fraction of window for normalisn windw
-            float sysError = psMetadataLookupF32(NULL, subRecipe, "SYS.ERR"); // Relative systematic error in images
-            float skyErr = psMetadataLookupF32(NULL, subRecipe, "SKY.ERR"); // Additional error in sky
-            float covarFrac = psMetadataLookupF32(NULL, subRecipe, "COVAR.FRAC"); // Fraction for covariance calculation
-
-            const char *typeStr = psMetadataLookupStr(NULL, subRecipe, "KERNEL.TYPE"); // Kernel type
-            pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
-            psVector *widths = psMetadataLookupPtr(NULL, subRecipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
-            psVector *orders = psMetadataLookupPtr(NULL, subRecipe, "ISIS.ORDERS"); // ISIS Polynomial orders
-            int inner = psMetadataLookupS32(NULL, subRecipe, "INNER"); // Inner radius
-            int ringsOrder = psMetadataLookupS32(NULL, subRecipe, "RINGS.ORDER"); // RINGS polynomial order
-            int binning = psMetadataLookupS32(NULL, subRecipe, "SPAM.BINNING"); // Binning for SPAM kernel
-            float badFrac = psMetadataLookupF32(NULL, subRecipe, "BADFRAC"); // Maximum bad fraction
-            bool optimum = psMetadataLookupBool(&mdok, subRecipe, "OPTIMUM"); // Derive optimum parameters?
-            float optMin = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MIN"); // Minimum width for search
-            float optMax = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MAX"); // Maximum width for search
-            float optStep = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.STEP"); // Step for search
-            float optThresh = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.TOL"); // Tolerance for search
-            int optOrder = psMetadataLookupS32(&mdok, subRecipe, "OPTIMUM.ORDER"); // Order for search
-            float poorFrac = psMetadataLookupF32(&mdok, subRecipe, "POOR.FRACTION"); // Fraction for "poor"
-
-            bool scale = psMetadataLookupBool(NULL, subRecipe, "SCALE");        // Scale kernel parameters?
-            float scaleRef = psMetadataLookupF32(NULL, subRecipe, "SCALE.REF"); // Reference for scaling
-            float scaleMin = psMetadataLookupF32(NULL, subRecipe, "SCALE.MIN"); // Minimum for scaling
-            float scaleMax = psMetadataLookupF32(NULL, subRecipe, "SCALE.MAX"); // Maximum for scaling
-            if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
-                psError(PPSTACK_ERR_CONFIG, false,
-                        "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.",
-                        scaleRef, scaleMin, scaleMax);
-                return false;
-            }
-
-
-            // These values are specified specifically for stacking
-            const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename
-
-            psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
-            if (optimum) {
-                optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
-            }
-
-            pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
-
-            psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background
-            psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
-            if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
-                psError(PPSTACK_ERR_DATA, false, "Can't measure background for image.");
-                psFree(fake);
-                psFree(optWidths);
-                psFree(conv);
-                psFree(bg);
-                psFree(rng);
-                return false;
-            }
-            float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_REGION);
+            regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
+        }
+        psFree(iter);
+    }
+
+    psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels
+    {
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_UNKNOWN);
+            pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
+            kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
+        }
+        psFree(iter);
+    }
+    psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match");
+
+    // Record chi^2
+    {
+        double sum = 0.0;           // Sum of chi^2
+        int num = 0;                // Number of measurements of chi^2
+        psString regex = NULL;      // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_UNKNOWN);
+            pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
+            sum += kernels->mean;
+            num++;
+        }
+        psFree(iter);
+        options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
+    }
+
+    return true;
+}
+
+// Kernel normalisation for convolved readout
+bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index) {
+
+    double sum = 0.0;           // Sum of chi^2
+    int num = 0;                // Number of measurements of chi^2
+    psString regex = NULL;      // Regular expression
+    psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
+    psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
+    psFree(regex);
+    psMetadataItem *item = NULL;// Item from iteration
+    while ((item = psMetadataGetAndIncrement(iter))) {
+        assert(item->type == PS_TYPE_F32);
+        float norm = item->data.F32; // Normalisation
+        sum += norm;
+        num++;
+    }
+    psFree(iter);
+    float conv = sum/num;       // Mean normalisation from convolution
+    float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
+    float renorm =  stars / conv; // Renormalisation to apply
+    psLogMsg("psphotStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", index, renorm, conv, stars);
+
+    psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
+    psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
+    return true;
+}
+
+// adjust scaling for readout (remove background, ..., determine weighting)
+bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index) {
+
+    psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
+    psAssert(stackRecipe, "We've thrown an error on this before.");
+
+    // Look up appropriate values from the ppSub recipe
+    psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    psAssert(subRecipe, "recipe missing");
+
+    psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
+    psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
+
+    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+    // Ensure the background value is zero
+    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+
+    // XXX why is this in config->arguments and not recipe?
+    if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
+        if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
+            psAbort("Can't measure background for image.");
+            // XXX we used to clear error: why is this acceptable? psErrorClear(); 
+        }
+
+        float value = psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN);
+        float stdev = psStatsGetValue(bg, PS_STAT_ROBUST_STDEV);
+
+        psLogMsg("psphotStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", value, stdev);
+        psBinaryOp(readout->image, readout->image, "-", psScalarAlloc(value, PS_TYPE_F32));
+    }
+
+    if (!stackRenormaliseReadout(config, readout)) {
+        psFree(rng);
+        psFree(bg);
+        return false;
+    }
+
+    // Measure the variance level for the weighting
+    if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) {
+        if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
+            psError(PSPHOT_ERR_DATA, false, "Can't measure mean variance for image.");
             psFree(rng);
             psFree(bg);
-
-            // For the sake of stamps, remove nearby sources
-            psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index],
-                                                       footprint); // Filtered list of sources
-
-            bool oldThreads = pmReadoutFakeThreads(true); // Old threading state
-            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
-                                          stampSources, SOURCE_MASK, NULL, NULL, options->psf,
-                                          minFlux, footprint + size, false, true)) {
-                psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF.");
-                psFree(fake);
-                psFree(optWidths);
-                psFree(conv);
-                return false;
-            }
-            pmReadoutFakeThreads(oldThreads);
-
-            fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);
-
-            // Add the background into the target image
-            psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
-            psBinaryOp(fake->image, fake->image, "+", bgImage);
-            psFree(bgImage);
-
-            dumpImage();
-
-            if (threads > 0) {
-                pmSubtractionThreadsInit();
-            }
-
-            // Do the image matching
-            pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readout->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
-            if (kernel) {
-                if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis,
-                                               stride, kernelError, covarFrac, maskVal, maskBad, maskPoor,
-                                               poorFrac, badFrac)) {
-                    psError(psErrorCodeLast(), false, "Unable to convolve images.");
-                    psFree(fake);
-                    psFree(optWidths);
-                    psFree(stampSources);
-                    psFree(conv);
-                    if (threads > 0) {
-                        pmSubtractionThreadsFinalize();
-                    }
-                    return false;
-                }
-            } else {
-                // Scale the input parameters
-                psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
-                if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy,
-                                                       options->inputSeeing->data.F32[index],
-                                                       options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
-                    psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
-                    psFree(fake);
-                    psFree(optWidths);
-                    psFree(stampSources);
-                    psFree(conv);
-                    psFree(widthsCopy);
-                    if (threads > 0) {
-                        pmSubtractionThreadsFinalize();
-                    }
-                    return false;
-                }
-
-                if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
-                                        threshold, stampSources, stampsName, type, size, order, widthsCopy,
-                                        orders, inner, ringsOrder, binning, penalty,
-                                        optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
-                                        sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor,
-                                        poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
-                    psError(psErrorCodeLast(), false, "Unable to match images.");
-                    psFree(fake);
-                    psFree(optWidths);
-                    psFree(stampSources);
-                    psFree(conv);
-                    psFree(widthsCopy);
-                    if (threads > 0) {
-                        pmSubtractionThreadsFinalize();
-                    }
-                    return false;
-                }
-                psFree(widthsCopy);
-            }
-
-            dumpImage2();
-
-            psFree(fake);
-            psFree(optWidths);
-            psFree(stampSources);
-
-            if (threads > 0) {
-                pmSubtractionThreadsFinalize();
-            }
-
-            // Replace original images with convolved
-            psFree(readout->image);
-            psFree(readout->mask);
-            psFree(readout->variance);
-            psFree(readout->covariance);
-            readout->image  = psMemIncrRefCounter(conv->image);
-            readout->mask   = psMemIncrRefCounter(conv->mask);
-            readout->variance = psMemIncrRefCounter(conv->variance);
-            readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC);
-
-}
-
-bool saveMatchData () {
-       // Extract the regions and solutions used in the image matching
-        // This stops them from being freed when we iterate back up the FPA
-        psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions
-        {
-            psString regex = NULL;          // Regular expression
-            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
-            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
-            psFree(regex);
-            psMetadataItem *item = NULL;// Item from iteration
-            while ((item = psMetadataGetAndIncrement(iter))) {
-                assert(item->type == PS_DATA_REGION);
-                regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
-            }
-            psFree(iter);
+            return false;
         }
-        psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels
-        {
-            psString regex = NULL;          // Regular expression
-            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
-            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
-            psFree(regex);
-            psMetadataItem *item = NULL;// Item from iteration
-            while ((item = psMetadataGetAndIncrement(iter))) {
-                assert(item->type == PS_DATA_UNKNOWN);
-                pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
-                kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
-            }
-            psFree(iter);
-        }
-        psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match");
-}
-
-bool saveChiSquare() {
-        // Record chi^2
-        {
-            double sum = 0.0;           // Sum of chi^2
-            int num = 0;                // Number of measurements of chi^2
-            psString regex = NULL;      // Regular expression
-            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
-            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
-            psFree(regex);
-            psMetadataItem *item = NULL;// Item from iteration
-            while ((item = psMetadataGetAndIncrement(iter))) {
-                assert(item->type == PS_DATA_UNKNOWN);
-                pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
-                sum += kernels->mean;
-                num++;
-            }
-            psFree(iter);
-            options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
-        }
-
-}
-
-bool renormKernel() {
-        // Kernel normalisation
-        {
-            double sum = 0.0;           // Sum of chi^2
-            int num = 0;                // Number of measurements of chi^2
-            psString regex = NULL;      // Regular expression
-            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
-            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
-            psFree(regex);
-            psMetadataItem *item = NULL;// Item from iteration
-            while ((item = psMetadataGetAndIncrement(iter))) {
-                assert(item->type == PS_TYPE_F32);
-                float norm = item->data.F32; // Normalisation
-                sum += norm;
-                num++;
-            }
-            psFree(iter);
-            float conv = sum/num;       // Mean normalisation from convolution
-            float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
-            float renorm =  stars / conv; // Renormalisation to apply
-            psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n",
-                     index, renorm, conv, stars);
-            psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
-            psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
-        }
-
-}
+        options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance));
+    } else {
+        options->weightings->data.F32[index] = 1.0;
+    }
+    psLogMsg("psphotStack", PS_LOG_INFO, "Weighting for image %d is %f\n", index, options->weightings->data.F32[index]);
+
+    psFree(rng);
+    psFree(bg);
+    return true;
+}
+
+# define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
+# define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
+
+// generate a fake readout against which to PSF match
+pmReadout *makeFakeReadout(pmConfig *config, pmReadout *readoutRaw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize) {
+
+    pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
+
+    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+    if (!psImageBackground(bg, NULL, readoutRaw->image, readoutRaw->mask, maskVal, rng)) {
+        psError(PSPHOT_ERR_DATA, false, "Can't measure background for image.");
+        psFree(fake);
+        psFree(bg);
+        psFree(rng);
+        return NULL;
+    }
+    float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image
+    psFree(rng);
+    psFree(bg);
+
+    bool oldThreads = pmReadoutFakeThreads(true); // Old threading state
+    if (!pmReadoutFakeFromSources(fake, readoutRaw->image->numCols, readoutRaw->image->numRows, sources, SOURCE_MASK, NULL, NULL, psf, minFlux, fullSize, false, true)) {
+        psError(PSPHOT_ERR_DATA, false, "Unable to generate fake image with target PSF.");
+        psFree(fake);
+        return NULL;
+    }
+    pmReadoutFakeThreads(oldThreads);
+
+    fake->mask = psImageCopy(NULL, readoutRaw->mask, PS_TYPE_IMAGE_MASK);
+
+    // Add the background into the target image
+    psImage *bgImage = stackBackgroundModel(readoutRaw, config); // Image of background
+    psBinaryOp(fake->image, fake->image, "+", bgImage);
+    psFree(bgImage);
+
+    return fake;
+}
+
+// set the widths 
+psVector *SetOptWidths (bool *optimum, psMetadata *recipe) {
+
+    bool status;
+
+    *optimum = psMetadataLookupBool(&status, recipe, "OPTIMUM"); // Derive optimum parameters?
+    psAssert (status, "missing recipe value %s", "OPTIMUM");
+
+    psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
+
+    if (*optimum) {
+        float optMin = psMetadataLookupF32(&status,  recipe, "OPTIMUM.MIN"); // Minimum width for search
+        psAssert (status, "missing recipe value %s", "OPTIMUM.MIN");
+        
+        float optMax = psMetadataLookupF32(&status,  recipe, "OPTIMUM.MAX"); // Maximum width for search
+        psAssert (status, "missing recipe value %s", "OPTIMUM.MAX");
+        
+        float optStep = psMetadataLookupF32(&status, recipe, "OPTIMUM.STEP"); // Step for search
+        psAssert (status, "missing recipe value %s", "OPTIMUM.STEP");
+
+        optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
+    }
+
+    return optWidths;
+}