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

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

#include "ppSub.h"

#define WCS_TOLERANCE 0.001             // Tolerance for WCS
//#define TESTING                         // For test output

#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
                     PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
                     PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources


// Copy every instance of a single keyword from one metadata to another
static bool metadataCopySingle(psMetadata *target, psMetadata *source, const char *name)
{
    PS_ASSERT_METADATA_NON_NULL(target, false);
    PS_ASSERT_METADATA_NON_NULL(source, false);
    PS_ASSERT_STRING_NON_EMPTY(name, false);

    psString regex = NULL;      // Regular expression
    psStringAppend(&regex, "^%s$", name);
    psMetadataIterator *iter = psMetadataIteratorAlloc(source, PS_LIST_HEAD, regex); // Iterator
    psFree(regex);
    psMetadataItem *item;       // Item from iteration
    while ((item = psMetadataGetAndIncrement(iter))) {
        psMetadataAddItem(target, item, PS_LIST_TAIL, PS_META_DUPLICATE_OK);
    }
    psFree(iter);

    return true;
}


bool ppSubReadout(pmConfig *config, psMetadata *stats, const pmFPAview *view)
{
    psTimerStart("PPSUB_MATCH");
    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
    pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
    pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources
    pmReadout *inConv = pmReadoutAlloc(NULL); // Convolved version of input
    pmReadout *refConv = pmReadoutAlloc(NULL); // Convolved version of reference
    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout: subtraction
    pmFPA *outFPA = outCell->parent->parent; // Output FPA
    pmHDU *outHDU = outFPA->hdu; // Output HDU
    if (!outHDU->header) {
        outHDU->header = psMetadataAlloc();
    }

    psImage *input = inRO->image;       // Input image
    psImage *reference = refRO->image;  // Reference image
    PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);

    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
    if (threads > 0) {
        if (!psThreadPoolInit(threads)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", threads);
            return false;
        }
        pmSubtractionThreadsInit(inRO, refRO);
    }

    // Look up recipe values
    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    psAssert(recipe, "We checked this earlier, so it should be here.");

    const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type
    psAssert(typeStr, "We put it here in ppSubArguments.c");
    pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel
    if (type == PM_SUBTRACTION_KERNEL_NONE) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", typeStr);
        psFree(inConv);
        psFree(refConv);
        psFree(outRO);
        return false;
    }

    bool mdok;                          // Status of MD lookup
    int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size
    int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order
    float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs
    float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
    int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size
    float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
    int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches
    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
    float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
    float sys = psMetadataLookupF32(NULL, recipe, "SYS"); // Relative systematic error in kernel
    bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction?
    psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
    psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders
    int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius
    int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order
    int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
    float penalty = psMetadataLookupF32(NULL, recipe, "PENALTY"); // Penalty for wideness
    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
    psImageMaskType 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
    psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
    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
    float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
    const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps

    bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
    float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
    float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
    float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
    float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
    int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
    bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution?

    // Statistics for renormalisation
    psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "RENORM.MEAN"));
    psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "RENORM.STDEV"));
    if (renormMean == PS_STAT_NONE || renormStdev == PS_STAT_NONE) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "Unable to parse renormalisation statistics from recipe.");
        return false;
    }
    int renormNum = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); // Number of samples for renormalisation
    float renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width of Gaussian phot

    psString interpModeStr = psMetadataLookupStr(&mdok, recipe, "INTERPOLATION"); // Interpolation mode
    psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp
    if (interpMode == PS_INTERPOLATE_NONE) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unknown interpolation mode: %s", interpModeStr);
        return false;
    }
    float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"

    pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode

    // Generate masks if they don't exist
    int numCols = input->numCols, numRows = input->numRows; // Image dimensions
    if (!inRO->mask) {
        if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) {
            pmReadoutSetMask(inRO, pmConfigMaskGet("SAT", config), pmConfigMaskGet("BAD", config));
        } else {
            inRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
            psImageInit(inRO->mask, 0);
        }
    }
    if (!refRO->mask) {
        if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) {
            pmReadoutSetMask(refRO, pmConfigMaskGet("SAT", config), pmConfigMaskGet("BAD", config));
        } else {
            refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
            psImageInit(refRO->mask, 0);
        }
    }

    if (!pmReadoutMaskNonfinite(inRO, pmConfigMaskGet("SAT", config))) {
        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in input.");
        return false;
    }
    if (!pmReadoutMaskNonfinite(refRO, pmConfigMaskGet("SAT", config))) {
        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference.");
        return false;
    }

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

    psArray *sources = NULL;            // Sources in image
    if (sourcesRO) {
        sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
    }

