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

#include <stdio.h>
#include <pslib.h>
#include <psmodules.h>

#include "ppStack.h"

#define ARRAY_BUFFER 16                 // Number to add to array at a time
#define MAG_IGNORE 50                   // Ignore magnitudes fainter than this --- they're not real!
#define FAKE_SIZE 1                     // Size of fake convolution kernel
#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
                     PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
#define FAINT_SOURCE_FRAC 1.0e-4         // Set minimum flux to this fraction of faintest source flux

#define TESTING                         // Enable debugging output

#ifdef TESTING
// Read a FITS image
static bool readImage(psImage **target, // Target for image
                      const char *name, // Name of FITS file
                      const pmConfig *config // Configuration
    )
{
    psString resolved = pmConfigConvertFilename(name, config, false, false); // Resolved filename
    psFits *fits = psFitsOpen(resolved, "r");
    psFree(resolved);
    if (!fits) {
        psError(PS_ERR_IO, false, "Unable to open previously produced image: %s", name);
        return false;
    }
    psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
    if (!image) {
        psError(PS_ERR_IO, false, "Unable to read previously produced image: %s", name);
        psFitsClose(fits);
        return false;
    }
    psFitsClose(fits);

    psFree(*target);
    *target = image;

    return true;
}
#endif

static psArray *stackSourcesFilter(psArray *sources, // Source list to filter
                                   int exclusion // Exclusion zone, pixels
    )
{
    psAssert(sources && sources->n > 0, "Require array of sources");
    if (exclusion <= 0) {
        return psMemIncrRefCounter(sources);
    }

    int num = sources->n;               // Number of sources
    psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates
    int numGood = 0;                    // Number of good sources
    for (int i = 0; i < num; i++) {
        pmSource *source = sources->data[i]; // Source of interest
        if (!source) {
            continue;
        }
        x->data.F32[numGood] = source->peak->xf;
        y->data.F32[numGood] = source->peak->yf;
        numGood++;
    }
    x->n = y->n = numGood;

    psTree *tree = psTreePlant(2, 2, x, y); // kd tree

    psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources
    psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source
    int numFiltered = 0;                // Number of filtered sources
    for (int i = 0; i < num; i++) {
        pmSource *source = sources->data[i]; // Source of interest
        if (!source) {
            continue;
        }
        coords->data.F64[0] = source->peak->xf;
        coords->data.F64[1] = source->peak->yf;

        long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
        psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
                coords->data.F64[0], coords->data.F64[1], numWithin);
        if (numWithin == 1) {
            // Only itself inside the exclusion zone
            filtered = psArrayAdd(filtered, filtered->n, source);
        } else {
            numFiltered++;
        }
    }
    psFree(coords);
    psFree(tree);

    psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);

    return filtered;
}

#define BG_SIZE 64                      // Large box half-size in which to measure background

// Generate a background model of the readout we're matching
// We only need to set the background around the sources, not everywhere.
psImage *stackBackgroundModel(pmReadout *ro, psMaskType maskVal, const psArray *sources, int size)
{
    psImage *image = ro->image, *mask = ro->mask; // Image and mask of readout
    int numCols = image->numCols, numRows = image->numRows; // Size of image

    psImage *bgImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Background image
    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for measuring background
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator

    for (int i = 0; i < sources->n; i++) {
        pmSource *source = sources->data[i]; // Source of interest
        if (!source) {
            continue;
        }

        int x = source->peak->xf + 0.5, y = source->peak->yf + 0.5; // Coordinates of source

        int xMin = PS_MAX(x - BG_SIZE, 0), xMax = PS_MIN(x + BG_SIZE, numCols);
        int yMin = PS_MAX(y - BG_SIZE, 0), yMax = PS_MIN(y + BG_SIZE, numCols);

        psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
        psImage *subImage = psImageSubset(image, region); // Subimage
        psImage *subMask = psImageSubset(mask, region); // Subimage mask
        if (!psImageBackground(stats, NULL, subImage, subMask, maskVal, rng)) {
            // Can't do anything about it
            psErrorClear();
            continue;
        }
        psFree(subImage);
        psFree(subMask);

        float value = stats->robustMedian; // Value to set background to

        int uMin = PS_MAX(x - size, 0), uMax = PS_MIN(x + size, numCols);
        int vMin = PS_MAX(y - size, 0), vMax = PS_MIN(y + size, numCols);

        for (int v = vMin; v < vMax; v++) {
            for (int u = uMin; u < uMax; u++) {
                bgImage->data.F32[v][u] = value;
            }
        }
    }

    psFree(stats);
    psFree(rng);

    return bgImage;
}

bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
                  psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
{
    assert(readout);
    assert(regions && !*regions);
    assert(kernels && !*kernels);
    assert(config);

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

    // Look up appropriate values from the ppSub recipe
    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size

    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
    psMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
    psMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels

    bool mdok;                          // Status of MD lookup
    bool renorm = psMetadataLookupBool(&mdok, recipe, "RENORM"); // Renormalise variances?
    float penalty = psMetadataLookupF32(NULL, ppsub, "PENALTY"); // Penalty for wideness
    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads

    if (!pmReadoutMaskNonfinite(readout, maskVal)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in readout.");
        return false;
    }

    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily

    static int numInput = -1;            // Index of input file
    numInput++;
#ifdef TESTING
    // Read previously produced kernel
    if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
        const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
        assert(outName);
        // Read convolution kernel
        {
            psString filename = NULL;   // Output filename
            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
            psFree(filename);
            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
            psFree(resolved);
            if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
                psFitsClose(fits);
                return false;
            }
            psFitsClose(fits);

            // Add in variance factor
            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
            if (!isfinite(vf)) {
                vf = 1.0;
            }
            if (isfinite(vfItem->data.F32)) {
                vfItem->data.F32 *= vf;
            } else {
                vfItem->data.F32 = vf;
            }
        }

        // Read image, mask, weight
        const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for image
        const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for mask
        const char *tempWeight = psMetadataLookupStr(NULL, recipe, "TEMP.WEIGHT"); // Suffix for weight map
        psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images
        psStringAppend(&imageName, "%s.%d.%s", outName, numInput, tempImage);
        psStringAppend(&maskName, "%s.%d.%s", outName, numInput, tempMask);
        psStringAppend(&weightName, "%s.%d.%s", outName, numInput, tempWeight);

        if (!readImage(&readout->image, imageName, config) || !readImage(&readout->mask, maskName, config) ||
            !readImage(&readout->weight, weightName, config)) {
            psError(PS_ERR_IO, false, "Unable to read previously produced image.");
            psFree(imageName);
            psFree(maskName);
            psFree(weightName);
            return false;
        }
        psFree(imageName);
        psFree(maskName);
        psFree(weightName);
    } else {
#endif

        // Normal operations here
        if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) {
            assert(psf);
            assert(sources);

            int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
            float regionSize = psMetadataLookupF32(NULL, ppsub, "REGION.SIZE"); // Size of iso-kernel regs
            float spacing = psMetadataLookupF32(NULL, ppsub, "STAMP.SPACING"); // Typical stamp spacing
            int footprint = psMetadataLookupS32(NULL, ppsub, "STAMP.FOOTPRINT"); // Stamp half-size
            float threshold = psMetadataLookupF32(NULL, ppsub, "STAMP.THRESHOLD"); // Threshold for stmps
            int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches
            int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations
            float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
            float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel
            pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(
                psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE")); // Kernel type
            psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths
            psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders
            int inner = psMetadataLookupS32(NULL, ppsub, "INNER"); // Inner radius
            int ringsOrder = psMetadataLookupS32(NULL, ppsub, "RINGS.ORDER"); // RINGS polynomial order
            int binning = psMetadataLookupS32(NULL, ppsub, "SPAM.BINNING"); // Binning for SPAM kernel
            float badFrac = psMetadataLookupF32(NULL, ppsub, "BADFRAC"); // Maximum bad fraction
            bool optimum = psMetadataLookupBool(&mdok, ppsub, "OPTIMUM"); // Derive optimum parameters?
            float optMin = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MIN"); // Minimum width for search
            float optMax = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MAX"); // Maximum width for search
            float optStep = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.STEP"); // Step for search
            float optThresh = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.TOL"); // Tolerance for search
            int optOrder = psMetadataLookupS32(&mdok, ppsub, "OPTIMUM.ORDER"); // Order for search
            float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor"

            // These values are specified specifically for stacking
            const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename

            psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
            if (optimum) {
                optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
            }

            // Generate a fake image to match to
            float minFlux = INFINITY;       // Minimum flux for fake image
            {
                psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for bg
                if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal, rng)) {
                    psWarning("Can't measure background for image.");
                    psErrorClear();
                } else {
                    minFlux = bg->robustStdev;
                }
                psFree(bg);
            }

