Changeset 16605

Timestamp:

Feb 22, 2008, 9:28:00 AM (18 years ago)

Author:

Paul Price

Message:

Merging in branch development. ppStack now works with incremental reads. Needs some cleanup, but it works.

Location:

trunk/ppStack/src

Files:

• Makefile.am (modified) (1 diff)
• ppStack.h (modified) (2 diffs)
• ppStackArguments.c (modified) (2 diffs)
• ppStackCamera.c (modified) (3 diffs)
• ppStackLoop.c (modified) (2 diffs)
• ppStackMatch.c (modified) (5 diffs)
• ppStackPSF.c (added)
• ppStackPhotometry.c (added)
• ppStackReadout.c (modified) (12 diffs)

trunk/ppStack/src/Makefile.am

@@ -9 +9 @@
         ppStackCamera.c         \
         ppStackLoop.c           \
+        ppStackPSF.c            \
         ppStackReadout.c        \
+        ppStackPhotometry.c     \
         ppStackVersion.c        \
         ppStackMatch.c

trunk/ppStack/src/ppStack.h

@@ -19 +19 @@
     );
 
+// Determine target PSF for input images
+pmPSF *ppStackPSF(const pmConfig *config, // Configuration
+                  int numCols, int numRows, // Size of image
+                  const psList *list    // List of input PSFs
+    );
+
 // Perform stacking on a readout
-bool ppStackReadout(pmConfig *config,   // Configuration
-                    const pmFPAview *view // View for readout
+bool ppStackReadout(const pmConfig *config,   // Configuration
+                    pmReadout *outRO,   // Output readout
+                    const psArray *readouts, // Input readouts
+                    const psArray *regions, // Array with array of regions used in each PSF matching
+                    const psArray *kernels // Array with array of kernels used in each PSF matching
+    );
+
+// Perform photometry on stack
+bool ppStackPhotometry(pmConfig *config, // Configuration
+                       const pmReadout *readout, // Readout to be photometered
+                       const pmFPAview *view // View to readout
     );
 
@@ -35 +50 @@
 
 /// Convolve image to match specified seeing
-bool ppStackMatch(pmReadout *output,    ///< Convolved readout
-                  const pmReadout *input, ///< Readout to be convolved
+bool ppStackMatch(pmReadout *readout, ///< Readout to be convolved; replaced with output
+                  psArray **regions, // Array of regions used in each PSF matching, returned
+                  psArray **kernels, // Array of kernels used in each PSF matching, returned
                   const pmReadout *sourcesRO, ///< Readout with sources
                   const pmPSF *psf,     ///< Target PSF

trunk/ppStack/src/ppStackArguments.c

@@ -116 +116 @@
     psMetadataAddU8(arguments,  PS_LIST_TAIL, "-mask-blank", 0, "Mask value for blank region", 0);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-rows", 0, "Rows to read at once", 128);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Do photometry on stacked image?", false);
     psMetadataAddS32(arguments, PS_LIST_TAIL, "-psf-instances", 0, "Number of instances for PSF generation", 5);
@@ -160 +161 @@
     VALUE_ARG_RECIPE_MASK("-mask-blank",      "MASK.BLANK");
     VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32);
+    VALUE_ARG_RECIPE_INT("-rows",             "ROWS",           S32, 0);
 
     VALUE_ARG_RECIPE_INT("-psf-instances", "PSF.INSTANCES", S32, 0);

trunk/ppStack/src/ppStackCamera.c

@@ -10 +10 @@
 
 #include "ppStack.h"
+
+
+#if 0
+// Define an output convolved image file
+static pmFPAfile *defineOutputConvolved(const char *name, // FPA file name
+                                        pmFPA *fpa, // FPA to bind
+                                        const pmConfig *config, // Configuration
+                                        pmFPAfileType type // Expected type
+    )
+{
+    pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, name);
+    if (!file) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to define output convolved file %s", name);
+        return NULL;
+    }
+    if (file->type != PM_FPA_FILE_IMAGE) {
+        psError(PS_ERR_IO, true, "PPSTACK.OUTCONV is not of type %s", pmFPAfileStringFromType(type));
+        return NULL;
+    }
+
+    return file;
+}
+
+// Define an input convolved image file, using the output as a basis
+static pmFPAfile *defineInputConvolved(const char *inputName, // Input FPA file name
+                                       pmFPAfile *outFile, // Corresponding output FPA file
+                                       pmConfig *config, // Configuration
+                                       pmFPAfileType type // Expected type
+    )
+{
+    pmFPAview *view = pmFPAviewAlloc(0); // View into sky cells
+    view->chip = view->cell = view->readout = 0;
+
+    psString imageName = pmFPANameFromRule(outFile->filerule, outFile->fpa, view);
+    psArray *imageNames = psArrayAlloc(1);
+    imageNames->data[0] = imageName;
+    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "INCONV.FILENAMES", PS_META_REPLACE,
+                       "Filenames of input convolved image files", imageNames);
+    psFree(imageNames);
+    bool found = false;                 // Found the file?
+    pmFPAfile *imageFile = pmFPAfileDefineFromArgs(&found, config, "PPSTACK.INCONV",
+                                                   "INCONV.FILENAMES");
+    psMetadataRemoveKey(config->arguments, "INCONV.FILENAMES");
+    if (!imageFile || !found) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to define %s file", inputName);
+        return NULL;
+    }
+    if (imageFile->type != type) {
+        psError(PS_ERR_IO, true, "PPSTACK.INCONV is not of type %s",
+                pmFPAfileStringFromType(type));
+        return NULL;
+    }
+
+    return imageFile;
+}
+#endif
+
+
 
 bool ppStackCamera(pmConfig *config)
