#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <unistd.h>
#include <string.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.WEIGHT", NULL };

// Output files for the combination
static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", 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;
}

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


bool ppStackLoop(pmConfig *config)
{
    assert(config);

    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    psAssert(recipe, "We've thrown an error on this before.");

    bool mdok;                          // Status of MD lookup
    bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files?
    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 *tempWeight = psMetadataLookupStr(NULL, recipe, "TEMP.WEIGHT"); // Suffix for temp weight maps
    if (!tempImage || !tempMask || !tempWeight) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.WEIGHT in recipe");
        return false;
    }

    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
    int numScans = psMetadataLookupS32(NULL, recipe, "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

    if (!pmConfigMaskSetBits(NULL, NULL, config)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
        return false;
    }

    memDump("start");

    // Preparation iteration: Load the sources, and get a target PSF model
    psTrace("ppStack", 1, "Determining target PSF....\n");
    psArray *globalSources = NULL;      // Global list of sources for matching (haveSources = TRUE)
    psArray *indSources = psArrayAlloc(num); // Individual lists of sources for matching (haveSources = FALSE)
    pmPSF *targetPSF = NULL;            // Target PSF
    bool haveSources = psMetadataLookupBool(NULL, config->arguments, "HAVE.SOURCES"); // Global sources?
    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
        psArray *sourceLists = NULL;    // Source lists for merging sources from multiple readouts
        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(globalSources);
                psFree(indSources);
                psFree(sourceLists);
                psFree(fileIter);
                psFree(psfs);
                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(globalSources);
                psFree(indSources);
                psFree(sourceLists);
                psFree(fileIter);
                psFree(psfs);
                return false;
            }
            if (numCols == 0 && numRows == 0) {
                numCols = naxis1;
                numRows = naxis2;
            }

            if (!haveSources) {
                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;
                }
                indSources->data[index] = psMemIncrRefCounter(sources);
#ifdef TESTING
                ppStackSourcesPrint(sources);
#endif
                sourceLists = ppStackSourceListAdd(sourceLists, sources, config);
                if (!sourceLists) {
                    psError(PS_ERR_UNKNOWN, false, "Unable to add sources to list.");
                    psFree(view);
                    psFree(globalSources);
                    psFree(indSources);
                    psFree(fileIter);
                    psFree(psfs);
                    return false;
                }
            }

            index++;
        }
        psFree(fileIter);

        targetPSF = ppStackPSF(config, numCols, numRows, psfs);
        psFree(psfs);
        if (!targetPSF) {
            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
            psFree(globalSources);
            psFree(indSources);
            psFree(view);
            return false;
        }

        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;

        if (haveSources) {
            // We want to hang on to the 'sources' even when its host FPA is blown away
            pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.INPUT.SOURCES");
            globalSources = psMetadataLookupPtr(NULL, sourcesRO->analysis, "PSPHOT.SOURCES");
            psMemIncrRefCounter(globalSources);
            if (!globalSources) {
                psError(PS_ERR_UNKNOWN, true, "Unable to find sources.");
                psFree(view);
                return false;
            }
        } else {
            globalSources = ppStackSourceListCombine(sourceLists, config);
            psFree(sourceLists);
            if (!globalSources) {
                psError(PS_ERR_UNKNOWN, false, "Unable to add sources to list.");
                psFree(view);
                psFree(psfs);
                return false;
            }
#ifdef TESTING
            ppStackSourcesPrint(globalSources);
#endif
        }

        psFree(view);
        filesIterateUp(config);
    }

    memDump("psf");

    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.%d.%s", outName, i, tempImage);
        psStringAppend(&maskName, "%s.%d.%s", outName, i, tempMask);
        psStringAppend(&weightName, "%s.%d.%s", outName, i, tempWeight);
        imageNames->data[i] = imageName;
        maskNames->data[i] = maskName;
        weightNames->data[i] = weightName;
    }

    // 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, 0); // Random number generator
    int numGood = 0;                    // Number of good frames
    int numCols = 0, numRows = 0;       // Size of image
    for (int i = 0; i < num; i++) {
        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(globalSources);
            psFree(indSources);
            psFree(targetPSF);
            psFree(rng);
            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(globalSources);
            psFree(indSources);
            psFree(targetPSF);
            psFree(rng);
            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
        psArray *sources = haveSources ? globalSources : indSources->data[i]; // Sources for matching
        if (!ppStackMatch(readout, &regions, &kernels, sources, targetPSF, rng, config)) {
            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
            psErrorClear();
            continue;
        }

        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);
        writeImage(maskNames->data[i],   hdu->header, readout->mask, config);
        writeImage(weightNames->data[i], hdu->header, readout->weight, config);

        cells->data[i] = psMemIncrRefCounter(readout->parent);
        filesIterateUp(config);
        numGood++;

        memDump("match");
    }
    psFree(globalSources);
    psFree(indSources);
    psFree(targetPSF);
    psFree(rng);

    if (numGood == 0) {
        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
        return false;
    }