#if 0
            float maxMag = -INFINITY;   // Maximum magnitude of sources
            for (int i = 0; i < sources->n; i++) {
                pmSource *source = sources->data[i]; // Source of interest
                if (source->psfMag > maxMag && source->psfMag <= MAG_IGNORE &&
                    !(source->mode & SOURCE_MASK)) {
                    maxMag = source->psfMag;
                }
            }
            minFlux = PS_MIN(FAINT_SOURCE_FRAC * powf(10.0, -0.4 * maxMag), minFlux);
#endif


#if 0
            // Testing the normalisation of the fake image
            {
                pmReadout *test = pmReadoutAlloc(NULL); // Test readout
                psArray *sources = psArrayAlloc(1);     // Array of sources
                pmSource *source = pmSourceAlloc();     // Source
                sources->data[0] = source;
                source->peak = pmPeakAlloc(50, 50, 10000, PM_PEAK_LONE);
                source->psfMag = -13.0;
                pmReadoutFakeFromSources(test, 100, 100, sources, NULL, NULL, psf, 0.1, 0, false, true);
                float sum = 0.0;
                for (int y = 0; y < test->image->numRows; y++) {
                    for (int x = 0; x < test->image->numCols; x++) {
                        sum += test->image->data.F32[y][x];
                    }
                }
                fprintf(stderr, "Photometry: %f --> %f = -13.0 ???\n", sum, -2.5*log10(sum));

                psFits *fits = psFitsOpen("testphot.fits", "w");
                psFitsWriteImage(fits, NULL, test->image, 0, NULL);
                psFitsClose(fits);
                exit(0);
            }
#endif

            pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF

            // For the sake of stamps, remove nearby sources
            psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources

            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
                                          stampSources, NULL, NULL, psf, minFlux, 0, false, true)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
                psFree(fake);
                psFree(optWidths);
                psFree(output);
                return false;
            }

            // Add the background into the target image
            psImage *bgImage = stackBackgroundModel(readout, maskVal, stampSources, footprint); // Image of background
            psBinaryOp(fake->image, fake->image, "+", bgImage);
            psFree(bgImage);

#ifdef TESTING
            {
                psString name = NULL;
                psStringAppend(&name, "fake_%03d.fits", numInput);
                psFits *fits = psFitsOpen(name, "w");
                psFree(name);
                psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
                psFitsClose(fits);
            }
            {
                psString name = NULL;
                psStringAppend(&name, "real_%03d.fits", numInput);
                psFits *fits = psFitsOpen(name, "w");
                psFree(name);
                psFitsWriteImage(fits, NULL, readout->image, 0, NULL);
                psFitsClose(fits);
            }
#endif

            // Renormalise the variances if desired
            if (renorm) {
                // Statistics for renormalisation
                psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
                                                                                        "RENORM.MEAN"));
                psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
                                                                                         "RENORM.STDEV"));

                if (!pmReadoutWeightRenormPixels(readout, maskBad, renormMean, renormStdev, rng)) {
                    psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
                    psFree(output);
                    return false;
                }
            }

            if (threads > 0) {
                pmSubtractionThreadsInit(readout, fake);
            }

            // Do the image matching
            if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, stride, regionSize, spacing,
                                    threshold, stampSources, stampsName, type, size, order, widths, orders,
                                    inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder,
                                    optThresh, iter, rej, sysError, maskVal, maskBad, maskPoor, poorFrac,
                                    badFrac, PM_SUBTRACTION_MODE_1)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
                psFree(fake);
                psFree(optWidths);
                psFree(stampSources);
                psFree(output);
                return false;
            }

