Changeset 27876

Timestamp:

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

Author:

eugene

Message:

working on psphotStack

Location:

branches/eam_branches/psphot.20100506

Files:

• doc/stack.txt (modified) (1 diff)
• src/Makefile.am (modified) (1 diff)
• src/psphot.h (modified) (1 diff)
• src/psphotErrorCodes.dat (modified) (1 diff)
• src/psphotSetMaskBits.c (modified) (1 diff)
• src/psphotStackArguments.c (modified) (3 diffs)
• src/psphotStackImageLoop.c (modified) (4 diffs)
• src/psphotStackMatchPSFs.c (modified) (4 diffs)
• src/psphotStackMatchPSFsPrepare.c (added)
• src/psphotStackMatchPSFsUtils.c (modified) (12 diffs)
• src/psphotStackOptions.c (added)
• src/psphotStackParseCamera.c (modified) (3 diffs)

branches/eam_branches/psphot.20100506/doc/stack.txt

@@ +1 @@
+
+20100506:
+
+  we may load RAW and/or CNV images.
+
+  * options->convolve:
+    * if only one of RAW or CNV exists, convolve it to match target PSF
+    * if both exist, convolve desired one
+    
+    * if matching PSF exists, use it (unless told to re-measure)
+    
+    * output goes to OUT (which is then used by psphot analysis routines)
+
+  * detect
+
+    * if (RAW) ? RAW : OUT
+    
+  *   
+
 
 20100503:

branches/eam_branches/psphot.20100506/src/Makefile.am

@@ -95 +95 @@
         psphotFitSourcesLinearStack.c \
         psphotSourceMatch.c           \
+        psphotStackMatchPSFs.c       \
+        psphotStackMatchPSFsUtils.c  \
         psphotCleanup.c
 

branches/eam_branches/psphot.20100506/src/psphot.h

@@ -367 +367 @@
 int pmPhotObjSortByX (const void **a, const void **b);
 
+/// Options for stacking process
+typedef struct {
+    // Setup
+    
+    int num;                            // Number of inputs
+    bool convolve;                      // Convolve images?
+    // psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images
+
+    // bool matchZPs;                      // Adjust relative fluxes based on transparency analysis?
+    // bool photometry;                    // Perform photometry?
+    // psMetadata *stats;                  // Statistics for output
+    // FILE *statsFile;                    // File to which to write statistics
+    // psArray *origImages, *origMasks, *origVariances; // Filenames of the original images
+    // psArray *origCovars;                // Original covariances matrices
+    // int quality;                        // Bad data quality flag
+
+    // Prepare
+    pmPSF *psf;                         // Target PSF
+    psVector *inputSeeing;              // Input seeing FWHMs
+    psVector *inputMask;                // Mask for inputs
+
+    float targetSeeing;                 // Target seeing FWHM
+    psArray *sourceLists;               // Individual lists of sources for matching
+    psVector *norm;                     // Normalisation for each image
+    psArray *psfs;
+
+    // psVector *exposures;                // Exposure times
+    // float sumExposure;                  // Sum of exposure times
+    // float zp;                           // Zero point for output
+    // psVector *inputMask;                // Mask for inputs
+    // psArray *sources;                   // Matched sources
+
+    // Convolve
+    psArray *kernels;                   // PSF-matching kernels --- required in the stacking
+    psArray *regions;                   // PSF-matching regions --- required in the stacking
+    psVector *matchChi2;                // chi^2 for stamps from matching
+    psVector *weightings;               // Combination weightings for images (1/noise^2)
+    // psArray *cells;                     // Cells for convolved images --- a handle for reading again
+    // int numCols, numRows;               // Size of image
+    // psArray *convCovars;                // Convolved covariance matrices
+
+    // Combine initial
+    // pmReadout *outRO;                   // Output readout
+    // pmReadout *expRO;                   // Exposure readout
+    // psArray *inspect;                   // Array of arrays of pixels to inspect
+
+    // Rejection
+    // psArray *rejected;                  // Rejected pixels
+} psphotStackOptions;
+
+/*** psphotStackMatchPSF prototypes ***/
+bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view);
+bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
+
+// psphotStackMatchPSFsUtils
+psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
+pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
+bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index);
+bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index);
+bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index);
+bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index);
+bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
+bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
+bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
+bool stackRenormaliseReadout(const pmConfig *config, pmReadout *readout);
+psImage *stackBackgroundModel(pmReadout *ro, const pmConfig *config);
+psArray *stackSourcesFilter(psArray *sources, int exclusion);
+void coordsFromSource(float *x, float *y, const pmSource *source);
+bool readImage(psImage **target, const char *name, const pmConfig *config);
+
 #endif

