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

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


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

    psTimerStart("PPSTACK_TOTAL");

    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 *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");
        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.");
        return false;
    }

    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");
        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

    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
    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);
                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);
                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);
                    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);
                    return false;
                }

                psFree(view);
                if (!filesIterateUp(config)) {
                    psFree(sourceLists);
                    psFree(targetPSF);
                    return false;
                }
                pmFPAfileActivate(config->files, false, NULL);
                fileActivation(config, prepareFiles, true);
            }

            index++;
        }
        psFree(fileIter);

        // Zero point calibration
        sumExposure = ppStackSourcesTransparency(sourceLists, view, config);
        if (!isfinite(sumExposure) || sumExposure <= 0) {
            psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
            psFree(sourceLists);
            psFree(targetPSF);
            return false;
        }

        // Generate target PSF
        targetPSF = ppStackPSF(config, numCols, numRows, psfs);
        psFree(psfs);
        if (!targetPSF) {
            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
            psFree(sourceLists);
            psFree(view);
            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, 0); // Random number generator
    int numGood = 0;                    // Number of good frames
    int numCols = 0, numRows = 0;       // Size of image
    psVector *inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
    psVectorInit(inputMask, 0);
    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++) {
        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);
            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);
#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);
                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);
                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;
            }

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

    memDump("final");

#ifdef TESTING
    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);
        return false;
    }
    if (!hdu->header) {
        hdu->header = psMetadataAlloc();
    }
    ppStackVersionMetadata(hdu->header);

    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)) {
        return false;
    }

    psLogMsg("ppStack", PS_LOG_INFO, "Stage 9 : Final Output : %f sec", psTimerClear("PPSTACK_STEPS"));

    memDump("finish");

    return true;
}