#ifdef TESTING
    psTraceSetLevel("psModules.imcombine", 7);
#endif

    memDump("preinitial");

    // Stack the convolved files
    psTrace("ppStack", 1, "Initial stack of convolved images....\n");
    {
        pmFPAfileActivate(config->files, false, NULL);
        fileActivation(config, combineFiles, 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++) {
            pmCell *cell = cells->data[i]; // Cell of interest
            if (!cell) {
                // Bad image
                continue;
            }
            pmReadout *ro = cell->readouts->data[0]; // Readout of interest
            if (!ro) {
                ro = pmReadoutAlloc(cell);
            }
            readouts->data[i] = ro; // 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++) {
            if (!readouts->data[i]) {
                // Bad image
                continue;
            }
            // Resolved names
            psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
            psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
            psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false);
            imageFits->data[i] = psFitsOpen(imageResolved, "r");
            maskFits->data[i] = psFitsOpen(maskResolved, "r");
            weightFits->data[i] = psFitsOpen(weightResolved, "r");
            psFree(imageResolved);
            psFree(maskResolved);
            psFree(weightResolved);
            if (!imageFits->data[i] || !maskFits->data[i] || !weightFits->data[i]) {
                psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s",
                        (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)weightNames->data[i]);
                psFree(subKernels);
                psFree(subRegions);
                return false;
            }
        }

        // Read convolutions by chunks
        bool more = true;               // More to read?
        for (int numChunk = 0; more; numChunk++) {
            psTrace("ppStack", 2, "Initial stack of chunk %d....\n", numChunk);
            for (int i = 0; i < num; i++) {
                pmReadout *readout = readouts->data[i];
                if (!readout) {
                    // Bad image
                    continue;
                }

                if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, overlap,
                                        config) ||
                    !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, overlap,
                                            config) ||
                    !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, overlap,
                                              config)) {
                    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;
                }
            }

#ifndef PS_NO_TRACE
            {
                pmReadout *ro = NULL;   // Representative readout
                for (int i = 0; i < num && !ro; i++) {
                    ro = readouts->data[i];
                }
                psAssert(ro, "There should be a readout here.");
                psTrace("ppStack", 4, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols,
                        ro->row0, ro->row0 + ro->image->numRows);
            }
#endif

            if (!ppStackReadoutInitial(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;
            }

            for (int i = 0; i < num && more; i++) {
                pmReadout *readout = readouts->data[i];
                if (!readout) {
                    // Bad image
                    continue;
                }
                more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans, config);
                more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans, config);
                more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans, config);
            }
        }

        memDump("initial");

        // Reset for the second read
        // Extract the rejection lists
        psArray *rejected = psArrayAlloc(num); // Pixels to inspect
        int numRejected = 0;            // Number of inputs rejected completely
        for (int i = 0; i < num; i++) {
            pmReadout *ro = readouts->data[i]; // Readout of interest
            if (!ro) {
                // Bad image
                numRejected++;
                continue;
            }
            psPixels *inspect = psPixelsAllocEmpty(PIXELS_BUFFER); // Inspection list for this readout
            psMetadataIterator *iter = psMetadataIteratorAlloc(ro->analysis, PS_LIST_HEAD,
                                                               "^" PPSTACK_INSPECT_PIXELS "$"); // Iterator
            psMetadataItem *item;
            while ((item = psMetadataGetAndIncrement(iter))) {
                psPixels *pixels = item->data.V; // Rejected pixels
                if (!pixels) {
                    continue;
                }
                psTrace("ppStack", 5, "Adding %ld pixels to inspect to image %d", pixels->n, i);
                inspect = psPixelsConcatenate(inspect, pixels);
            }
            psFree(iter);
            psMetadataRemoveKey(ro->analysis, PPSTACK_INSPECT_PIXELS);
            pmReadoutFreeData(ro);

            psLogMsg("ppStack", PS_LOG_INFO, "%ld total pixels to inspect from image %d", inspect->n, i);

#ifdef TESTING
            {
                psImage *mask = psPixelsToMask(NULL, inspect, psRegionSet(0, numCols - 1, 0, numRows - 1),
                                               0xff); // Mask image
                psString name = NULL;           // Name of image
                psStringAppend(&name, "inspect_%03d.fits", i);
                psFits *fits = psFitsOpen(name, "w");
                psFree(name);
                psFitsWriteImage(fits, NULL, mask, 0, NULL);
                psFree(mask);
                psFitsClose(fits);
            }
#endif

            psPixels *reject = pmStackReject(inspect, numCols, numRows, threshold, poorFrac,
                                             subRegions->data[i], subKernels->data[i]); // Pixels to reject

#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);
                psFits *fits = psFitsOpen(name, "w");
                psFree(name);
                psFitsWriteImage(fits, NULL, mask, 0, NULL);
                psFree(mask);
                psFitsClose(fits);
            }