branches/eam_branches/psphot.20100506/src/psphotErrorCodes.dat

@@ -11 +11 @@
 APERTURE                Problem with aperture photometry
 SKY                     Problem in sky determination
+IO                      Problem in data I/O
 # these errors correspond to standard exit conditions
 ARGUMENTS               Incorrect arguments

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

@@ -37 +37 @@
     return true;
 }
+
+// XXX should these be in config->analysis or somewhere else besides 'recipe'?

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

@@ -11 +11 @@
     if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(stdout);
     if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(stdout);
+
+    if (psArgumentGet(argc, argv, "-template")) dumpTemplate();
 
     // load config data from default locations
@@ -84 +86 @@
           "\n"
           "where INPUTS.mdc contains various METADATAs, each with:\n"
-          "\tIMAGE(STR):     Image filename\n"
-          "\tMASK(STR):      Mask filename\n"
-          "\tVARIANCE(STR):  Variance map filename\n"
+          "\tIMAGE    : Image filename\n"
+          "\tMASK     : Mask filename\n"
+          "\tVARIANCE : Variance map filename\n"
+          "(use -template to generate a sample input.mdc file)\n"
           "OUTROOT is the 'root name' for output files\n"
           "\n"
@@ -144 +147 @@
 }
 
+static void dumpTemplate(void) {
+
+    fprintf (stdout, "# this line is required for multiple INPUT blocks to be accepted\n");
+    fprintf (stdout, "INPUT MULTI\n\n");
+
+    fprintf (stdout, "# copy and repeat the following block as needed (one per input image set)\n");
+    fprintf (stdout, "# either RAW (unconvolved) or CNV (convolved) input images are required\n");
+    fprintf (stdout, "# if both are supplied, by default RAW is used for detection, CNV is convolved (further) to match target PSF\n");
+    fprintf (stdout, "# if MASK or VARIANCE images are not supplied, they will be generated\n");
+    fprintf (stdout, "# if MASK or VARIANCE images are not supplied, they will be generated\n");
+    fprintf (stdout, "# PSF may be supplied for the convolution target\n");
+    fprintf (stdout, "INPUT METADATA\n");
+    fprintf (stdout, "  RAW:IMAGE     STR   file.im.fits   # signal image filename\n");
+    fprintf (stdout, "  RAW:MASK      STR   file.mk.fits   # mask image filename\n");
+    fprintf (stdout, "  RAW:VARIANCE  STR   file.wt.fits   # variance image filename\n");
+    fprintf (stdout, "  RAW:PSF       STR   file.psf.fits  # psf from input unconvolved image\n");
+
+    fprintf (stdout, "  CNV:IMAGE     STR   file.im.fits   # signal image filename\n");
+    fprintf (stdout, "  CNV:MASK      STR   file.mk.fits   # mask image filename\n");
+    fprintf (stdout, "  CNV:VARIANCE  STR   file.wt.fits   # variance image filename\n");
+    fprintf (stdout, "  CNV:PSF       STR   file.psf.fits  # psf from input convolved image\n");
+
+    fprintf (stdout, "  SOURCES       STR   file.cmf       # measured source positions\n");
+    fprintf (stdout, "END\n");
+
+    psLibFinalize();
+    exit(PS_EXIT_SUCCESS);
+}

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