#if 0
    // Interpolate over bad pixels, so the bad pixels don't explode
    if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image.");
        return false;
    }
    if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image.");
        return false;
    }
    maskVal |= maskBad;
#endif

    // Match the PSFs
    if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, spacing, threshold,
                            sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
                            binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, sys,
                            maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
        psFree(inConv);
        psFree(refConv);
        psFree(outRO);
        return false;
    }
    psFree(optWidths);

    pmSubtractionThreadsFinalize(inRO, refRO);

    // Add kernel descrption to header
    pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis,
                                                        PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel
    if (!kernels) {
        kernels = psMetadataLookupPtr(&mdok, refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
    }
    if (!kernels) {
        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SUBTRACTION.KERNEL");
        psFree(inConv);
        psFree(refConv);
        psFree(outRO);
        return false;
    }
    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0,
                     "Subtraction kernel", kernels->description);

    {
        if (psMetadataLookup(inConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL)) {
            outRO->analysis = psMetadataCopy(outRO->analysis, inConv->analysis);
        } else if (psMetadataLookup(refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL)) {
            outRO->analysis = psMetadataCopy(outRO->analysis, refConv->analysis);
        } else {
            psWarning("Unable to find subtraction kernel in analysis metadata.");
        }

        // Update variance factors
        // It's not possible to do this perfectly, since we're combining different images:
        // vf_sum sigma_sum^2 = vf_1 * sigma_1^2 + vf_2 * sigma_2^2
        // It's not possible to write vf_sum as a function of vf_1 and vf_2 with no reference to the sigmas.
        // Instead, we're going to cheat.
        bool mdok;                      // Status of MD lookup
        pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, outRO->analysis,
                                                            PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
        if (!mdok) {
            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find subtraction kernels.");
            psFree(inConv);
            psFree(refConv);
            psFree(outRO);
            return false;
        }
        float vfIn = psMetadataLookupF32(NULL, inRO->parent->concepts,
                                         "CELL.VARFACTOR"); // Variance factor for input
        if (!isfinite(vfIn)) {
            vfIn = 1.0;
        }
        float vfRef = psMetadataLookupF32(NULL, refRO->parent->concepts,
                                          "CELL.VARFACTOR"); // Variance factor for ref
        if (!isfinite(vfRef)) {
            vfRef = 1.0;
        }
        if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
            vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
        }
        if (kernels->mode == PM_SUBTRACTION_MODE_2) {
            vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
        } else if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
            vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, true);
        }

        psStats *vfStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
        if (!psImageBackground(vfStats, NULL, inConv->weight, inConv->mask, maskVal | maskBad, rng)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved input");
            psFree(inConv);
            psFree(refConv);
            psFree(outRO);
            psFree(vfStats);
            psFree(rng);
            return false;
        }
        float inMeanVar = vfStats->robustMedian; // Mean variance of input
        if (!psImageBackground(vfStats, NULL, refConv->weight, refConv->mask, maskVal | maskBad, rng)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved reference");
            psFree(inConv);
            psFree(refConv);
            psFree(outRO);
            psFree(vfStats);
            psFree(rng);
            return false;
        }
        float refMeanVar = vfStats->robustMedian; // Mean variance of reference
        psFree(rng);
        psFree(vfStats);

        float vf = (vfIn * inMeanVar + vfRef * refMeanVar) / (inMeanVar + refMeanVar); // Resulting var factor
        psMetadataItem *vfItem = psMetadataLookup(outRO->parent->concepts, "CELL.VARFACTOR");
        vfItem->data.F32 = vf;

        // Statistics on the matching
        if (stats) {
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY);
            metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY);

            psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs",
                             psTimerClear("PPSUB_MATCH"));
        }
    }