@@ -138 +196 @@
             }
         }
+
+#if 0
+        // Output convolved files
+        pmFPAfile *outconvImage  = defineOutputConvolved("PPSTACK.OUTCONV", imageFile->fpa, config,
+                                                         PM_FPA_FILE_IMAGE);
+        pmFPAfile *outconvMask   = defineOutputConvolved("PPSTACK.OUTCONV.MASK", imageFile->fpa, config,
+                                                         PM_FPA_FILE_MASK);
+        pmFPAfile *outconvWeight = defineOutputConvolved("PPSTACK.OUTCONV.WEIGHT", imageFile->fpa, config,
+                                                         PM_FPA_FILE_WEIGHT);
+        if (!outconvImage || !outconvMask || !outconvWeight) {
+            return false;
+        }
+
+        // Input convolved files
+        pmFPAfile *inconvImage  = defineInputConvolved("PPSTACK.INCONV", outconvImage, config,
+                                                       PM_FPA_FILE_IMAGE);
+        pmFPAfile *inconvMask   = defineInputConvolved("PPSTACK.INCONV.MASK", outconvMask, config,
+                                                       PM_FPA_FILE_MASK);
+        pmFPAfile *inconvWeight = defineInputConvolved("PPSTACK.INCONV.WEIGHT", outconvWeight, config,
+                                                       PM_FPA_FILE_WEIGHT);
+        if (!inconvImage || !inconvMask || !inconvWeight) {
+            return false;
+        }
+#endif
 
         psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
@@ -241 +323 @@
     }
 
+
     // Output PSF
     if (havePSFs) {

trunk/ppStack/src/ppStackLoop.c

@@ -11 +11 @@
 #include "ppStack.h"
 
+// Here follows lists of files for activation/deactivation at various stages.  Each must be NULL-terminated.
+
+#if 0
+// All files in the system
+static char *allFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT",
+                            "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
+                            "PSPHOT.PSF.SAVE", "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES",
+                            "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", "PSPHOT.BACKMDL.STDEV",
+                            "PSPHOT.BACKGND", "PSPHOT.BACKSUB", "SOURCE.PLOT.MOMENTS",
+                            "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", "PSPHOT.INPUT.CMF",
+                            0 };
+#endif
+
+// Files required in preparation for convolution
+static char *prepareFiles[] = { "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", 0 };
+
+// Files required for the convolution
+static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", 0 };
+
+// Output files for the combination
+static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", 0 };
+
+// Files for photometry
+static char *photFiles[] = { "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL",
+                             "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
+                             "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID",
+                             "PSPHOT.INPUT.CMF", 0 };
+
+//#define CONVOLVED_ALREADY               // Already have the convolution products --- testing
+
+
+// Activate/deactivate a list of files
+static void fileActivation(pmConfig *config, // Configuration
+                           char **files, // Files to turn on/off
+                           bool state   // Activation state
+    )
+{
+    assert(config);
+    assert(files);
+
+    for (int i = 0; files[i] != NULL; i++) {
+        pmFPAfileActivate(config->files, state, files[i]);
+    }
+    return;
+}
+
+// Activate/deactivate a single element for a list; return array of files
+static psArray *fileActivationSingle(pmConfig *config, // Configuration
+                                     char **files, // Files to turn on/off
+                                     bool state,   // Activation state
+                                     int num // Number of file in sequence
+                                     )
+{
+    assert(config);
+    assert(files);
+
+    psList *list = psListAlloc(NULL);   // List of files
+
+    for (int i = 0; files[i] != NULL; i++) {
+        pmFPAfile *file = pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file
+        psListAdd(list, PS_LIST_TAIL, file);
+    }
+
+    psArray *array = psListToArray(list);
+    psFree(list);
+
+    return array;
+}
+
+#if 0
+// Set the data level for a list of files
+static void fileSetDataLevel(pmConfig *config, // Configuration
+                             char **files, // Files for which to set level
+                             pmFPALevel level // Level to set
+                             )
+{
+    assert(config);
+    assert(files);
+
+    for (int i = 0; files[i] != NULL; i++) {
+        psArray *selected = pmFPAfileSelect(config->files, files[i]); // Selected files of interest
+        for (int j = 0; j < selected->n; j++) {
+            pmFPAfile *file = selected->data[j];
+            assert(file);
+            file->dataLevel = level;
+        }
+        psFree(selected);
+    }
+    return;
+}
+#endif
+
+// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
+static pmFPAview *filesIterateDown(pmConfig *config // Configuration
+                                  )
+{
+    assert(config);
+
+    pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    view->chip = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    view->cell = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    view->readout = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    return view;
+}
+
+// Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells
+static bool filesIterateUp(pmConfig *config // Configuration
+                           )
+{
+    assert(config);
+
+    pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
+    view->chip = view->cell = view->readout = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    view->readout = -1;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    view->cell = -1;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    view->chip = -1;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    return true;
+}
+
+#ifndef CONVOLVED_ALREADY
+// Write an image to a FITS file
+static bool writeImage(const char *name, // Name of image
+                       psMetadata *header, // Header
+                       const psImage *image // Image
+                       )
+{
+    assert(name);
+    assert(image);
+
+    psFits *fits = psFitsOpen(name, "w");
+    if (!fits) {
+        psError(PS_ERR_IO, false, "Unable to open FITS file to write image.");
+        return false;
+    }
+    if (!psFitsWriteImage(fits, header, image, 0, NULL)) {
+        psError(PS_ERR_IO, false, "Unable to write FITS image.");
+        return false;
+    }
+    psFitsClose(fits);
+    return true;
+}
+#endif
+
+
 bool ppStackLoop(pmConfig *config)
 {
+    assert(config);
+
     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
 
@@ -32 +202 @@
     }
 