@@ -14 +14 @@
     pmReadout *readout;
 
-    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
-    if (!status) {
+    pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
+    pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
+
+    pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+
+    if (!input) {
         psError(PSPHOT_ERR_PROG, false, "Can't find input data!");
         return false;
@@ -29 +33 @@
         psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
         if (! chip->process || ! chip->file_exists) { continue; }
-        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphotStack.");
+        // if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphotStack.");
 
         // there is now only a single chip (multiple readouts?). loop over it and process
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
             psLogMsg ("psphot", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
-            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Cell in psphotStack.");
+            // if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Cell in psphotStack.");
 
             // process each of the readouts
@@ -41 +45 @@
                 if (! readout->data_exists) { continue; }
 
-# if (0)                
-                // uncomment to generate matched psfs
+                // PSF matching
                 if (!psphotStackMatchPSFs (config, view)) {
                     psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
@@ -48 +51 @@
                     return false;
                 }
-# endif
 
                 // XXX for now, we assume there is only a single chip in the PHU:

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

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view, bool firstPass)
+bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view)
 {
     bool status = true;
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PSPHOT"); 
+    psAssert(recipe, "We've thrown an error on this before.");
 
     int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
-    // loop over the available readouts
+    // 'options' carries info needed to perform the stack matching
+    psphotStackOptions *options = psphotStackOptionsAlloc(num);
+
+    options->convolve = psMetadataLookupBool (&status, recipe, "PSPHOT.STACK.MATCH.PSF");
+    psAssert (status, "PSPHOT.STACK.MATCH.PSF not in recipe");
+
+    if (options->convolve) {
+        char *convolveSource = psMetadataLookupStr (&status, recipe, "PSPHOT.STACK.MATCH.PSF.SOURCE");
+        options->convolveSource = psphotStackConvolveSourceFromString (convolveSource);
+        if (options->convolveSource == PSPHOT_CNV_SRC_NONE) {
+            psError (PSPHOT_ERR_CONFIG, true, "stack convolution source not defined in recipe");
+            return false;
+        }
+    }
+
+    // loop over the available readouts (ignore chisq image)
     for (int i = 0; i < num; i++) {
-        if (!psphotStackMatchPSFsReadout (config, view, i)) {
+        if (!psphotStackMatchPSFsPrepare (config, view, options, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to set PSF matching options for entry %d", i);
+            return false;
+        }
+    }
+
+    // Generate target PSF
+    if (options->convolve) {
+        options->psf = ppStackPSF(config, options->numCols, options->numRows, options->psfs, options->inputMask);
+        if (!options->psf) {
+            psError(psErrorCodeLast(), false, "Unable to determine output PSF.");
+            return false;
+        }
+        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
+        options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target
+        psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
+
+        // XXX is this needed to supply the psf to psphot??
+        // pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip
+        // psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "Target PSF", options->psf);
+        // outChip->data_exists = true;
+    }
+
+    // loop over the available readouts (ignore chisq image)
+    for (int i = 0; i < num; i++) { 
+        if (!psphotStackMatchPSFsReadout (config, view, options, i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
             return false;
@@ -19 +62 @@
 
 // convolve the image to match desired PSF
-bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, int index) {
-
-    bool status;
-    int pass;
-    float NSIGMA_PEAK = 25.0;
-    int NMAX = 0;
+bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) {
 
     // find the currently selected readout
-    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
+    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", index); // File of interest
     psAssert (fileRaw, "missing file?");
 
-    pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT.CONV", index); // File of interest
+    pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", index); // File of interest
     psAssert (fileCnv, "missing file?");
 
@@ -39 +77 @@
     psAssert (readoutCnv, "missing readout?");
 
-    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
-    psAssert (maskVal, "missing mask value?");
-
-    /***** set up recipe options *****/
-
-    psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // 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");
-
-    int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size
-
-    float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
-
-    psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
-    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
-    psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor
-    psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
-    psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
-    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
-
-    bool mdok;                          // Status of MD lookup
-    float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness
-    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
-
-    if (!pmReadoutMaskNonfinite(readout, maskVal)) {
+    // set NAN pixels to 'SAT'
+    psImageMaskType maskVal = pmConfigMaskGet("SAT", config);
+    if (!pmReadoutMaskNonfinite(readoutRaw, maskVal)) {
         psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout.");
         return false;
@@ -74 +86 @@
     // Image Matching (PSFs or just flux)
     if (options->convolve) {
-      // Full match of PSFs
-        pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily
-
-        // Read previously produced kernel
-        if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
-          loadKernel();
-        } else {
-          matchKernel();
-        } // !DEBUG.STACK
-
-        saveMatchData();
-
-        saveChiSquare();
-
-        renormKernel();
-
-        // Reject image completely if the maximum deconvolution fraction exceeds the limit
-        float deconv = psMetadataLookupF32(NULL, conv->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);
-            psFree(conv);
-            return NULL;
-        }
-
-        readout->analysis = psMetadataCopy(readout->analysis, conv->analysis);
-
-        psFree(conv);
+        matchKernel(config, readoutCnv, readoutRaw, options, index);
+        saveMatchData(readoutCnv, options, index);
+        // renormKernel(readoutCnv, options, index);
     } else {
-        // only match the flux
-        float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
-        psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
-        psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
+        // only match the flux (NO! not for multi-filter, at least!)
+        // XXX do not generate readoutCnv in this case?
+        // float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
+        // psBinaryOp(readoutRaw->image, readoutRaw->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
+        // psBinaryOp(readoutRaw->variance, readoutRaw->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
     }
 
-    // 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
-    if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
-        psWarning("Can't measure background for image.");
-        psErrorClear();
-    } else {
-        if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
-            psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
-                     psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
-            (void)psBinaryOp(readout->image, readout->image, "-",
-                             psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
-        }
-    }
+    rescaleData(readoutCnv, config, options, index);
 
-    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(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image.");
-            psFree(rng);
-            psFree(bg);
-            return false;
-        }
-        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("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n",
-             index, options->weightings->data.F32[index]);
-
-    psFree(rng);
-    psFree(bg);
-
-    dumpImage3();
+    dumpImage(readoutCnv, readoutRaw, index, "convolved");
 
     return true;
 }
+
+
+# if (0)
+// Read previously produced kernel
+if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
+    loadKernel(config, readoutCnv, options, index);
+} else {
+    matchKernel(config, readoutCnv, readoutRaw, options, index);
+}
+# endif

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;
+}

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

@@ -25 +25 @@
         psMetadata *input = item->data.md; // The input metadata of interest
 
-        // look for 'IMAGE', 'MASK', 'VARIANCE' in folder (only 'IMAGE' is required)
-
-        psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image
-        if (!image || strlen(image) == 0) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name);
+        pmFPAfile *rawInputFile = NULL;
+        pmFPAfile *cnvInputFile = NULL;
+
+        // RAW (unconvolved) input data (RAW:IMAGE, RAW:MASK, RAW:VARIANCE, RAW:PSF)
+        psString rawImage = psMetadataLookupStr(NULL, input, "RAW:IMAGE"); // Name of image
+        if (rawImage && strlen(rawImage) > 0) {
+            pmFPAfile *rawInputFile = defineFile(config, NULL, "PSPHOT.STACK.INPUT.RAW", rawImage, PM_FPA_FILE_IMAGE); // File for image
+            if (!rawInputFile) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, rawImage);
+                return false;
+            }
+            psString mask = psMetadataLookupStr(&status, input, "RAW:MASK"); // Name of mask
+            if (mask && strlen(mask) > 0) {
+                if (!defineFile(config, rawInputFile, "PSPHOT.STACK.MASK.RAW", mask, PM_FPA_FILE_MASK)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
+                    return false;
+                }
+            }
+            psString variance = psMetadataLookupStr(&status, input, "RAW:VARIANCE"); // Name of variance map
+            if (variance && strlen(variance) > 0) {
+                if (!defineFile(config, rawInputFile, "PSPHOT.STACK.VARIANCE.RAW", variance, PM_FPA_FILE_VARIANCE)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
+                    return false;
+                }
+            }
+            psString psf = psMetadataLookupStr(&status, input, "RAW:PSF"); // Name of mask
+            if (psf && strlen(psf) > 0) {
+                if (!defineFile(config, rawInputFile, "PSPHOT.STACK.PSF.RAW", psf, PM_FPA_FILE_PSF)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf);
+                    return false;
+                }
+            }
+        }
+
+        // CNV (convolved) input data (CNV:IMAGE, CNV:MASK, CNV:VARIANCE, CNV:PSF)
+        psString cnvImage = psMetadataLookupStr(NULL, input, "CNV:IMAGE"); // Name of image
+        if (cnvImage && strlen(cnvImage) > 0) {
+            pmFPAfile *cnvInputFile = defineFile(config, NULL, "PSPHOT.STACK.INPUT.CNV", cnvImage, PM_FPA_FILE_IMAGE); // File for image
+            if (!cnvInputFile) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, cnvImage);
+                return false;
+            }
+            psString mask = psMetadataLookupStr(&status, input, "CNV:MASK"); // Name of mask
+            if (mask && strlen(mask) > 0) {
+                if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.MASK.CNV", mask, PM_FPA_FILE_MASK)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
+                    return false;
+                }
+            }
+            psString variance = psMetadataLookupStr(&status, input, "CNV:VARIANCE"); // Name of variance map
+            if (variance && strlen(variance) > 0) {
+                if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.VARIANCE.CNV", variance, PM_FPA_FILE_VARIANCE)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
+                    return false;
+                }
+            }
+            psString psf = psMetadataLookupStr(&status, input, "CNV:PSF"); // Name of mask
+            if (psf && strlen(psf) > 0) {
+                if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.PSF.CNV", psf, PM_FPA_FILE_PSF)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf);
+                    return false;
+                }
+            }
+        }
+
+        if (!rawInputFile && !cnvInputFile) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s (%d) lacks IMAGE of type STR", item->name, i);
             return false;
         }