#ifdef TESTING
    psImage *kernelImage = psMetadataLookupPtr(&mdok, inConv->analysis,
                                               "SUBTRACTION.KERNEL.IMAGE"); // Image of the kernels
    if (!kernelImage) {
        kernelImage = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.KERNEL.IMAGE");
    }
    psFits *fits = psFitsOpen("kernel.fits", "w");
    psFitsWriteImage(fits, NULL, kernelImage, 0, NULL);
    psFitsClose(fits);
#endif

    // Subtraction is: minuend - subtrahend
    pmReadout *minuend = inConv;
    pmReadout *subtrahend = refConv;

    if (reverse) {
        pmReadout *temp = subtrahend;
        subtrahend = minuend;
        minuend = temp;
    }

#ifdef TESTING
    {
        pmReadoutMaskApply(minuend, maskVal);
        psFits *fits = psFitsOpen("minuend.fits", "w");
        psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
        psFitsClose(fits);
    }
    {
        pmReadoutMaskApply(subtrahend, maskVal);
        psFits *fits = psFitsOpen("subtrahend.fits", "w");
        psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
            psFitsClose(fits);
    }
#endif

    // Photometry is to be performed in two stages:
    // 1. Measure the PSF using the PSF-matched images
    // 2. Find and measure sources on the subtracted image
    psTimerStart("PPSUB_PHOT");

    // Photometry stage 1: measure the PSF
    pmPSF *psf = NULL;                  // PSF for photometry
    if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
        outRO->image = psImageCopy(outRO->image, minuend->image, PS_TYPE_F32);
        if (minuend->weight) {
            outRO->weight = psImageCopy(outRO->weight, minuend->weight, PS_TYPE_F32);
        }
        if (minuend->mask) {
            outRO->mask = psImageCopy(outRO->mask, minuend->mask, PS_TYPE_IMAGE_MASK);
        }
        outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;

        if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
            psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
            if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth,
                                           renormMean, renormStdev, NULL)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
                psFree(outRO);
                return false;
            }
        }

        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
        pmFPACopy(photFile->fpa, outRO->parent->parent->parent);

        // We have a list of sources, so we could use those to redetermine the PSF model.
        // Until Gene makes the necessary adaptations to psphot, we will simply redetermine the PSF model from
        // scratch.

        if (psMetadataLookupBool(&mdok, recipe, "PSPHOT.VISUAL")) {
            psphotSetVisual(true);
        }

        // Need to ensure aperture residual is not calculated
        psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
        if (!psphotRecipe) {
            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find %s recipe.", PSPHOT_RECIPE);
            return false;
        }
        const char *breakpoint = psMetadataLookupStr(&mdok, psphotRecipe, "BREAK_POINT"); // Current break point
        psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE,
                         "ALTERED break point for psphot operations", "PSFMODEL");

        // set maskValue and markValue in the psphot recipe
        psImageMaskType maskValue = maskVal;
        psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking
        psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue);
        psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for marking",
                        markValue);

        if (!psphotReadout(config, view)) {
            psWarning("Unable to perform photometry on subtracted image.");
            psErrorStackPrint(stderr, "Error stack from photometry:");
            psErrorClear();
        }

        if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
            !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
            !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
            psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
            return false;
        }

        if (breakpoint) {
            psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE,
                             "RESTORED break point for psphot operations", breakpoint);
        }

        // Blow away the sources psphot found --- they're irrelevant for the subtraction
        pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with sources
        psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
        psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER");

        psf = psMemIncrRefCounter(psMetadataLookupPtr(NULL, photRO->parent->parent->analysis, "PSPHOT.PSF"));
        if (!psf) {
            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF from psphot");
            return false;
        }

        // Record the FWHM
        pmHDU *hdu = pmHDUFromCell(outRO->parent); // HDU with header
        psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MAJ");
        psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MIN");


        pmCell *photCell = pmFPAfileThisCell(config->files, view, "PSPHOT.INPUT");
        pmCellFreeReadouts(photCell);
    }

    // Do the actual subtraction
    outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
    if (inConv->weight && refConv->weight) {
        outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight);
    }
    outRO->mask = (psImage*)psBinaryOp(outRO->mask, inConv->mask, "|", refConv->mask);

    psFree(inConv);
    psFree(refConv);

    outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;

    // Copy astrometry over
    // It should get into the output images and photometry
    pmFPA *refFPA = refRO->parent->parent->parent; // Reference FPA
    pmHDU *refHDU = refFPA->hdu;        // Reference HDU
    if (!outHDU || !refHDU) {
        psWarning("Unable to find HDU at FPA level to copy astrometry.");
    } else if (!pmAstromReadWCS(outFPA, outCell->parent, refHDU->header, 1.0)) {
        psWarning("Unable to read WCS astrometry from reference FPA.");
        psErrorClear();
    } else if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
        psWarning("Unable to write WCS astrometry to output FPA.");
        psErrorClear();
    }