-    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Token input file
-    if (!input) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find input data!\n");
-        return false;
-    }
     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file
     if (!output) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!\n");
-        return false;
-    }
-
-    pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
-    pmHDU *lastHDU = NULL;              // Last HDU that was updated
-
-    // Iterate over the FPA hierarchy
-    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-        return false;
-    }
-
-    pmChip *chip;                       // Chip of interest
-    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-            return false;
-        }
-
-        pmCell *cell;                // Cell of interest
-        while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!");
+        return false;
+    }
+    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    int numScans = psMetadataLookupS32(NULL, config->arguments, "ROWS"); // Number of scans to read at once
+    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    int overlap = 2 * psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans
+
+    // Preparation iteration: Load the sources, and get a target PSF model
+    pmReadout *sources = NULL;          // Readout with sources to use for PSF matching
+    pmPSF *targetPSF = NULL;            // Target PSF
+    {
+        pmFPAfileActivate(config->files, false, NULL);
+        fileActivation(config, prepareFiles, true);
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            return false;
+        }
+
+        // We want to hang on to the 'sources' even when its host FPA is blown away
+        sources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"));
+        if (!sources) {
+            psError(PS_ERR_UNKNOWN, true, "Unable to find sources.");
+            psFree(view);
+            return false;
+        }
+
+        // Generate list of PSFs
+        psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
+                                                               "^PPSTACK.INPUT$"); // Iterator
+        psMetadataItem *fileItem; // Item from iteration
+        psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope
+        int numCols = 0, numRows = 0;   // Size of image
+        while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
+            assert(fileItem->type == PS_DATA_UNKNOWN);
+            pmFPAfile *inputFile = fileItem->data.V; // An input file
+            pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
+            pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
+            if (!psf) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
+                psFree(view);
+                psFree(sources);
+                psFree(fileIter);
+                psFree(psfList);
                 return false;
             }
-
-            // Put version information into the header
+            psListAdd(psfList, PS_LIST_TAIL, psf);
+
+            pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
             pmHDU *hdu = pmHDUFromCell(cell);
