Changeset 23341

Timestamp:

Mar 17, 2009, 10:00:02 AM (17 years ago)

Author:

Paul Price

Message:

Reorganising ppStackLoop: splitting into multiple separate functions which are joined using a common data carrier, ppStackOptions (which aren't really 'options', but that's the word I chose...).

Location:

trunk/ppStack/src

Files:

• Makefile.am (modified) (1 diff)
• ppStack.c (modified) (1 diff)
• ppStack.h (modified) (5 diffs)
• ppStackCleanup.c (added)
• ppStackCombineFinal.c (added)
• ppStackCombineInitial.c (added)
• ppStackConvolve.c (added)
• ppStackFiles.c (added)
• ppStackFinish.c (added)
• ppStackLoop.c (modified) (2 diffs)
• ppStackLoop.h (added)
• ppStackOptions.c (added)
• ppStackOptions.h (added)
• ppStackPhotometry.c (modified) (5 diffs)
• ppStackPrepare.c (added)
• ppStackReadout.c (modified) (3 diffs)
• ppStackReject.c (added)
• ppStackSetup.c (added)
• ppStackThread.c (modified) (4 diffs)
• ppStackThread.h (added)

trunk/ppStack/src/Makefile.am

@@ -17 +17 @@
         ppStackArguments.c      \
         ppStackCamera.c         \
+        ppStackFiles.c          \
         ppStackLoop.c           \
         ppStackPSF.c            \
         ppStackReadout.c        \
-        ppStackPhotometry.c     \
         ppStackVersion.c        \
         ppStackMatch.c          \
         ppStackSources.c        \
-        ppStackThread.c
+        ppStackThread.c         \
+        ppStackOptions.c        \
+        ppStackSetup.c          \
+        ppStackPrepare.c        \
+        ppStackConvolve.c       \
+        ppStackCombineInitial.c \
+        ppStackReject.c         \
+        ppStackCombineFinal.c   \
+        ppStackCleanup.c        \
+        ppStackPhotometry.c     \
+        ppStackFinish.c
 
 noinst_HEADERS =                \
-        ppStack.h
+        ppStack.h               \
+        ppStackLoop.h           \
+        ppStackOptions.h        \
+        ppStackThread.h
 
 CLEANFILES = *~

trunk/ppStack/src/ppStack.c

@@ -9 +9 @@
 
 #include "ppStack.h"
+#include "ppStackLoop.h"
 
 #define TIMER_NAME "PPSTACK"            // Name of timer

trunk/ppStack/src/ppStack.h

@@ -1 @@
-#ifndef PP_STACK_H
-#define PP_STACK_H
+#ifndef PPSTACK_H
+#define PPSTACK_H
 
 #define PPSTACK_RECIPE "PPSTACK"        // Name of the recipe
@@ -18 +18 @@
 } ppStackMask;
 
-// Thread for stacking chunks
-//
-// Each input file contributes a readout, into which is read a chunk from that file
-typedef struct {
-    psArray *readouts;                  // Input readouts to read and stack
-    bool read;                          // Has the scan been read?
-    bool busy;                          // Is the scan being processed?
-    int firstScan;                      // First row of the chunk to be read for this group
-    int lastScan;                       // Last row of the chunk to be read for this group
-} ppStackThread;
+// List of files
+typedef enum {
+    PPSTACK_FILES_PREPARE,              // Files for preparation
+    PPSTACK_FILES_CONVOLVE,             // Files for convolution
+    PPSTACK_FILES_COMBINE,              // Files for combination
+    PPSTACK_FILES_PHOT                  // Files for photometry
+} ppStackFileList;
 
-// Allocator
-ppStackThread *ppStackThreadAlloc(psArray *readouts // Inputs readouts to read and stack
-    );
-
-// Data for threads
-typedef struct {
-    psArray *threads;                   // Threads doing stacking
-    int lastScan;                       // Last row that's been read
-    psArray *imageFits;                 // FITS file pointers for images
-    psArray *maskFits;                  // FITS file pointers for masks
-    psArray *varianceFits;                // FITS file pointers for variances
-} ppStackThreadData;
-
-// Set up thread data
-ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, // Array of input cells
-                                          const psArray *imageNames, // Names of images to read
-                                          const psArray *maskNames, // Names of masks to read
-                                          const psArray *varianceNames, // Names of variance maps to read
-                                          const psArray *covariances, // Covariance matrices (already read)
-                                          const pmConfig *config // Configuration
-    );
-
-// Read chunk into the first available file thread
-ppStackThread *ppStackThreadRead(bool *status, // Status of read
-                                 ppStackThreadData *stack, // Stacks available for reading
-                                 pmConfig *config, // Configuration
-                                 int numChunk, // Chunk number (only for interest)
-                                 int overlap // Overlap between subsequent scans
-    );
-
-// Initialise the threads
-void ppStackThreadInit(void);
 
 // Setup command-line arguments
@@ -73 +38 @@
 // Parse cameras
 bool ppStackCamera(pmConfig *config     // Configuration
-    );
-
-// Loop over the inputs, doing the combination
-bool ppStackLoop(pmConfig *config       // Configuration
     );
 
@@ -127 +88 @@
     );
 
-// Perform photometry on stack
-bool ppStackPhotometry(pmConfig *config, // Configuration
-                       const pmReadout *readout, // Readout to be photometered
-                       const pmFPAview *view // View to readout
-    );
-
 // Return software version
 psString ppStackVersion(void);
@@ -171 +126 @@
     );
 
+/// Dump memory debugging information
+void ppStackMemDump(
+    const char *name                    ///< Stage name, for inclusion in the output file name
+    );
+
+/// Activate/deactivate a list of files
+void ppStackFileActivation(
+    pmConfig *config,                   // Configuration
+    ppStackFileList list,               // Files to turn on/off
+    bool state                          // Activation state
+    );
+
+// Activate/deactivate a single element for a list
+void ppStackFileActivationSingle(
+    pmConfig *config,                   // Configuration
+    ppStackFileList list,               // Files to turn on/off
+    bool state,                         // Activation state
+    int num                             // Number of file in sequence
+    );
+
+/// Iterate down the hierarchy, loading files
+///
+/// We can get away with this simplistic treatment of the FPA hierarchy because we're working on skycells.
+pmFPAview *ppStackFilesIterateDown(
+    pmConfig *config                    // Configuration
+    );
+
+/// Iterate up the hierarchy, writing files
+///
+/// We can get away with this simplistic treatment of the FPA hierarchy because we're working on skycells.
+bool ppStackFilesIterateUp(
+    pmConfig *config                    // Configuration
+    );
+
+/// Write an image to a FITS file
+bool ppStackWriteImage(
+    const char *name,                   // Name of image
+    psMetadata *header,                 // Header
+    const psImage *image,               // Image
+    pmConfig *config                    // Configuration
+    );
+
+
 #endif