#endif

            psFree(inspect);
            if (!reject) {
                psWarning("Rejection on image %d didn't work --- reject entire image.", i);
                numRejected++;
            } else {
                float frac = reject->n / (float)(outRO->image->numCols * outRO->image->numRows); // Pixel frac
                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;
                    numRejected++;
                }
            }
            rejected->data[i] = reject;

            memDump("reject");

        }
        psFree(subKernels);
        psFree(subRegions);

        if (numRejected == num) {
            psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed.");
            psFree(readouts);
            psFree(rejected);
            psFree(outRO);
            psFree(view);
            return false;
        }

        memDump("prefinal");

        // Read convolutions by chunks
        psTrace("ppStack", 2, "Final stack of convolved images....\n");
        more = true;
        for (int numChunk = 0; more; numChunk++) {
            psTrace("ppStack", 2, "Final stack of chunk %d....\n", numChunk);
            for (int i = 0; i < num; i++) {
                if (!rejected->data[i]) {
                    continue;
                }
                pmReadout *readout = readouts->data[i];
                assert(readout);

                if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, 0, config) ||
                    !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, 0, config) ||
                    !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, 0,
                                              config)) {
                    psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
                    psFree(readouts);
                    psFree(rejected);
                    psFree(outRO);
                    psFree(view);
                    return false;
                }
            }

#ifndef PS_NO_TRACE
            {
                pmReadout *ro = NULL;
                for (int i = 0; i < num && !ro; i++) {
                    if (!rejected->data[i]) {
                        continue;
                    }
                    ro = readouts->data[i];
                }
                if (ro && ro->image) {
                    psTrace("ppStack", 4, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols,
                            ro->row0, ro->row0 + ro->image->numRows);
                }
            }
#endif

            if (!ppStackReadoutFinal(config, outRO, readouts, rejected)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
                psFree(readouts);
                psFree(rejected);
                psFree(outRO);
                psFree(view);
                return false;
            }

            for (int i = 0; i < num && more; i++) {
                pmReadout *readout = readouts->data[i];
                if (!readout) {
                    continue;
                }
                more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans, config);
                more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans, config);
                more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans, config);
            }
        }

        psFree(rejected);
        for (int i = 0; i < readouts->n; i++) {
            pmReadout *ro = readouts->data[i]; // Readout of interest
            if (!ro) {
                continue;
            }
            pmReadoutFreeData(ro);
        }
        psFree(readouts);

        // Ensure masked regions really look masked
        {
            psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
            psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
            if (!pmReadoutMaskApply(outRO, maskBad)) {
                psWarning("Unable to apply mask");
            }
        }

        // Close up
        bool wcsDone = false;           // Have we done the WCS?
        for (int i = 0; i < num; i++) {
            if (!readouts->data[i]) {
                continue;
            }
            if (!wcsDone) {
                // Copy astrometry over
                wcsDone = true;
                pmReadout *inRO = readouts->data[i]; // Template readout
                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;
                        }
                    }
                }
            }

            psFitsClose(imageFits->data[i]);
            psFitsClose(maskFits->data[i]);
            psFitsClose(weightFits->data[i]);
            imageFits->data[i] = NULL;
            maskFits->data[i] = NULL;
            weightFits->data[i] = NULL;
            if (tempDelete) {
                psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
                psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
                psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false);
                if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
                    unlink(weightResolved) == -1) {
                    psWarning("Unable to delete temporary files for image %d", i);
                }
                psFree(imageResolved);
                psFree(maskResolved);
                psFree(weightResolved);
            }
        }
        psFree(imageNames);
        psFree(maskNames);
        psFree(weightNames);
        psFree(imageFits);
        psFree(maskFits);
        psFree(weightFits);

        memDump("final");

        if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {

            fileActivation(config, combineFiles, false);
            fileActivation(config, photFiles, true);
            pmFPAview *photView = filesIterateDown(config);

            psTrace("ppStack", 1, "Photometering stacked image....\n");
            if (!ppStackPhotometry(config, outRO, view)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
                psFree(outRO);
                psFree(photView);
                psFree(view);
                return false;
            }
            psFree(photView);

            fileActivation(config, combineFiles, true);
        }

        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
            psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels

            ppStatsFPA(stats, outCell->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.");
            return false;
        }
        if (!hdu->header) {
            hdu->header = psMetadataAlloc();
        }
        ppStackVersionMetadata(hdu->header);

        // Write out the output files
        filesIterateUp(config);

        psFree(outRO);
        psFree(view);
    }

    // Write out summary statistics
    if (stats) {
        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);
    }

    memDump("finish");

    return true;
}