-            if (hdu && hdu != lastHDU) {
-                if (!hdu->header) {
-                    hdu->header = psMetadataAlloc();
-                }
-                ppStackVersionMetadata(hdu->header);
-                lastHDU = hdu;
-            }
-
-            pmReadout *readout;         // Readout of interest
-            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1))) {
-                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            assert(hdu && hdu->header);
+            int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns
+            int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows
+            if (naxis1 <= 0 || naxis2 <= 0) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF.");
+                psFree(view);
+                psFree(sources);
+                psFree(fileIter);
+                psFree(psfList);
+                return false;
+            }
+            if (numCols == 0 && numRows == 0) {
+                numCols = naxis1;
+                numRows = naxis2;
+            }
+        }
+        psFree(fileIter);
+        psFree(view);
+
+        targetPSF = ppStackPSF(config, numCols, numRows, psfList);
+        psFree(psfList);
+        if (!targetPSF) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            psFree(sources);
+            return false;
+        }
+
+        filesIterateUp(config);
+    }
+
+    const char *suffix = "conv.fits";   // Suffix for convolved images; ultimately this will be from recipe
+    const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
+    assert(outName);
+    psArray *imageNames = psArrayAlloc(num);
+    psArray *maskNames = psArrayAlloc(num);
+    psArray *weightNames = psArrayAlloc(num);
+    for (int i = 0; i < num; i++) {
+        psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images
+        psStringAppend(&imageName, "%s.im-%d.%s", outName, i, suffix);
+        psStringAppend(&maskName, "%s.mk-%d.%s", outName, i, suffix);
+        psStringAppend(&weightName, "%s.wt-%d.%s", outName, i, suffix);
+        imageNames->data[i] = imageName;
+        maskNames->data[i] = maskName;
+        weightNames->data[i] = weightName;
+    }
+
+    // Generate convolutions and write them to disk
+    psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
+    psArray *subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking
+    psArray *subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking
+    for (int i = 0; i < num; i++) {
+        pmFPAfileActivate(config->files, false, NULL);
+        psArray *files = fileActivationSingle(config, convolveFiles, true, i);
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            psFree(sources);
+            psFree(targetPSF);
+            return false;
+        }
+
+        pmReadout *readout = pmFPAviewThisReadout(view, ((pmFPAfile*)files->data[0])->fpa); // Input readout
+        psFree(view);
+
+#ifndef CONVOLVED_ALREADY
+        // Background subtraction, scaling and normalisation is performed automatically by the image matching
+        psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
+        if (!ppStackMatch(readout, &regions, &kernels, sources, targetPSF, config)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to match image %d --- ignoring.", i);
+            psFree(sources);
+            psFree(targetPSF);
+            return false;
+        }
+
+        subRegions->data[i] = regions;
+        subKernels->data[i] = kernels;
+
+        // Write the temporary convolved files
+        pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
+        assert(hdu);
+        writeImage(imageNames->data[i],  hdu->header, readout->image);
+        writeImage(maskNames->data[i],   hdu->header, readout->mask);
+        writeImage(weightNames->data[i], hdu->header, readout->weight);
+#endif
+
+        cells->data[i] = psMemIncrRefCounter(readout->parent);
+        filesIterateUp(config);
+    }
+    psFree(sources);
+    psFree(targetPSF);
+
+    // Stack the convolved files
+    {
+        pmFPAfileActivate(config->files, false, NULL);
+        fileActivation(config, combineFiles, true);
+        if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY")) {
+            fileActivation(config, photFiles, true);
+        }
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            psFree(cells);
+            psFree(subKernels);
+            psFree(subRegions);
+            return false;
+        }
+        pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
+        pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
+
+        psArray *readouts = psArrayAlloc(num); // Readouts for convolved images
+        for (int i = 0; i < num; i++) {
+            readouts->data[i] = pmReadoutAlloc(cells->data[i]); // Readout into which to read
+        }
+        psFree(cells);
+
+        // FITS files
+        psArray *imageFits  = psArrayAlloc(num);
+        psArray *maskFits   = psArrayAlloc(num);
+        psArray *weightFits = psArrayAlloc(num);
+        for (int i = 0; i < num; i++) {
+            imageFits->data[i] = psFitsOpen(imageNames->data[i], "r");
+            maskFits->data[i] = psFitsOpen(maskNames->data[i], "r");
+            weightFits->data[i] = psFitsOpen(weightNames->data[i], "r");
+        }
+
+        // Read convolutions by chunks
+        bool more = true;               // More to read?
+        for (int numChunk = 0; more; numChunk++) {
+            for (int i = 0; i < num; i++) {
+                pmReadout *readout = readouts->data[i];
+                assert(readout);
+
+                if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, overlap) ||
+                    !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, overlap) ||
+                    !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, overlap)) {
+                    psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
+                    psFree(readouts);
+                    psFree(subKernels);
+                    psFree(subRegions);
+                    psFree(outRO);
+                    psFree(view);
                     return false;
                 }
-
-                // Perform the analysis
-                if (!ppStackReadout(config, view)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
-                    return false;
-                }
-
-                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-                    return false;
-                }
-            }
-
-            // Perform statistics on the cell
-            if (stats) {
-                ppStatsFPA(stats, output->fpa, view, maskBlank, config);
-            }
-
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            }
+
+#ifdef TESTING
+            {
+                pmReadout *ro = readouts->data[0];
+                psTrace("ppStack", 1, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols,
+                        ro->row0, ro->row0 + ro->image->numRows);
+            }
+#endif
+
+            if (!ppStackReadout(config, outRO, readouts, subRegions, subKernels)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
+                psFree(readouts);
+                psFree(subKernels);
+                psFree(subRegions);
+                psFree(outRO);
+                psFree(view);
                 return false;
             }