trunk/ppStack/src/ppStackLoop.c

@@ -4 +4 @@
 
 #include <stdio.h>
-#include <unistd.h>
-#include <string.h>
-#include <libgen.h>
 #include <pslib.h>
 #include <psmodules.h>
-#include <ppStats.h>
 
 #include "ppStack.h"
-
-//#define TESTING
-
-#define WCS_TOLERANCE 0.001             // Tolerance for WCS
-#define PIXELS_BUFFER 1024              // Initial size of pixel lists
-
-// Here follows lists of files for activation/deactivation at various stages.  Each must be NULL-terminated.
-
-// Files required in preparation for convolution
-static char *prepareFiles[] = { "PPSTACK.INPUT.PSF", "PPSTACK.INPUT.SOURCES", "PPSTACK.TARGET.PSF", NULL };
-
-// Files required for the convolution
-static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL };
-
-// Output files for the combination
-static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE",
-                                "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL };
-
-// Files for photometry
-static char *photFiles[] = { "PSPHOT.INPUT", "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", NULL };
-
-static void memDump(const char *name)
-{
-    return;
-
-    static int num = 0;                 // Counter, to make files unique and give an idea of sequence
-
-    psString filename = NULL;           // Name of file
-    psStringAppend(&filename, "memdump_%s_%03d.txt", name, num);
-    FILE *memFile = fopen(filename, "w");
-    psFree(filename);
-
-    psMemBlock **leaks = NULL;
-    int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true);
-    fprintf(memFile, "# MemBlock Size Source\n");
-    unsigned long total = 0;            // Total memory used
-    for (int i = 0; i < numLeaks; i++) {
-        psMemBlock *mb = leaks[i];
-        fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize,
-                mb->file, mb->lineno);
-        total += mb->userMemorySize;
-    }
-    fclose(memFile);
-    psFree(leaks);
-
-    fprintf(stderr, "Memdump %s %d: Memory use: %ld, sbrk: %p\n", name, num, total, sbrk(0));
-    num++;
-}
-
-
-
-// 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
-static void 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);
-
-    for (int i = 0; files[i] != NULL; i++) {
-        pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file
-    }
-    return;
-}
-
-// 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;
-    }
-    psFree(view);
-    return true;
-}
-
-// Write an image to a FITS file
-static bool writeImage(const char *name, // Name of image
-                       psMetadata *header, // Header
-                       const psImage *image, // Image
-                       pmConfig *config // Configuration
-                       )
-{
-    assert(name);
-    assert(image);
-
-    psString resolved = pmConfigConvertFilename(name, config, true, true); // Resolved file name
-    psFits *fits = psFitsOpen(resolved, "w");
-    if (!fits) {
-        psError(PS_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved);
-        psFree(resolved);
-        return false;
-    }
-    if (!psFitsWriteImage(fits, header, image, 0, NULL)) {
-        psError(PS_ERR_IO, false, "Unable to write FITS image %s.", resolved);
-        psFitsClose(fits);
-        psFree(resolved);
-        return false;
-    }
-    psFitsClose(fits);
-    psFree(resolved);
-    return true;
-}
-
+#include "ppStackLoop.h"
 
 bool ppStackLoop(pmConfig *config)
@@ -183 +15 @@
 
     psTimerStart("PPSTACK_TOTAL");
+    psTimerStart("PPSTACK_STEPS");
 
-    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
-    psAssert(recipe, "We've thrown an error on this before.");
+    ppStackOptions *options = ppStackOptionsAlloc(); // Options for stacking
 