#if 0
    pmReadoutMaskApply(outRO, maskBad);
#endif

#ifdef TESTING
    for (int y = 0; y < outRO->image->numRows; y++) {
        for (int x = 0; x < outRO->image->numCols; x++) {
            if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
                printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x]);
            }
        }
    }
#endif

    if (!pmFPACopyConcepts(outCell->parent->parent, inRO->parent->parent->parent)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
        psFree(outRO);
        return false;
    }

#if 0
    // XXX Under development

    // QA: photometry of known sources with measured PSF
    if (sourcesRO && psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
        sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");



        // For each source, need to:
        // * generate appropriate model
        // * select pixels
        // * define fit radius
        // Don't need to measure sky (psphotFitSourcesLinear does that)

        psArray *dummy = psArrayAlloc(1); // Dummy array of sources, for psphotFitSourcesLinear with 1 source
        for (long i = 0; i < sources->n; i++) {
            pmSource *source = sources->data[i];
            if (!source || (source->mode & SOURCE_MASK)) {
                continue;
            }

            float origMag = source->psfMag; // Original magnitude

            // Coordinates of source
            float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
            float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];

            model = pmModelFromPSFforXY(psf, x, y, NAN);
            if (!fakeModel) {
                psWarning("Can't generate model");
                continue;
            }
            if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) {
                psErrorClear();
                continue;
            }
pmModel *pmModelFromPSFforXY (
    const pmPSF *psf,
    float Xo,
    float Yo,
    float Io
    );

bool pmSourceDefinePixels(
    pmSource *mySource,                 ///< source to be re-defined
    const pmReadout *readout,  ///< base the source on this readout
    psF32 x,                            ///< center coords of source
    psF32 y,                            ///< center coords of source
    psF32 Radius                        ///< size of box on source
);


bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {

    }