-        pmFPAfile *imageFile = defineFile(config, NULL, "PSPHOT.INPUT", image, PM_FPA_FILE_IMAGE); // File for image
-        if (!imageFile) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image);
-            return false;
-        }
-
-        psString mask = psMetadataLookupStr(&status, input, "MASK"); // Name of mask
-        if (mask && strlen(mask) > 0) {
-            if (!defineFile(config, imageFile, "PSPHOT.MASK", mask, PM_FPA_FILE_MASK)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
-                return false;
-            }
-        }
-
-        psString variance = psMetadataLookupStr(&status, input, "VARIANCE"); // Name of variance map
-        if (variance && strlen(variance) > 0) {
-            if (!defineFile(config, imageFile, "PSPHOT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
-                return false;
-            }
-        }
-        // the output sources are carried on the input->fpa structures
-        pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, imageFile, "PSPHOT.STACK.OUTPUT");
-        if (!outsources) {
-            psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT");
-            return false;
-        }
-        outsources->save = true;
-        outsources->fileID = i;         // this is used to generate output names
+
+        psString sources = psMetadataLookupStr(&status, input, "SOURCES"); // Name of mask
+        pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile;
+        if (sources && strlen(sources) > 0) {
+            if (!defineFile(config, srcInputFile, "PSPHOT.STACK.SOURCES.CNV", sources, PM_FPA_FILE_CMF)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)", i, sources);
+                return false;
+            }
+        }
+
+        // generate an pmFPAimage for the output convolved image
+        // XXX output of these files should be optional
+        {
+            pmFPAfile *outputImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.STACK.OUTPUT.IMAGE");
+            if (!outputImage) {
+                psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.STACK.OUTPUT.IMAGE");
+                return false;
+            }
+            outputImage->save = true;
+            outsources->fileID = i;             // this is used to generate output names
+
+            pmFPAfile *outputMask = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.MASK");
+            if (!outputMask) {
+                psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.STACK.OUTPUT.MASK"));
+                return NULL;
+            }
+            if (outputMask->type != PM_FPA_FILE_MASK) {
+                psError(PS_ERR_IO, true, "PSPHOT.STACK.OUTPUT.MASK is not of type MASK");
+                return NULL;
+            }
+            outputMask->save = true;
+
+            pmFPAfile *outputVariance = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.VARIANCE");
+            if (!outputVariance) {
+                psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.STACK.OUTPUT.VARIANCE"));
+                return NULL;
+            }
+            if (outputVariance->type != PM_FPA_FILE_VARIANCE) {
+                psError(PS_ERR_IO, true, "PSPHOT.STACK.OUTPUT.VARIANCE is not of type VARIANCE");
+                return NULL;
+            }
+            outputVariance->save = true;
+
+            // the output sources are carried on the outputImage->fpa structures
+            pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, outputImage, "PSPHOT.STACK.OUTPUT");
+            if (!outsources) {
+                psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT");
+                return false;
+            }
+            outsources->save = true;
+            outsources->fileID = i;             // this is used to generate output names
+        }
     }
     psMetadataRemoveKey(config->arguments, "FILENAMES");