-    bool mdok;                          // Status of MD lookup
-    bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files?
-    const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
-    if (!tempDir) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe");
+    // Setup
+    psTrace("ppStack", 1, "Setup....\n");
+    if (!ppStackSetup(options, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to setup.");
+        psFree(options);
+        return false;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 0: Setup: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("setup");
+
+
+    // Preparation for stacking
+    psTrace("ppStack", 1, "Preparation for stacking: merging sources, determining target PSF....\n");
+    if (!ppStackPrepare(options, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to prepare for stacking.");
+        psFree(options);
+        return false;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 1: Load Sources and Generate Target PSF: %f sec",
+             psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("prepare");
+
+
+    // Convolve inputs
+    psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
+    if (!ppStackConvolve(options, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
+        psFree(options);
+        return false;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec",
+             psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("convolve");
+
+
+    // Start threading
+    ppStackThreadInit();
+    ppStackThreadData *stack = ppStackThreadDataSetup(options, config);
+    if (!stack) {
+        psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
+        psFree(options);
         return false;
     }
 
-    psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments, "OUTPUT")); // Name for temporary files
-    const char *tempName = basename(outputName);
-    if (!tempName) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files.");
+    // Initial combination
+    psTrace("ppStack", 1, "Initial stack of convolved images....\n");
+    if (!ppStackCombineInitial(stack, options, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to perform initial combination.");
+        psFree(stack);
+        psFree(options);
         return false;
     }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec",
+             psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("convolve");
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 4: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
 
-    const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images
-    const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks
-    const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for temp variance maps
-    if (!tempImage || !tempMask || !tempVariance) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
+
+    // Pixel rejection
+    psTrace("ppStack", 1, "Reject pixels....\n");
+    if (!ppStackReject(options, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to reject pixels.");
+        psFree(stack);
+        psFree(options);
         return false;
     }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Pixel Rejection: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("reject");
 
-    float threshold = psMetadataLookupF32(NULL, recipe, "THRESHOLD.MASK"); // Threshold for mask deconvolution
-    float imageRej = psMetadataLookupF32(NULL, recipe, "IMAGE.REJ"); // Maximum fraction of image to reject
-                                                                     // before rejecting entire image
-    float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
-
-    const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics
-    psMetadata *stats = NULL;           // Container for statistics
-    FILE *statsFile = NULL;             // File stream for statistics
-    if (statsName && strlen(statsName) > 0) {
-        psString resolved = pmConfigConvertFilename(statsName, config, true, true); // Resolved filename
-        statsFile = fopen(resolved, "w");
-        if (!statsFile) {
-            psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
-            psFree(resolved);
-            return false;
-        } else {
-            stats = psMetadataAlloc();
-        }
-        psFree(resolved);
-    }
-
-    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file
-    if (!output) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!");
-        return false;
-    }
-    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
-
-    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
-    int overlap = 2 * psMetadataLookupS32(NULL, ppsub,
-                                          "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans
-
-    if (!pmConfigMaskSetBits(NULL, NULL, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
-        return false;
-    }
-
-    memDump("start");
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 0 : Initialization and Configuration : %f sec", psTimerClear("PPSTACK_STEPS"));
-
-    // Preparation iteration: Load the sources, and get a target PSF model
-    psTrace("ppStack", 1, "Determining target PSF....\n");
-    psArray *sourceLists = psArrayAlloc(num); // Individual lists of sources for matching
-    pmPSF *targetPSF = NULL;            // Target PSF
-    float sumExposure = NAN;            // Sum of exposure times
-    psVector *inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
-    psVectorInit(inputMask, 0);
-    if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
-        pmFPAfileActivate(config->files, false, NULL);
-        fileActivation(config, prepareFiles, true);
-        pmFPAview *view = filesIterateDown(config);
-        if (!view) {
-            return false;
-        }
-
-        // Generate list of PSFs
-        psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
-                                                               "^PPSTACK.INPUT$"); // Iterator
-        psMetadataItem *fileItem; // Item from iteration
-        psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
-        int numCols = 0, numRows = 0;   // Size of image
-        int index = 0;                  // Index for file
-        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(sourceLists);
-                psFree(fileIter);
-                psFree(psfs);
-                psFree(inputMask);
-                return false;
-            }
-            psfs->data[index] = psMemIncrRefCounter(psf);
-            psMetadataRemoveKey(chip->analysis, "PSPHOT.PSF");
-
-            pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
-            pmHDU *hdu = pmHDUFromCell(cell);
-            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(sourceLists);
-                psFree(fileIter);
-                psFree(psfs);
-                psFree(inputMask);
-                return false;
-            }
-            if (numCols == 0 && numRows == 0) {
-                numCols = naxis1;
-                numRows = naxis2;
-            }
-
-            pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
-            psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
-            if (!sources) {
-                psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
-                return NULL;
-            }
-            sourceLists->data[index] = psMemIncrRefCounter(sources);
-
-            // Re-do photometry if we don't trust the source lists
-            if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
-                psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
-                pmFPAfileActivate(config->files, false, NULL);
-                fileActivationSingle(config, convolveFiles, true, index);
-                pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
-                pmFPAview *view = filesIterateDown(config);
-                if (!view) {
-                    psFree(sourceLists);
-                    psFree(targetPSF);
-                    psFree(inputMask);
-                   return false;
-                }
-
-                pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
-
-                if (!ppStackInputPhotometry(ro, sources, config)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
-                    psFree(sourceLists);
-                    psFree(targetPSF);
-                    psFree(inputMask);
-                    return false;
-                }
-
-                psFree(view);
-                if (!filesIterateUp(config)) {
-                    psFree(sourceLists);
-                    psFree(targetPSF);
-                    psFree(inputMask);
-                    return false;
-                }
-                pmFPAfileActivate(config->files, false, NULL);
-                fileActivation(config, prepareFiles, true);
-            }
-
-            index++;
-        }
-        psFree(fileIter);
-
-        // Zero point calibration
-        sumExposure = ppStackSourcesTransparency(sourceLists, inputMask, view, config);
-        if (!isfinite(sumExposure) || sumExposure <= 0) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
-            psFree(sourceLists);
-            psFree(targetPSF);
-            psFree(inputMask);
-            return false;
-        }
-
-        // Generate target PSF
-        targetPSF = ppStackPSF(config, numCols, numRows, psfs, inputMask);
-        psFree(psfs);
-        if (!targetPSF) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
-            psFree(sourceLists);
-            psFree(view);
-            psFree(inputMask);
-            return false;
-        }
-        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
-                         "Target PSF for stack", targetPSF);
-
-        pmChip *outChip = pmFPAviewThisChip(view, output->fpa); // Output chip
-        psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
-                         "Target PSF", targetPSF);
-        outChip->data_exists = true;
-
-        psFree(view);
-        if (!filesIterateUp(config)) {
-            return false;
-        }
-    }
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 1 : Load Sources and Generate Target PSF : %f sec", psTimerClear("PPSTACK_STEPS"));
-    memDump("psf");
-
-    psArray *imageNames = psArrayAlloc(num);
-    psArray *maskNames = psArrayAlloc(num);
-    psArray *varianceNames = psArrayAlloc(num);
-    for (int i = 0; i < num; i++) {
-        psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
-        psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage);
-        psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask);
-        psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
-        psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
-        imageNames->data[i] = imageName;
-        maskNames->data[i] = maskName;
-        varianceNames->data[i] = varianceName;
-    }
-    // Free the outputName that we grabbed above to set up tempName.
-    psFree(outputName);
-
-    // Generate convolutions and write them to disk
-    psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
-    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
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
-    int numGood = 0;                    // Number of good frames
-    int numCols = 0, numRows = 0;       // Size of image
-    psVector *matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching
-    psVectorInit(matchChi2, NAN);
-    psVector *weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
-    psVectorInit(weightings, NAN);
-    psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
-    psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
-    psArray *covariances = psArrayAlloc(num); // Covariance matrices
-    for (int i = 0; i < num; i++) {
-        if (inputMask->data.U8[i]) {
-            continue;
-        }
-        psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num);
-        pmFPAfileActivate(config->files, false, NULL);
-        fileActivationSingle(config, convolveFiles, true, i);
-        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
-        pmFPAview *view = filesIterateDown(config);
-        if (!view) {
-            psFree(sourceLists);
-            psFree(targetPSF);
-            psFree(rng);
-            psFree(inputMask);
-            psFree(matchChi2);
-            psFree(fpaList);
-            psFree(cellList);
-            psFree(covariances);
-            return false;
-        }
-
-        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout
-        psFree(view);
-
-        if (numCols == 0 && numRows == 0) {
-            numCols = readout->image->numCols;
-            numRows = readout->image->numRows;
-        } else if (numCols != readout->image->numCols || numRows != readout->image->numRows) {
-            psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d",
-                    readout->image->numCols, readout->image->numRows, numCols, numRows);
-            psFree(sourceLists);
-            psFree(targetPSF);
-            psFree(rng);
-            psFree(inputMask);
-            psFree(matchChi2);
-            psFree(fpaList);
-            psFree(cellList);
-            psFree(covariances);
-            return false;
-        }
-
-        // Background subtraction, scaling and normalisation is performed automatically by the image matching
-        psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
-        psTimerStart("PPSTACK_MATCH");
-
-        if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i],
-                          sourceLists->data[i], targetPSF, rng, config)) {
-            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
-            inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_MATCH;
-            psErrorClear();
-            continue;
-        }
-        covariances->data[i] = psMemIncrRefCounter(readout->covariance);
-
-        if (stats) {
-            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
-                                                                PM_SUBTRACTION_ANALYSIS_KERNEL);// Conv kernel
-            psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", PS_META_DUPLICATE_OK,
-                             "Time to match PSF", psTimerMark("PPSTACK_MATCH"));
-            psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK,
-                             "Mean deviation for stamps", kernels->mean);
-            psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK,
-                             "RMS deviation for stamps", kernels->rms);
-            psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK,
-                             "Number of stamps", kernels->numStamps);
-            float deconv = psMetadataLookupF32(NULL, readout->analysis,
-                                               PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Deconvolution fraction
-            psMetadataAddF32(stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK,
-                             "Deconvolution fraction for kernel", deconv);
-            psMetadataAddF32(stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK,
-                             "Weighting for image", weightings->data.F32[i]);
-        }
-        psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH"));
-
-        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, config);
-        psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
-        pmConfigMaskWriteHeader(config, maskHeader);
-        writeImage(maskNames->data[i], maskHeader, readout->mask, config);
-        psFree(maskHeader);
-        psImageCovarianceTransfer(readout->variance, readout->covariance);
-        writeImage(varianceNames->data[i], hdu->header, readout->variance, config);
-#ifdef TESTING
-        {
-            psString name = NULL;
-            psStringAppend(&name, "covariance_%d.fits", i);
-            writeImage(name, hdu->header, readout->covariance->image, config);
-            pmStackVisualPlotTestImage(readout->covariance->image, name);
-            psFree(name);
-        }
-#endif
-
-        pmCell *inCell = readout->parent; // Input cell
-
-        psListAdd(cellList, PS_LIST_TAIL, inCell);
-        psListAdd(fpaList, PS_LIST_TAIL, inCell->parent->parent);
-
-        cells->data[i] = psMemIncrRefCounter(inCell);
-        if (!filesIterateUp(config)) {
-            psFree(fpaList);
-            psFree(cellList);
-            psFree(covariances);
-            return false;
-        }
-        numGood++;
-
-        memDump("match");
-    }
-    psFree(sourceLists);
-    psFree(targetPSF);
-    psFree(rng);
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 2 : Generate Convolutions and Save : %f sec", psTimerClear("PPSTACK_STEPS"));
-
-    if (numGood == 0) {
-        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
-        psFree(subKernels);
-        psFree(subRegions);
-        psFree(cells);
-        psFree(inputMask);
-        psFree(matchChi2);
-        psFree(fpaList);
-        psFree(cellList);
-        psFree(covariances);
-        return false;
-    }
-
-
-    // Update concepts for output
-    {
-        pmFPAview view;                 // View for output
-        view.chip = view.cell = view.readout = 0;
-        pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell
-        pmFPA *outFPA = outCell->parent->parent; // Output FPA
-        pmConceptsAverageFPAs(outFPA, fpaList);
-        pmConceptsAverageCells(outCell, cellList, NULL, NULL, true);
-        psFree(fpaList);
-        psFree(cellList);
-
-        // Update the value of a concept
-#define UPDATE_CONCEPT(SOURCE, NAME, VALUE) { \
-            psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
-            psAssert(item, "Concept should be present"); \
-            psAssert(item->type == PS_TYPE_F32, "Concept should be F32"); \
-            item->data.F32 = VALUE; \
-        }
-
-        UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",  sumExposure);
-        UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", sumExposure);
-        UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
-    }
-
-
-    // Reject images out-of-hand on the basis of their match chi^2
-    {
-        psVector *values = psVectorAllocEmpty(num, PS_TYPE_F32); // Values to sort
-        for (int i = 0; i < num; i++) {
-            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
-                continue;
-            }
-            values->data.F32[values->n++] = matchChi2->data.F32[i];
-        }
-        assert(values->n == numGood);
-        if (!psVectorSortInPlace(values)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to sort vector.");
-            psFree(subKernels);
-            psFree(subRegions);
-            psFree(cells);
-            psFree(inputMask);
-            psFree(matchChi2);
-            psFree(values);
-            psFree(covariances);
-            return false;
-        }
-        float median = numGood % 2 ? values->data.F32[numGood / 2] :
-            0.5 * (values->data.F32[numGood / 2 - 1] + values->data.F32[numGood / 2]);
-
-        float rms = 0.74 * (values->data.F32[numGood * 3 / 4] -
-                            values->data.F32[numGood / 4]); // Estimated RMS from interquartile range
-
-        psFree(values);
-
-        float rej = psMetadataLookupF32(NULL, recipe, "MATCH.REJ"); // Rejection threshold (stdevs)
-        if (isfinite(rej)) {
-            float thresh = median + rej * rms; // Threshold for rejection
-            psLogMsg("ppStack", PS_LOG_INFO, "chi^2 rejection threshold = %f + %f * %f = %f",
-                     median, rej, rms, thresh);
-
-            int numRej = 0;                 // Number rejected
-            numGood = 0;                    // Number of good images
-            for (int i = 0; i < num; i++) {
-              if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
-                    continue;
-                }
-                if (matchChi2->data.F32[i] > thresh) {
-                    numRej++;
-                    inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CHI2;
-                    psLogMsg("ppStack", PS_LOG_INFO, "Rejecting image %d because of large matching chi^2: %f",
-                             i, matchChi2->data.F32[i]);
-                } else {
-                    psLogMsg("ppStack", PS_LOG_INFO, "Image %d has matching chi^2: %f",
-                             i, matchChi2->data.F32[i]);
-                    numGood++;
-                }
-            }
-        }
-    }
-
-    if (numGood == 0) {
-        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
-        psFree(subKernels);
-        psFree(subRegions);
-        psFree(cells);
-        psFree(inputMask);
-        psFree(matchChi2);
-        psFree(covariances);
-        return false;
-    }
-
-#ifdef TESTING
-    psTraceSetLevel("psModules.imcombine", 7);
-#endif
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 3 : Basic Image Rejection  : %f sec", psTimerClear("PPSTACK_STEPS"));
-
-    // Start threading
-    ppStackThreadInit();
-    ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, varianceNames,
-                                                      covariances, config);
-    if (!stack) {
-        psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
-        psFree(subKernels);
-        psFree(subRegions);
-        psFree(inputMask);
-        psFree(matchChi2);
-        psFree(cells);
-        psFree(covariances);
-        return false;
-    }
-
-    psTimerStart("PPSTACK_INITIAL");
-
-    memDump("preinitial");
-
-    // Stack the convolved files
-    psTrace("ppStack", 1, "Initial stack of convolved images....\n");
-    pmReadout *outRO = NULL;            // Output readout
-    pmFPAview *view = NULL;             // View to readout
-    psArray *inspect = NULL;            // Array of arrays of pixels to inspect
-    {
-        int row0, col0;                 // Offset for readout
-        int numCols, numRows;           // Size of readout
-        ppStackThread *thread = stack->threads->data[0]; // Representative thread
-        if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) {
-            psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
-            psFree(subKernels);
-            psFree(subRegions);
-            psFree(stack);
-            psFree(inputMask);
-            psFree(matchChi2);
-            psFree(cells);
-            psFree(covariances);
-            return false;
-        }
-
-        pmFPAfileActivate(config->files, false, NULL);
-        fileActivation(config, combineFiles, true);
-        view = filesIterateDown(config);
-        if (!view) {
-            psFree(subKernels);
-            psFree(subRegions);
-            psFree(stack);
-            psFree(inputMask);
-            psFree(matchChi2);
-            psFree(cells);
-            psFree(covariances);
-            return false;
-        }
-
-        pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
-        outRO = pmReadoutAlloc(outCell); // Output readout
-
-        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
-        psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
-        if (!pmReadoutStackDefineOutput(outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to prepare output.");
-            psFree(subKernels);
-            psFree(subRegions);
-            psFree(stack);
-            psFree(inputMask);
-            psFree(matchChi2);
-            psFree(view);
-            psFree(outRO);
-            psFree(cells);
-            psFree(covariances);
-            return false;
-        }
-
-        psFree(cells);
-
-        bool status;                    // Status of read
-        int numChunk;                   // Number of chunks
-        for (numChunk = 0; true; numChunk++) {
-            ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, overlap);
-            if (!status) {
-                // Something went wrong
-                psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
-                psFree(subKernels);
-                psFree(subRegions);
-                psFree(stack);
-                psFree(inputMask);
-                psFree(matchChi2);
-                psFree(view);
-                psFree(outRO);
-                psFree(covariances);
-                return false;
-            }
-            if (!thread) {
-                // Nothing more to read
-                break;
-            }
-
-            // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)
-            psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start
-            psArrayAdd(job->args, 1, thread);
-            psArrayAdd(job->args, 1, config);
-            psArrayAdd(job->args, 1, outRO);
-            psArrayAdd(job->args, 1, inputMask);
-            psArrayAdd(job->args, 1, weightings);
-            psArrayAdd(job->args, 1, matchChi2);
-            if (!psThreadJobAddPending(job)) {
-                psFree(job);
-                psFree(subKernels);
-                psFree(subRegions);
-                psFree(stack);
-                psFree(inputMask);
-                psFree(weightings);
-                psFree(matchChi2);
-                psFree(view);
-                psFree(outRO);
-                psFree(covariances);
-                return false;
-            }
-            psFree(job);
-        }
-
-        if (!psThreadPoolWait(false)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination.");
-            psFree(subKernels);
-            psFree(subRegions);
-            psFree(stack);
-            psFree(inputMask);
-            psFree(matchChi2);
-            psFree(view);
-            psFree(outRO);
-            psFree(covariances);
-            return false;
-        }
-
-        // Harvest the jobs, gathering the inspection lists
-        inspect = psArrayAlloc(num);
-        for (int i = 0; i < num; i++) {
-            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-                continue;
-            }
-            inspect->data[i] = psArrayAllocEmpty(numChunk);
-        }
-        psThreadJob *job;               // Completed job
-        while ((job = psThreadJobGetDone())) {
-            psAssert(strcmp(job->type, "PPSTACK_INITIAL_COMBINE") == 0,
-                     "Job has incorrect type: %s", job->type);
-            psArray *results = job->results; // Results of job
-            for (int i = 0; i < num; i++) {
-                if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-                    continue;
-                }
-                inspect->data[i] = psArrayAdd(inspect->data[i], 1, results->data[i]);
-            }
-            psFree(job);
-        }
-
-        memDump("initial");
-    }
-
-#ifdef TESTING
-    writeImage("combined_initial.fits", NULL, outRO->image, config);
-    pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits");
-#endif
-
-    if (stats) {
-        psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_INITIAL", 0,
-                         "Time to make initial stack", psTimerMark("PPSTACK_INITIAL"));
-    }
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 4 : Make Initial Stack : %f sec", psTimerClear("PPSTACK_STEPS"));
-    psLogMsg("ppStack", PS_LOG_INFO, "Time to make initial stack: %f sec", psTimerClear("PPSTACK_INITIAL"));
-
-    psTimerStart("PPSTACK_REJECT");
-
-    // Pixel rejection
-    psArray *rejected = psArrayAlloc(num);
-    {
-        int numRejected = 0;        // Number of inputs rejected completely
-
-        // Count images rejected out of hand
-        for (int i = 0; i < num; i++) {
-            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-                numRejected++;
-            }
-        }
-
-        for (int i = 0; i < num; i++) {
-            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-                continue;
-            }
-
-            psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start
-            psArrayAdd(job->args, 1, inspect);
-            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
-            if (!psThreadJobAddPending(job)) {
-                psFree(job);
-                psFree(subKernels);
-                psFree(subRegions);
-                psFree(stack);
-                psFree(inputMask);
-                psFree(view);
-                psFree(outRO);
-                psFree(inspect);
-                psFree(rejected);
-                psFree(covariances);
-                psFree(matchChi2);
-                return false;
-            }
-            psFree(job);
-        }
-
-        if (!psThreadPoolWait(true)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to concatenate inspection lists.");
-            psFree(subKernels);
-            psFree(subRegions);
-            psFree(stack);
-            psFree(inputMask);
-            psFree(view);
-            psFree(outRO);
-            psFree(inspect);
-            psFree(rejected);
-            psFree(covariances);
-            psFree(matchChi2);
-            return false;
-        }
-
-        if (psMetadataLookupS32(NULL, config->arguments, "-threads") > 0) {
-            pmStackRejectThreadsInit();
-        }
-
-        int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches
-
-        // Reject bad pixels
-        for (int i = 0; i < num; i++) {
-            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-                continue;
-            }
-            psTimerStart("PPSTACK_REJECT");
-
-#ifdef TESTING
-            {
-                psImage *mask = psPixelsToMask(NULL, inspect->data[i],
-                                               psRegionSet(0, numCols - 1, 0, numRows - 1),
-                                               0xff); // Mask image
-                psString name = NULL;           // Name of image
-                psStringAppend(&name, "inspect_%03d.fits", i);
-                pmStackVisualPlotTestImage(mask, name);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, NULL, mask, 0, NULL);
-                psFree(mask);
-                psFitsClose(fits);
-            }
-#endif
-
-            psPixels *reject = pmStackReject(inspect->data[i], numCols, numRows, threshold, poorFrac, stride,
-                                             subRegions->data[i], subKernels->data[i]); // Rejected pixels
-
-#ifdef TESTING
-            {
-                psImage *mask = psPixelsToMask(NULL, reject, psRegionSet(0, numCols - 1, 0, numRows - 1),
-                                               0xff); // Mask image
-                psString name = NULL;           // Name of image
-                psStringAppend(&name, "reject_%03d.fits", i);
-                pmStackVisualPlotTestImage(mask, name);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, NULL, mask, 0, NULL);
-                psFree(mask);
-                psFitsClose(fits);
-            }
-#endif
-
-            psFree(inspect->data[i]);
-            inspect->data[i] = NULL;
-
-            if (!reject) {
-                psWarning("Rejection on image %d didn't work --- reject entire image.", i);
-                numRejected++;
-                inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_REJECT;
-            } else {
-                float frac = reject->n / (float)(numCols * numRows); // Pixel fraction
-                psLogMsg("ppStack", PS_LOG_INFO, "%ld pixels rejected from image %d (%.1f%%)",
-                         reject->n, i, frac * 100.0);
-                if (frac > imageRej) {
-                    psWarning("Image %d rejected completely because rejection fraction (%.3f) "
-                              "exceeds limit (%.3f)", i, frac, imageRej);
-                    psFree(reject);
-                    // reject == NULL means reject image completely
-                    reject = NULL;
-                    inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_BAD;
-                    numRejected++;
-                }
-            }
-
-            // Images without a list of rejected pixels (the list may be empty) are rejected completely
-            rejected->data[i] = reject;
-
-            if (stats) {
-                psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_REJECT", PS_META_DUPLICATE_OK,
-                                 "Time to perform rejection", psTimerMark("PPSTACK_REJECT"));
-                psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_PIXELS", PS_META_DUPLICATE_OK,
-                                 "Number of pixels rejected", reject ? reject->n : 0);
-            }
-            psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i,
-                     psTimerClear("PPSTACK_REJECT"));
-        }
-
-        psFree(inspect);
-        psFree(subKernels);
-        psFree(subRegions);
-
-        if (numRejected >= num - 1) {
-            psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed.");
-            psFree(stack);
-            psFree(rejected);
-            psFree(inputMask);
-            psFree(view);
-            psFree(outRO);
-            psFree(covariances);
-            psFree(matchChi2);
-            return false;
-        }
-
-        if (stats) {
-            psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_IMAGES", 0,
-                             "Number of images rejected completely", numRejected);
-        }
-    }
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 5 : Pixel Rejection : %f sec", psTimerClear("PPSTACK_STEPS"));
-    psTimerStart("PPSTACK_FINAL");
-
-    memDump("reject");
 
     // Final combination
     psTrace("ppStack", 2, "Final stack of convolved images....\n");