#endif

    // Photometry stage 2: find and measure sources on the subtracted image
    if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
        // The PSF should already be stored for the readout
        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
        pmFPACopy(photFile->fpa, outFPA);

        pmReadout *psfRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.PSF.LOAD");
        if (!psfRO) {
            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find readout for file PSPHOT.PSF.LOAD");
            return false;
        }
        psMetadataAddPtr(psfRO->parent->parent->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
                         "PSF from matched addition", psf);
        psFree(psf);

        if (psMetadataLookupBool(&mdok, recipe, "PSPHOT.VISUAL")) {
            psphotSetVisual(true);
        }

        if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
            psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
            if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth,
                                           renormMean, renormStdev, NULL)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
                psFree(outRO);
                return false;
            }
        }

        // Need to ensure aperture residual is not calculated
        psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
        psMetadataItem *item = psMetadataLookup(psphotRecipe, "MEASURE.APTREND"); // Item determining aptrend
        if (!item) {
            psWarning("Unable to find MEASURE.APTREND in psphot recipe");
            psErrorClear();
        } else {
            item->data.B = false;
        }

        if (!psphotReadout(config, view)) {
            psWarning("Unable to perform photometry on subtracted image.");
            psErrorStackPrint(stderr, "Error stack from photometry:");
            psErrorClear();
        }

        if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
            !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
            !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
            psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
            return false;
        }

        pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
        pmFPAfileActivate(config->files, false, "PSPHOT.LOAD.PSF");

        if (stats) {
            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 ? sources->n : 0);
            psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0, "Time to do photometry",
                             psTimerClear("PPSUB_PHOT"));
        }

#ifdef TESTING
        // Record data about sources: not everything gets into the output CMF files
        {
            pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
            psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
            FILE *sourceFile = fopen("sources.dat", "w"); // File for sources
            fprintf(sourceFile,
                    "# x y mag mag_err psf_chisq cr_nsigma ext_nsigma psf_qf flags m_x m_y m_xx m_xy m_yy\n");
            for (int i = 0; i < sources->n; i++) {
                pmSource *source = sources->data[i];
                if (!source) {
                    continue;
                }

                float x, y;             // Position of source
                float chi2;             // chi^2 for source
                if (source->modelPSF) {
                    x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
                    y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
                    chi2 = source->modelPSF->chisq;
                } else if (source->peak) {
                    x = source->peak->xf;
                    y = source->peak->yf;
                    chi2 = NAN;
                } else {
                    psWarning("No position available for source.");
                    continue;
                }

                float xMoment = NAN, yMoment = NAN, xxMoment = NAN, xyMoment = NAN, yyMoment = NAN;
                if (source->moments) {
                    xMoment = source->moments->Mx;
                    yMoment = source->moments->My;
                    xxMoment = source->moments->Mxx;
                    xyMoment = source->moments->Mxy;
                    yyMoment = source->moments->Myy;
                }

                fprintf(sourceFile, "%f %f %f %f %f %f %f %f %d %f %f %f %f %f\n",
                        x, y, source->psfMag, source->errMag, chi2, source->crNsigma, source->extNsigma,
                        source->pixWeight, source->mode, xMoment, yMoment, xxMoment, xyMoment, yyMoment);
            }
            fclose(sourceFile);
        }
#endif

    }

    // Higher order background subtraction using psphot
    if (!ppSubBackground(outRO, config, view)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
        psFree(outRO);
        return false;
    }

    // Renormalising for pixels, because that's what magic desires
    if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
        psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
        if (!pmReadoutWeightRenormPixels(outRO, maskValue, renormMean, renormStdev, NULL)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
            psFree(outRO);
            return false;
        }
    }

    // Add additional data to the header
    pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF"); // Reference file
    pmFPAfile *inFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT"); // Input file
    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.REFERENCE", 0,
                     "Subtraction reference", refFile->filename);
    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.INPUT", 0,
                     "Subtraction input", inFile->filename);


    // 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, "PPSUB.OUTPUT.JPEG1");
        pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG2");

        pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
        if (!pmReadoutRebin(ro1, outRO, maskBad, 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);
    }

#ifdef TESTING
    // Significance image
    {
        psImage *sig = (psImage*)psBinaryOp(NULL, outRO->image, "*", outRO->image);
        psBinaryOp(sig, sig, "/", outRO->weight);
        psFits *fits = psFitsOpen("significance.fits", "w");
        psFitsWriteImage(fits, NULL, sig, 0, NULL);
        psFitsClose(fits);
        psFree(sig);
    }
#endif

    psFree(outRO);

    return true;
}