@@ -71 +156 @@
 
     // generate an pmFPAimage for the chisqImage
-    pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE");
-    if (!chisqImage) {
-        psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE");
-        return false;
-    }
-    chisqImage->save = true;
-
-    pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK");
-    if (!chisqMask) {
-        psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK"));
-        return NULL;
-    }
-    if (chisqMask->type != PM_FPA_FILE_MASK) {
-        psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK");
-        return NULL;
-    }
-    chisqMask->save = true;
-
-    pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE");
-    if (!chisqVariance) {
-        psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE"));
-        return NULL;
-    }
-    if (chisqVariance->type != PM_FPA_FILE_VARIANCE) {
-        psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE");
-        return NULL;
-    }
-    chisqVariance->save = true;
-
-# if (0)    
-    // define the additional input/output files associated with psphot
-    // XXX figure out which files are needed by psphotStack
-    if (false && !psphotDefineFiles (config, input)) {
-        psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files");
-        return false;
-    }
-# endif
+    // XXX output of these files should be optional
+    {
+        pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE");
+        if (!chisqImage) {
+            psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE");
+            return false;
+        }
+        chisqImage->save = true;
+
+        pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK");
+        if (!chisqMask) {
+            psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK"));
+            return NULL;
+        }
+        if (chisqMask->type != PM_FPA_FILE_MASK) {
+            psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK");
+            return NULL;
+        }
+        chisqMask->save = true;
+
+        pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE");
+        if (!chisqVariance) {
+            psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE"));
+            return NULL;
+        }
+        if (chisqVariance->type != PM_FPA_FILE_VARIANCE) {
+            psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE");
+            return NULL;
+        }
+        chisqVariance->save = true;
+    }
 
     psTrace("psphot", 1, "Done with psphotStackParseCamera...\n");
@@ -145 +224 @@
 }
 
+
+/***
+ *
+ *  psphotStack :
+
+ *    * inputs:
+ *      * unconvolved images
+ *      * raw convolved images
+ *      * psfs (unconvolved or convolved?)
+ *      * sources
+ 
+ * optionally convolve the unconvolved or the raw inputs
+ * optionally perform no convolutions
+ * optionally save the psf-matched images
+
+ */
+
+    
+# if (0)    
+    // define the additional input/output files associated with psphot
+    // XXX figure out which files are needed by psphotStack
+    if (false && !psphotDefineFiles (config, input)) {
+        psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files");
+        return false;
+    }
+# endif
+