-        }
-
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-            return false;
-        }
-    }
-
-    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-        return false;
-    }
-
-    psFree(view);
+
+            for (int i = 0; i < num && more; i++) {
+                pmReadout *readout = readouts->data[i];
+                assert(readout);
+                more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans);
+                more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans);
+                more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans);
+            }
+        }
+
+        psFree(readouts);
+        psFree(subKernels);
+        psFree(subRegions);
+        for (int i = 0; i < num; i++) {
+            psFitsClose(imageFits->data[i]);
+            psFitsClose(maskFits->data[i]);
+            psFitsClose(weightFits->data[i]);
+        }
+
+        if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY") &&
+            !ppStackPhotometry(config, outRO, view)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
+            psFree(outRO);
+            psFree(view);
+            return false;
+        }
+
+        // Statistics on output
+        if (stats) {
+            ppStatsFPA(stats, outCell->parent->parent, view, maskBlank, config);
+        }
+
+        // Put version information into the header
+        pmHDU *hdu = pmHDUFromCell(outRO->parent);
+        if (!hdu) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
+            return false;
+        }
+        if (!hdu->header) {
+            hdu->header = psMetadataAlloc();
+        }
+        ppStackVersionMetadata(hdu->header);
+
+        // Write out the output files
+        fileActivation(config, combineFiles, true);
+        filesIterateUp(config);
+
+        psFree(outRO);
+        psFree(view);
+    }
+
 
     // Write out summary statistics

trunk/ppStack/src/ppStackMatch.c

@@ -9 +9 @@
 #include "ppStack.h"
 
+#define ARRAY_BUFFER 16                 // Number to add to array at a time
+
+
 //#define TESTING
 
-bool ppStackMatch(pmReadout *output, const pmReadout *input, const pmReadout *sourcesRO,
-                  const pmPSF *psf, const pmConfig *config)
+bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels,
+                  const pmReadout *sourcesRO, const pmPSF *psf, const pmConfig *config)
 {
+    assert(readout);
+    assert(regions && !*regions);
+    assert(kernels && !*kernels);
+    assert(sourcesRO);
+    assert(psf);
+    assert(config);
+
     // Look up appropriate values from the ppSub recipe
     bool mdok;                          // Status of MD lookup
@@ -68 +78 @@
     pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
 
-    if (!pmReadoutFakeFromSources(fake, input->image->numCols, input->image->numRows, sources, NULL, NULL,
+    if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL, NULL,
                                   psf, powf(10.0, -0.4 * maxMag), 0, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
@@ -85 +95 @@
 
     // Do the image matching
-    if (!pmSubtractionMatch(output, NULL, input, fake, footprint, regionSize, spacing, threshold,
+    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
+    if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold,
                             sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
                             binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
@@ -92 +103 @@
         psFree(fake);
         psFree(optWidths);
+        psFree(output);
         return false;
     }
@@ -97 +109 @@
     psFree(optWidths);
 
+    // Replace original images with convolved
+    psFree(readout->image);
+    psFree(readout->mask);
+    psFree(readout->weight);
+    readout->image  = psMemIncrRefCounter(output->image);
+    readout->mask   = psMemIncrRefCounter(output->mask);
+    readout->weight = psMemIncrRefCounter(output->weight);
+
+    // Extract the regions and solutions used in the image matching
+    // This stops them from being freed when we iterate back up the FPA
+    *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
+    {
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
+        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);
+    }
+    *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels
+    {
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_UNKNOWN);
+            // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
+            // images by scans.
+            pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
+            kernel->numCols = readout->image->numCols;
+            kernel->numRows = readout->image->numRows;
+
+            *kernels = psArrayAdd(*kernels, ARRAY_BUFFER, kernel);
+        }
+        psFree(iter);
+    }
+    assert((*regions)->n == (*kernels)->n);
+
+
+    psFree(output);
+
     return true;
 }

trunk/ppStack/src/ppStackReadout.c

@@ -10 +10 @@
 #include "ppStack.h"
 
-#define ARRAY_BUFFER 16                 // Number to add to array at a time
 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
 
-//#define NO_CONVOLUTION                  // Don't perform convolutions?
-//#define CONVOLUTION_FILES               // Write convolutions?
 //#define REJECTION_FILES                 // Write rejection mask?
 //#define INSPECTION_FILES                // Write inspection mask?
 