-    {
-        stack->lastScan = 0;            // Reset read
-        bool status;                    // Status of read
-        for (int numChunk = 0; true; numChunk++) {
-            ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0);
-            if (!status) {
-                // Something went wrong
-                psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
-                psFree(stack);
-                psFree(rejected);
-                psFree(inputMask);
-                psFree(view);
-                psFree(outRO);
-                psFree(covariances);
-                psFree(matchChi2);
-                return false;
-            }
-            if (!thread) {
-                // Nothing more to read
-                break;
-            }
+    if (!ppStackCombineFinal(stack, options, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination.");
+        psFree(stack);
+        psFree(options);
+        return false;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("final");
 
-            // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
-            psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
-            psArrayAdd(job->args, 1, thread);
-            psArrayAdd(job->args, 1, config);
-            psArrayAdd(job->args, 1, outRO);
-            psArrayAdd(job->args, 1, inputMask);
-            psArrayAdd(job->args, 1, rejected);
-            psArrayAdd(job->args, 1, weightings);
-            psArrayAdd(job->args, 1, matchChi2);
-            if (!psThreadJobAddPending(job)) {
-                psFree(job);
-                psFree(stack);
-                psFree(rejected);
-                psFree(inputMask);
-                psFree(view);
-                psFree(outRO);
-                psFree(covariances);
-                psFree(matchChi2);
-                return false;
-            }
-            psFree(job);
-        }
 
-        if (!psThreadPoolWait(true)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to do final combination.");
-            psFree(stack);
-            psFree(inputMask);
-            psFree(rejected);
-            psFree(view);
-            psFree(outRO);
-            psFree(covariances);
-            psFree(matchChi2);
-            return false;
-        }
+    // Clean up
+    psTrace("ppStack", 2, "Cleaning up after combination....\n");
+    if (!ppStackCombineFinal(stack, options, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to clean up.");
+        psFree(stack);
+        psFree(options);
+        return false;
     }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("cleanup");
 
-    memDump("final");
-
-#ifdef TESTING
-    pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits");
-    writeImage("combined_final.fits", NULL, outRO->image, config);
-#endif
-
-    // Sum covariance matrices
-    for (int i = 0; i < num; i++) {
-        if (inputMask->data.U8[i]) {
-            psFree(covariances->data[i]);
-            covariances->data[i] = NULL;
-        }
-    }
-    outRO->covariance = psImageCovarianceSum(covariances);
-    psFree(covariances);
-    psFree(matchChi2);
-
-    if (stats) {
-        psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_FINAL", 0,
-                         "Time to make final stack", psTimerMark("PPSTACK_FINAL"));
-    }
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 6 : Final Stack : %f sec", psTimerClear("PPSTACK_STEPS"));
-    psLogMsg("ppStack", PS_LOG_INFO, "Time to make final stack: %f sec", psTimerClear("PPSTACK_FINAL"));
-
-    psTrace("ppStack", 2, "Cleaning up after combination....\n");
-
-#if 0
-    // Ensure masked regions really look masked
-    {
-        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
-        psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
-        if (!pmReadoutMaskApply(outRO, maskBad)) {
-            psWarning("Unable to apply mask");
-        }
-    }
-#endif
-
-    // Close up
-    bool wcsDone = false;           // Have we done the WCS?
-    for (int i = 0; i < num; i++) {
-        if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-            continue;
-        }
-
-        ppStackThread *thread = stack->threads->data[0]; // Representative stack
-        pmReadout *inRO = thread->readouts->data[i]; // Template readout
-        if (inRO && !wcsDone) {
-            // Copy astrometry over
-            wcsDone = true;
-            pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
-            pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
-            pmChip *outChip = outRO->parent->parent; // Output chip
-            pmFPA *outFPA = outChip->parent; // Output FPA
-            if (!outHDU || !inHDU) {
-                psWarning("Unable to find HDU at FPA level to copy astrometry.");
-            } else {
-                if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
-                    psErrorClear();
-                    psWarning("Unable to read WCS astrometry from input FPA.");
-                    wcsDone = false;
-                } else {
-                    if (!outHDU->header) {
-                        outHDU->header = psMetadataAlloc();
-                    }
-                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
-                        psErrorClear();
-                        psWarning("Unable to write WCS astrometry to output FPA.");
-                        wcsDone = false;
-                    }
-                }
-            }
-        }
-
-        if (tempDelete) {
-            psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
-            psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
-            psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false);
-            if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
-                unlink(varianceResolved) == -1) {
-                psWarning("Unable to delete temporary files for image %d", i);
-            }
-            psFree(imageResolved);
-            psFree(maskResolved);
-            psFree(varianceResolved);
-        }
-    }
-    psFree(imageNames);
-    psFree(maskNames);
-    psFree(varianceNames);
-
-    psFree(inputMask);
     psFree(stack);
 