#ifdef TESTING
            {
                psString name = NULL;
                psStringAppend(&name, "conv_%03d.fits", numInput);
                psFits *fits = psFitsOpen(name, "w");
                psFree(name);
                psFitsWriteImage(fits, NULL, output->image, 0, NULL);
                psFitsClose(fits);
            }
            {
                psString name = NULL;
                psStringAppend(&name, "diff_%03d.fits", numInput);
                psFits *fits = psFitsOpen(name, "w");
                psFree(name);
                psBinaryOp(fake->image, output->image, "-", fake->image);
                psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
                psFitsClose(fits);
            }
#endif

            psFree(fake);
            psFree(optWidths);
            psFree(stampSources);

            if (threads > 0) {
                pmSubtractionThreadsFinalize(readout, fake);
            }

            // Set the variance factor
            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
            float vf = psMetadataLookupF32(NULL, output->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
            if (!isfinite(vf)) {
                vf = 1.0;
            }
            if (isfinite(vfItem->data.F32)) {
                vfItem->data.F32 *= vf;
            } else {
                vfItem->data.F32 = vf;
            }

            // Replace original images with convolved
            psFree(readout->image);
            psFree(readout->mask);
            psFree(readout->weight);
            readout->image  = psMemIncrRefCounter(output->image);
            readout->mask   = psMemIncrRefCounter(output->mask);
            readout->weight = psMemIncrRefCounter(output->weight);
        } else {
            // Fake the convolution
            psRegion *region = psRegionAlloc(0, readout->image->numCols - 1, 0, readout->image->numRows - 1);
            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
                             PS_DATA_REGION | PS_META_DUPLICATE_OK, "Fake subtraction region", region);
            psFree(region);
            pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(FAKE_SIZE, 0, penalty,
                                                                     PM_SUBTRACTION_MODE_1);
            // Set solution to delta function
            kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64);
            psVectorInit(kernels->solution1, 0.0);
            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
            kernels->solution1->data.F64[normIndex] = 1.0;
            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
                             PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Fake subtraction kernel", kernels);
            psFree(kernels);
        }

#ifdef TESTING
        // Write convolution kernel
        {
            const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
            assert(outName);

            psString filename = NULL;   // Output filename
            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
            psString resolved = pmConfigConvertFilename(filename, config, true, false); // Resolved filename
            psFree(filename);
            psFits *fits = psFitsOpen(resolved, "w"); // FITS file for subtraction kernel
            psFree(resolved);
            pmReadoutWriteSubtractionKernels(output, fits);
            psFitsClose(fits);
        }
    }
#endif

    readout->analysis = psMetadataCopy(readout->analysis, output->analysis);