-bool ppStackReadout(pmConfig *config, const pmFPAview *view)
+static int sectionNum = 0;              // Section number; for debugging outputs
+
+
+bool ppStackReadout(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+                    const psArray *regions, const psArray *kernels)
 {
+    assert(config);
+    assert(outRO);
+    assert(readouts);
+    assert(regions);
+    assert(kernels);
+    assert(readouts->n == regions->n);
+    assert(regions->n == kernels->n);
+
+
     // Get the recipe values
-    bool mdok;                          // Status of MD lookup
     int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
     float combineRej = psMetadataLookupF32(NULL, config->arguments, "COMBINE.REJ"); // Combination threshold
@@ -27 +36 @@
     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
     float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution
-    bool photometry = psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY"); // Perform photometry?
-    int psfInstances = psMetadataLookupS32(&mdok, config->arguments, "PSF.INSTANCES"); // Number of instances for PSF
-    float psfRadius = psMetadataLookupF32(NULL, config->arguments, "PSF.RADIUS"); // Radius for PSF
-    const char *psfModel = psMetadataLookupStr(NULL, config->arguments, "PSF.MODEL"); // Model for PSF
-    int psfOrder = psMetadataLookupS32(&mdok, config->arguments, "PSF.ORDER"); // Spatial order for PSF
-
-    // Get the output target
-    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
-    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
-    pmFPA *outFPA = outCell->parent->parent; // Output FPA
-
-    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
-    assert(num > 0);
-
-    // Get the input sources
-    psArray *stack = psArrayAllocEmpty(ARRAY_BUFFER); // The stack of inputs
-    pmReadout *sources = pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"); // Sources for stamps
-    psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
-                                                           "^PPSTACK.INPUT$"); // Iterator over input files
-    psMetadataItem *fileItem;           // Item from iteration
-    int fileNum = 0;                    // Number of file
-    psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope
-    pmPSF *outPSF = NULL;               // Ouptut PSF
-    int numCols = 0, numRows = 0;       // Size of image
-    while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
-        assert(fileItem->type == PS_DATA_UNKNOWN);
-        pmFPAfile *inputFile = fileItem->data.V; // An input file
-        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
-        pmCell *cell = ro->parent;      // The cell
-        pmChip *chip = cell->parent;    // The chip: holds the PSF
-
-        bool mdok;                      // Status of MD lookup
-        pmPSF *psf = psMetadataLookupPtr(&mdok, chip->analysis, "PSPHOT.PSF");
-        if (mdok && psf) {
-            psListAdd(psfList, PS_LIST_TAIL, psf);
-            if (numCols == 0 && numRows == 0) {
-                numCols = ro->image->numCols;
-                numRows = ro->image->numRows;
-            }
-        }
-    }
-
-    bool convolve = false;              // Convolve the input images?
-    if (psfList->n > 0) {
-        psArray *psfArray = psListToArray(psfList); // Array of PSFs
-        outPSF = pmPSFEnvelope(numCols, numRows, psfArray, psfInstances, psfRadius, psfModel,
-                               psfOrder, psfOrder);
-        psFree(psfArray);
-        if (!outPSF) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
-            // XXX Free stuff
-            return false;
-        }
-        convolve = true;
-        PS_ASSERT_PTR_NON_NULL(sources, false);
-    }
-
-    // Iterate through again to get the convolved images (or not)
-
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
-    psRandom *rng = psRandomAlloc(0, PS_RANDOM_TAUS); // Random number generator
+
+    int num = readouts->n;              // Number of inputs
+    psArray *stack = psArrayAlloc(num); // Array for stacking
+
+    pmCell *outCell = outRO->parent;    // Output cell
+    pmChip *outChip = outCell->parent;  // Output chip
+    pmFPA *outFPA = outChip->parent;    // Output FPA
+
     float totExposure = 0.0;            // Total exposure time
     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
     psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
-
-    psMetadataIteratorSet(fileIter, PS_LIST_HEAD);
-    fileNum = 0;
-    while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
-        assert(fileItem->type == PS_DATA_UNKNOWN);
-        pmFPAfile *inputFile = fileItem->data.V; // An input file
-        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
+    for (int i = 0; i < num; i++) {
+        pmReadout *ro = readouts->data[i];
+        assert(ro);
+        pmFPA *fpa = ro->parent->parent->parent; // Parent FPA
 
         bool mdok;                      // Status of MD lookup
-        float weighting = psMetadataLookupF32(&mdok, inputFile->fpa->analysis,
-                                              "PPSTACK.WEIGHTING"); // Relative weighting
+        float weighting = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.WEIGHTING"); // Relative weight
         if (!mdok || !isfinite(weighting)) {
-            psWarning("No WEIGHTING supplied for image %d --- set to unity.", fileNum);
+            psWarning("No WEIGHTING supplied for image %d --- set to unity.", i);
             weighting = 1.0;
         }
-        float scale = psMetadataLookupF32(&mdok, inputFile->fpa->analysis, "PPSTACK.SCALE"); // Rel. scale
-        if (!mdok || !isfinite(scale)) {
-            psWarning("No SCALE supplied for image %d --- set to unity.", fileNum);
-            scale = 1.0;
-        }
 
         float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time
-#if 0
         if (!isfinite(exposure)) {
             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                     "CELL.EXPOSURE is not set for input file %ld", stack->n);
-            psFree(stats);
-            psFree(rng);
-            psFree(fileIter);
             psFree(fpaList);
             psFree(cellList);
-            psFree(outRO);
             psFree(stack);
             return false;
         }