-    memDump("cleanup");
 
-    // Generate binned JPEGs
-    {
-        int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level
-        int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level
-
-        // Target cells
-        pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG1");
-        pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG2");
-        psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
-
-        pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
-        if (!pmReadoutRebin(ro1, outRO, maskValue, bin1, bin1) || !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to bin output.");
-            psFree(ro1);
-            psFree(ro2);
-            psFree(outRO);
-            return false;
-        }
-        psFree(ro1);
-        psFree(ro2);
-    }
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7 : WCS & JPEGS : %f sec", psTimerClear("PPSTACK_STEPS"));
-
-    if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
-        psTrace("ppStack", 1, "Photometering stacked image....\n");
-
-        psTimerStart("PPSTACK_PHOT");
-
-        fileActivation(config, combineFiles, false);
-        fileActivation(config, photFiles, true);
-        pmFPAview *photView = filesIterateDown(config);
-        fileActivation(config, combineFiles, true);
-
-        if (!ppStackPhotometry(config, outRO, photView)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
-            filesIterateUp(config);
-            psFree(outRO);
-            psFree(photView);
-            return false;
-        }
-        psFree(photView);
-
-        if (stats) {
-            pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File
-            pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
-            psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
-            psMetadataAddS32(stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", sources->n);
-            psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0,
-                             "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
-        }
-        psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT"));
-    }
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8 : Photometry Analysis : %f sec", psTimerClear("PPSTACK_STEPS"));
-
-    psThreadPoolFinalize();
-
-    memDump("phot");
-
-    // Statistics on output
-    if (stats) {
-        psTrace("ppStack", 1, "Gathering statistics on stacked image....\n");
-        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
-        psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
-
-        ppStatsFPA(stats, outRO->parent->parent->parent, view, maskBad, config);
-    }
-
-    memDump("stats");
-
-    // Put version information into the header
-    pmHDU *hdu = pmHDUFromCell(outRO->parent);
-    if (!hdu) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
-        psFree(outRO);
-        psFree(view);
+    // Photometry
+    psTrace("ppStack", 1, "Photometering stacked image....\n");
+    if (!ppStackPhotometry(options, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
+        psFree(options);
         return false;
     }
-    if (!hdu->header) {
-        hdu->header = psMetadataAlloc();
-    }
-    ppStackVersionHeader(hdu->header);
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("photometry");
 
-    psFree(outRO);
-    psFree(view);
 
-    // Write out summary statistics
-    if (stats) {
-        psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_STACK", 0,
-                         "Time to do photometry", psTimerClear("PPSTACK_TOTAL"));
-
-        const char *statsMDC = psMetadataConfigFormat(stats);
-        if (!statsMDC || strlen(statsMDC) == 0) {
-            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
-        } else {
-            fprintf(statsFile, "%s", statsMDC);
-        }
-        psFree((void *)statsMDC);
-        fclose(statsFile);
-
-        psFree(stats);
-    }
-
-    // Write out the output files
-    if (!filesIterateUp(config)) {
+    // Finish up
+    psTrace("ppStack", 1, "Finishing up....\n");
+    if (!ppStackFinish(options, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to finish up.");
+        psFree(options);
         return false;
     }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 9: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("finish");
 
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 9 : Final Output : %f sec", psTimerClear("PPSTACK_STEPS"));
-
-    memDump("finish");
-
+    psFree(options);
     return true;
 }

trunk/ppStack/src/ppStackPhotometry.c

@@ -9 +9 @@
 
 #include "ppStack.h"
+#include "ppStackLoop.h"
 
-#define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
-                          PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
-                          PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
-
-bool ppStackInputPhotometry(const pmReadout *ro, const psArray *sources, const pmConfig *config)
+bool ppStackPhotometry(ppStackOptions *options, pmConfig *config)
 {
-    PM_ASSERT_READOUT_NON_NULL(ro, false);
-    PS_ASSERT_ARRAY_NON_NULL(sources, false);
+    psAssert(options, "Require options");
+    psAssert(config, "Require configuration");
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
@@ -23 +20 @@
 
     bool mdok;                          // Status of MD lookup
-
-    float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement
-    if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe");
-        return false;
-    }
-    float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry
-    if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe");
-        return false;
-    }
-    float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry
-    if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe");
-        return false;
+    if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
+        // Nothing to do
+        return true;
     }
 