// Extract the regions and solutions used in the image matching
// This stops them from being freed when we iterate back up the FPA
    *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
    {
        psString regex = NULL;          // Regular expression
        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
        psFree(regex);
        psMetadataItem *item = NULL;// Item from iteration
        while ((item = psMetadataGetAndIncrement(iter))) {
            assert(item->type == PS_DATA_REGION);
            *regions = psArrayAdd(*regions, ARRAY_BUFFER, item->data.V);
        }
        psFree(iter);
    }
    *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels
    {
        psString regex = NULL;          // Regular expression
        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
        psFree(regex);
        psMetadataItem *item = NULL;// Item from iteration
        while ((item = psMetadataGetAndIncrement(iter))) {
            assert(item->type == PS_DATA_UNKNOWN);
            // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
            // images by scans.
            pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
            kernel->numCols = readout->image->numCols;
            kernel->numRows = readout->image->numRows;

            *kernels = psArrayAdd(*kernels, ARRAY_BUFFER, kernel);
        }
        psFree(iter);
    }
    assert((*regions)->n == (*kernels)->n);

    // Record chi^2
    {
        *chi2 = 0.0;
        int num = 0;                    // Number of measurements of chi^2
        psString regex = NULL;          // Regular expression
        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
        psFree(regex);
        psMetadataItem *item = NULL;// Item from iteration
        while ((item = psMetadataGetAndIncrement(iter))) {
            assert(item->type == PS_DATA_UNKNOWN);
            pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
            *chi2 += kernels->mean;
            num++;
        }
        psFree(iter);
        *chi2 /= num;
    }

    // Renormalise the variances if desired
    if (renorm) {
        // Statistics for renormalisation
        psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
                                                                                "RENORM.MEAN"));
        psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
                                                                                 "RENORM.STDEV"));

        if (!pmReadoutWeightRenormPixels(readout, maskBad, renormMean, renormStdev, rng)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
            psFree(output);
            return false;
        }
    }

    // Ensure the background value is zero
    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
    if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
        psWarning("Can't measure background for image.");
        psErrorClear();
    } else {
        psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
                 psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
        (void)psBinaryOp(readout->image, readout->image, "-",
                         psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
        *weighting = 1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
                         "Weighting by 1/noise^2 for stack",
                         1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)));
    }
    psFree(bg);

#define RADIUS 10                       // Radius of photometry
#define MIN_ERR 0.05                    // Minimum photometric error, mag

    // Ensure the normalisation is correct
    // XXX Ideally, would like to do proper PSF-fit photometry, but will settle for aperture photometry
    int numSources = sources->n;        // Number of sources
    psVector *ratio = psVectorAlloc(numSources, PS_TYPE_F32); // Flux ratios for sources
    psVector *ratioMask = psVectorAlloc(numSources, PS_TYPE_MASK); // Mask for flux ratios
    psVectorInit(ratioMask, 0xFF);
    psImage *image = readout->image;    // Image of interest
    psImage *mask = readout->mask;      // Mask of interest
    int numCols = image->numCols, numRows = image->numRows; // Size of image
    for (int i = 0; i < numSources; i++) {
        pmSource *source = sources->data[i]; // Source of interest
        if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) ||
            source->errMag > MIN_ERR) {
            continue;
        }

        float xSrc = source->peak->xf, ySrc = source->peak->yf; // Source coordinates
        int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x
        int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y
        int numPix = 0;                 // Number of pixels
        float sum = 0.0;                // Sum of pixels
        for (int y = yMin; y <= yMax; y++) {
            for (int x = xMin; x <= xMax; x++) {
                if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) {
                    continue;
                }
                sum += image->data.F32[y][x];
                numPix++;
            }
        }
        if (sum >= 0 && numPix > 0) {
            float mag = -2.5 * log10(sum * M_PI * PS_SQR(RADIUS) / numPix); // Instrumental magnitude
            ratio->data.F32[i] = mag - source->psfMag;
            ratioMask->data.PS_TYPE_MASK_DATA[i] = 0;
        }
    }

    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
    if (!psVectorStats(stats, ratio, NULL, ratioMask, 0xFF)) {
        psWarning("Unable to measure normalisation --- assuming correct.");
    } else {
        psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image by %f (+/- %f) mag\n",
                 stats->robustMedian, stats->robustStdev);
        float norm = powf(10.0, -0.4 * stats->robustMedian); // Normalisation to apply
        psBinaryOp(image, image, "*", psScalarAlloc(norm, PS_TYPE_F32));
    }
    psFree(stats);
    psFree(ratio);
    psFree(ratioMask);

#ifdef TESTING
    {
        psString name = NULL;
        psStringAppend(&name, "convolved_%03d.fits", numInput);
        psFits *fits = psFitsOpen(name, "w");
        psFree(name);
        psFitsWriteImage(fits, NULL, output->image, 0, NULL);
        psFitsClose(fits);
    }
#endif

    psFree(output);

    return true;
}