-#endif
         totExposure += exposure;        // Total exposure time
 
-        // Generate convolved version of input
-        pmReadout *convolved;
-        if (convolve) {
-            convolved = pmReadoutAlloc(NULL); // Convolved version of input readout
-            // Background subtraction, scaling and normalisation is performed automatically by the image
-            // matching
-            if (!ppStackMatch(convolved, ro, sources, outPSF, config)) {
-                psWarning("Unable to match image %d --- ignoring.", fileNum);
-                psErrorClear();
-                psFree(convolved);
-                // XXX Free the bad image so it's not taking up good memory!
-                continue;
-            }
-
-#ifdef CONVOLUTION_FILES
-            if (convolved->image) {
-                psString name = NULL;           // Name of image
-                psStringAppend(&name, "convolved%03d_image.fits", fileNum);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, NULL, convolved->image, 0, NULL);
-                psFitsClose(fits);
-            }
-            if (convolved->mask) {
-                psString name = NULL;           // Name of image
-                psStringAppend(&name, "convolved%03d_mask.fits", fileNum);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, NULL, convolved->mask, 0, NULL);
-                psFitsClose(fits);
-            }
-            if (convolved->weight) {
-                psString name = NULL;           // Name of image
-                psStringAppend(&name, "convolved%03d_weight.fits", fileNum);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, NULL, convolved->weight, 0, NULL);
-                psFitsClose(fits);
-            }
-#endif // CONVOLUTION_FILES
-
-        } else {
-            // No convolution --- just use the unconvolved images as the "convolved" images, with some tweaks
-            convolved = psMemIncrRefCounter(ro);
-            sources = NULL;
-
-            // Brain-dead background subtraction
-            if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskBad, rng)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to get image background on image %ld", stack->n);
-                psFree(stats);
-                psFree(rng);
-                psFree(fileIter);
-                psFree(fpaList);
-                psFree(cellList);
-                psFree(stack);
-                psFree(outRO);
-                psFree(convolved);
-                return false;
-            }
-            psTrace("ppStack", 3, "Background for image %d is %f\n", fileNum, stats->robustMedian);
-            (void)psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(stats->robustMedian, PS_TYPE_F32));
-
-            // Apply scaling
-            (void)psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(1.0 / scale, PS_TYPE_F32));
-
-            // Normalise each input by the exposure time
-            float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE");// Exposure time
-            if (!isfinite(exposure)) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                        "CELL.EXPOSURE is not set for input file %ld", stack->n);
-                psFree(stats);
-                psFree(rng);
-                psFree(fileIter);
-                psFree(fpaList);
-                psFree(cellList);
-                psFree(outRO);
-                psFree(convolved);
-                psFree(stack);
-                return false;
-            }
-
-            (void)psBinaryOp(ro->image, ro->image, "/", psScalarAlloc(exposure, PS_TYPE_F32));
-            if (ro->weight) {
-                (void)psBinaryOp(ro->weight, ro->weight, "/", psScalarAlloc(exposure, PS_TYPE_F32));
-            }
-        }
-
-        if (fileNum == 0) {
+#if 0
+        if (i == 0) {
             // Copy astrometry over
-            pmFPA *fpa = ro->parent->parent->parent; // Template FPA
             pmHDU *hdu = fpa->hdu; // Template HDU
             pmHDU *outHDU = outFPA->hdu; // Output HDU
@@ -224 +78 @@
                 psWarning("Unable to find HDU at FPA level to copy astrometry.");
             } else {
-                if (!pmAstromReadWCS(outFPA, outCell->parent, hdu->header, 1.0)) {
+                if (!pmAstromReadWCS(outFPA, outChip, hdu->header, 1.0)) {
                     psErrorClear();
                     psWarning("Unable to read WCS astrometry from input FPA.");
@@ -231 +85 @@
                         outHDU->header = psMetadataAlloc();
                     }
-                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
+                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
                         psErrorClear();
                         psWarning("Unable to write WCS astrometry to output FPA.");
@@ -238 +92 @@
             }
         }
-
-        // Don't need the original any more!
-        psFree(ro);
+#endif
 
         // Ensure there is a mask, or pmStackCombine will complain
-        if (!convolved->mask) {
-            convolved->mask = psImageAlloc(convolved->image->numCols, convolved->image->numRows,
-                                           PS_TYPE_MASK);
-            psImageInit(convolved->mask, 0);
-        }
-
-        psListAdd(fpaList, PS_LIST_TAIL, inputFile->fpa);
+        if (!ro->mask) {
+            ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_MASK);
+            psImageInit(ro->mask, 0);
+        }
+
+        psListAdd(fpaList, PS_LIST_TAIL, fpa);
         psListAdd(cellList, PS_LIST_TAIL, ro->parent);
 