-    psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in
-    if (!mdok || !maskValStr) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe");
-        return false;
-    }
-    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
+    psTimerStart("PPSTACK_PHOT");
 
-    psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout
+    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, false);
+    ppStackFileActivation(config, PPSTACK_FILES_PHOT, true);
+    pmFPAview *photView = ppStackFilesIterateDown(config); // View to readout
+    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true);
 
-    // Measure background
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
-    if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
-        psFree(stats);
-        psFree(rng);
-        return false;
-    }
-    float bg = stats->robustMedian; // Background level
-    psFree(stats);
-    psFree(rng);
-
-    // Photometer sources
-    int numCols = image->numCols, numRows = image->numRows; // Size of image
-    int numSources = sources->n; // Number of sources
-    int numGood = 0;            // Number of good sources
-    float maxMag = -INFINITY;   // Maximum magnitude
-    for (int j = 0; j < numSources; j++) {
-        pmSource *source = sources->data[j]; // Source of interest
-        if (source->mode & PHOT_SOURCE_MASK || !isfinite(source->psfMag)) {
-            source->psfMag = NAN;
-            continue;
-        }
-        float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
-        int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel
-        // Range of integration
-        int xMin = PS_MAX(0, xCentral - zpRadius), xMax = PS_MIN(numCols - 1, xCentral + zpRadius);
-        int yMin = PS_MAX(0, yCentral - zpRadius), yMax = PS_MIN(numRows - 1, yCentral + zpRadius);
-
-        // Integrate
-        double sumImage = 0.0, sumKernel = 0.0; // Sums from integration
-        int numBadPix = 0;         // Number of bad pixels
-        for (int v = yMin; v <= yMax; v++) {
-            float dy2 = PS_SQR(y - v); // Distance from centroid
-            for (int u = xMin; u <= xMax; u++) {
-                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[v][u] & maskVal) {
-                    numBadPix++;
-                    continue;
-                }
-                float dx2 = PS_SQR(x - u); // Distance from centroid
-                double kernel = exp(-0.5 * (dx2 + dy2) / PS_SQR(zpSigma)); // Kernel value
-                sumImage += (image->data.F32[v][u] - bg) * kernel;
-                sumKernel += kernel;
-            }
-        }
-
-        if (numBadPix > zpFrac * M_PI * PS_SQR(zpRadius) || !isfinite(sumImage)) {
-            source->psfMag = NAN;
-        } else {
-            source->psfMag = - 2.5 * log10(sumImage / sumKernel * M_PI * PS_SQR(zpRadius));
-            if (isfinite(source->psfMag)) {
-                numGood++;
-                maxMag = PS_MAX(maxMag, source->psfMag);
-            }
-        }
-    }
-    psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag);
-
-    return true;
-}
-
-
-
-bool ppStackPhotometry(pmConfig *config, const pmReadout *readout, const pmFPAview *view)
-{
-    pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-    pmFPACopy(photFile->fpa, readout->parent->parent->parent);
+    pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File for photometry
+    pmFPACopy(photFile->fpa, options->outRO->parent->parent->parent);
 
     if (psMetadataLookupBool(NULL, config->arguments, "-visual")) {
@@ -122 +39 @@
     }
 
+#if 0
     // Need to ensure aperture residual is not calculated
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
@@ -131 +49 @@
         item->data.B = false;
     }