-        pmStackData *data = pmStackDataAlloc(convolved, weighting); // Data to stack
-        psFree(convolved);
-        psArrayAdd(stack, ARRAY_BUFFER, data);
-        psFree(data);                   // Drop reference
-
-        fileNum++;
-    }
-    psFree(fileIter);
-    psFree(stats);
-    psFree(rng);
-
-    if (fileNum == 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not enough good files to combine.");
-        psFree(fpaList);
-        psFree(cellList);
-        psFree(stack);
-        psFree(outRO);
-        return false;
+        stack->data[i] = pmStackDataAlloc(ro, weighting);
     }
 
@@ -277 +111 @@
         psFree(cellList);
         psFree(stack);
-        psFree(outRO);
         return false;
     }
@@ -283 +116 @@
 #ifdef INSPECTION_FILES
     {
-        psFits *fits = psFitsOpen("combined.fits", "w");
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "combined_%03d.fits", sectionNum);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
         psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
         psFitsClose(fits);
@@ -291 +127 @@
         pmStackData *data = stack->data[i]; // Data for this image
         psImage *inspected = psPixelsToMask(NULL, data->pixels,
-                                            psRegionSet(0, data->readout->image->numCols,
-                                                        0, data->readout->image->numRows),
+                                            psRegionSet(0, data->readout->image->numCols - 1,
+                                                        0, data->readout->image->numRows - 1),
                                             maskBlank);
         psString name = NULL;           // Name of image
-        psStringAppend(&name, "inspect%03d.fits", i);
+        psStringAppend(&name, "inspect_%03d_%03d.fits", sectionNum, i);
         psFits *fits = psFitsOpen(name, "w");
         psFree(name);
@@ -304 +140 @@
 #endif
 
-    for (int i = 0; i < stack->n; i++) {
+    // Reject pixels
+    for (int i = 0; i < num; i++) {
         pmStackData *data = stack->data[i]; // Data for this image
-        pmReadout *readout = data->readout; // Readout for this image
-
-        // Extract the regions and solutions used in the image matching
-        psArray *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
-        {
-            psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD,
-                                                               "^SUBTRACTION.REGION$"); // Iterator
-            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 *solutions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of solutions
-        {
-            psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD,
-                                                               "^SUBTRACTION.SOLUTION$"); // Iterator
-            psMetadataItem *item = NULL;// Item from iteration
-            while ((item = psMetadataGetAndIncrement(iter))) {
-                assert(item->type == PS_DATA_VECTOR);
-                solutions = psArrayAdd(solutions, ARRAY_BUFFER, item->data.V);
-            }
-            psFree(iter);
-        }
-        assert(regions->n == solutions->n);
-
-        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
-                                                            "SUBTRACTION.KERNEL"); // Kernels
-
-        psPixels *reject = pmStackReject(data->pixels, threshold, regions,
-                                         solutions, kernels); // List of pixels to reject
+        pmReadout *readout = data->readout; // Readout of interest
+        int col0 = readout->col0, row0 = readout->row0; // Offset for readout
+        int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
+
+        psRegion *valid = psRegionAlloc(col0, col0 + numCols, row0, row0 + numRows); // Valid region for rej
+        psPixels *reject = pmStackReject(data->pixels, valid, threshold, regions->data[i],
+                                         kernels->data[i]); // Pixels to reject
+        psFree(valid);
         psFree(data->pixels);
         data->pixels = reject;
-
-        psFree(solutions);
-        psFree(regions);
     }
 
@@ -352 +162 @@
         }
         psImage *rejected = psPixelsToMask(NULL, data->pixels,
-                                           psRegionSet(0, data->readout->image->numCols,
-                                                       0, data->readout->image->numRows),
+                                           psRegionSet(0, data->readout->image->numCols - 1,
+                                                       0, data->readout->image->numRows - 1),
                                            maskBlank);
         psString name = NULL;           // Name of image
-        psStringAppend(&name, "reject%03d.fits", i);
+        psStringAppend(&name, "reject_%03d_%03d.fits", sectionNum, i);
         psFits *fits = psFitsOpen(name, "w");
         psFree(name);
@@ -370 +180 @@
         psFree(cellList);
         psFree(stack);
-        psFree(outRO);
         return false;
-    }
-
-    if (!convolve) {
-        // Restore image to counts using the total exposure time
-        (void)psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
-        if (outRO->weight) {
-            (void)psBinaryOp(outRO->weight, outRO->weight, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
-        }
     }
 
@@ -397 +198 @@
     psFree(cellList);
 
-    if (photometry) {
-        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-        pmFPACopy(photFile->fpa, outRO->parent->parent->parent);
-
-        if (!psphotReadout(config, view)) {
-            psWarning("Unable to perform photometry on stacked image.");
-            psErrorStackPrint(stderr, "Error stack from photometry:");
-            psErrorClear();
-        }
-
-        pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
-    }
-
     psFree(stack);
-    psFree(outRO);                      // Drop reference
+
+    sectionNum++;
 
     return true;