+#endif
 
     // set maskValue and markValue in the psphot recipe
     psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
     psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking
-    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue);
-    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for mark", markValue);
+    psMetadataAddImageMask(recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE,
+                            "Bits to mask", maskValue);
+    psMetadataAddImageMask(recipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE,
+                           "Bits to use for mark", markValue);
 
-    // XXX replace with psphotReadoutKnownSources (config, view)
-    // XXX this function requires that we construct the input source list and place it on PSPHOT.INPUT.CMF
-    pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT");
-    psArray *inSources = psMetadataLookupPtr (NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
+    pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT");
+    psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
     if (!inSources) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find input sources");
+        psFree(photView);
         return false;
     }
 
-    if (!psphotReadoutKnownSources(config, view, inSources)) {
+    if (!psphotReadoutKnownSources(config, photView, inSources)) {
         // Clear the error, so that the output files are written.
         psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on stacked image.");
+        psFree(photView);
         return false;
     }
+    psFree(photView);
 
     if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
@@ -162 +84 @@
     pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
 
+    if (options->stats) {
+        pmReadout *photRO = pmFPAviewThisReadout(photView, photFile->fpa); // Readout with the sources
+        psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
+        psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0,
+                         "Number of sources detected", sources->n);
+        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0,
+                         "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT"));
+
     return true;
 }

trunk/ppStack/src/ppStackReadout.c

@@ -8 +8 @@
 #include <psphot.h>
 
+#include "ppStackThread.h"
 #include "ppStack.h"
 
@@ -16 +17 @@
     psArray *args = job->args;          // Arguments
     ppStackThread *thread = args->data[0]; // Thread
-    pmConfig *config = args->data[1];   // Configuration
-    pmReadout *outRO = args->data[2];   // Output readout
-    psVector *mask = args->data[3];     // Mask for inputs
-    psVector *weightings = args->data[4]; // Weightings (1/noise^2) for each image
-    psVector *addVariance = args->data[5]; // Additional variance when rejecting
+    ppStackOptions *options = args->data[1]; // Options
+    pmConfig *config = args->data[2];   // Configuration
+
+    pmReadout *outRO = options->outRO;  // Output readout
+    psVector *mask = options->inputMask; // Mask for inputs
+    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
+    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
 
     psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
@@ -37 +40 @@
     psArray *args = job->args;          // Arguments
     ppStackThread *thread = args->data[0]; // Thread
-    pmConfig *config = args->data[1];   // Configuration
-    pmReadout *outRO = args->data[2];   // Output readout
-    psVector *mask = args->data[3];     // Mask for inputs
-    psArray *rejected = args->data[4];  // Rejected pixels
-    psVector *weightings = args->data[5];  // Weighting (1/noise^2) for each image
-    psVector *addVariance = args->data[6];  // Weighting (1/noise^2) for each image
+    ppStackOptions *options = args->data[1]; // Options
+    pmConfig *config = args->data[2];   // Configuration
+
+    pmReadout *outRO = options->outRO;  // Output readout
+    psVector *mask = options->inputMask; // Mask for inputs
+    psArray *rejected = options->rejected; // Rejected pixels
+    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
+    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
 
     bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected,

trunk/ppStack/src/ppStackThread.c

@@ -9 +9 @@
 
 #include "ppStack.h"
+#include "ppStackOptions.h"
+#include "ppStackThread.h"
+
 
 #define THREAD_WAIT 10000               // Microseconds to wait if thread is not available
@@ -53 +56 @@
 }
 
-ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, const psArray *imageNames,
-                                          const psArray *maskNames, const psArray *varianceNames,
-                                          const psArray *covariances, const pmConfig *config)
-{
+ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config)
+{
+    psAssert(options, "Require options");
+    psAssert(config, "Require configuration");
+
+    const psArray *cells = options->cells; // Array of input cells
+    const psArray *imageNames = options->imageNames; // Names of images to read
+    const psArray *maskNames = options->maskNames; // Names of masks to read
+    const psArray *varianceNames = options->varianceNames; // Names of variance maps to read
+    const psArray *covariances = options->covariances; // Covariance matrices (already read)
+
     PS_ASSERT_ARRAY_NON_NULL(cells, NULL);
     PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL);
@@ -245 +255 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 6);
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 3);
         task->function = &ppStackReadoutInitialThread;
         psThreadTaskAdd(task);
@@ -259 +269 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7);
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 3);
         task->function = &ppStackReadoutFinalThread;
         psThreadTaskAdd(task);