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

#include <stdio.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


bool ppSubReadout(pmConfig *config, const pmFPAview *view)
{
    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);

    // Look up appropriate values
    bool mdok;                          // Status of MD lookup
    int size = psMetadataLookupS32(NULL, config->arguments, "KERNEL.SIZE"); // Kernel half-size
    int order = psMetadataLookupS32(NULL, config->arguments, "SPATIAL.ORDER"); // Spatial polynomial order
    float regionSize = psMetadataLookupF32(NULL, config->arguments, "REGION.SIZE"); // Size of iso-kernel regs
    float spacing = psMetadataLookupF32(NULL, config->arguments, "STAMP.SPACING"); // Typical stamp spacing
    int footprint = psMetadataLookupS32(NULL, config->arguments, "STAMP.FOOTPRINT"); // Stamp half-size
    float threshold = psMetadataLookupF32(NULL, config->arguments, "STAMP.THRESHOLD"); // Threshold for stmps
    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
    float rej = psMetadataLookupF32(NULL, config->arguments, "REJ"); // Rejection threshold
    pmSubtractionKernelsType type = psMetadataLookupS32(NULL, config->arguments,
                                                        "KERNEL.TYPE"); // Kernel type
    bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction?
    psVector *widths = psMetadataLookupPtr(NULL, config->arguments, "ISIS.WIDTHS"); // ISIS Gaussian widths
    psVector *orders = psMetadataLookupPtr(NULL, config->arguments, "ISIS.ORDERS"); // ISIS Polynomial orders
    int inner = psMetadataLookupS32(NULL, config->arguments, "INNER"); // Inner radius
    int ringsOrder = psMetadataLookupS32(NULL, config->arguments, "RINGS.ORDER"); // RINGS polynomial order
    int binning = psMetadataLookupS32(NULL, config->arguments, "SPAM.BINNING"); // Binning for SPAM kernel
    psMaskType maskBad = pmConfigMask(psMetadataLookupStr(NULL, config->arguments, "MASK.BAD"),
                                    config); // Value to mask
    psMaskType maskBlank = pmConfigMask(psMetadataLookupStr(NULL, config->arguments, "MASK.BLANK"),
                                        config); // Mask for blank reg.
    float badFrac = psMetadataLookupF32(NULL, config->arguments, "BADFRAC"); // Maximum bad fraction
    const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps

    bool optimum = psMetadataLookupBool(&mdok, config->arguments, "OPTIMUM"); // Derive optimum parameters?
    float optMin = psMetadataLookupF32(&mdok, config->arguments, "OPTIMUM.MIN"); // Minimum width for search
    float optMax = psMetadataLookupF32(&mdok, config->arguments, "OPTIMUM.MAX"); // Maximum width for search
    float optStep = psMetadataLookupF32(&mdok, config->arguments, "OPTIMUM.STEP"); // Step for search
    float optThresh = psMetadataLookupF32(&mdok, config->arguments, "OPTIMUM.TOL"); // Tolerance for search
    int optOrder = psMetadataLookupS32(&mdok, config->arguments, "OPTIMUM.ORDER"); // Order for search
    bool dual = psMetadataLookupBool(&mdok, config->arguments, "DUAL"); // Dual convolution?
    bool renorm = psMetadataLookupBool(&mdok, config->arguments, "RENORM"); // Renormalise weights?
    int renormWidth = psMetadataLookupS32(&mdok, config->arguments, "RENORM.WIDTH"); // Width for renormalise

    pmSubtractionMode mode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtraction mode

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

    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 (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, regionSize, spacing, threshold, sources,
                            stampsName, type, size, order, widths, orders, inner, ringsOrder,
                            binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
                            maskBlank, badFrac, mode)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
        psFree(inConv);
        psFree(refConv);
        psFree(outRO);
        return false;
    }
    psFree(optWidths);

    // Add kernel descrption to header
    pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis,
                                                        "SUBTRACTION.KERNEL"); // The subtraction kernels
    if (!kernels) {
        kernels = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.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);

#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

    // Do the subtraction
    {
        // Subtraction is: minuend - subtrahend
        pmReadout *minuend = inConv;
        pmReadout *subtrahend = refConv;

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

        outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
        outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
        if (minuend->weight && subtrahend->weight) {
            outRO->weight = (psImage*)psBinaryOp(outRO->weight, minuend->weight, "+", subtrahend->weight);
        }
        outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;

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

    psFree(inConv);
    psFree(refConv);

#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.U8[y][x] & maskBlank)) {
                printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.U8[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 (renorm) {
        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
        if (!pmReadoutWeightRenorm(outRO, maskBlank, PS_STAT_ROBUST_MEDIAN, PS_STAT_ROBUST_STDEV,
                                   renormWidth, rng)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to renormalise weights for photometry.");
            psFree(outRO);
            return false;
        }
        psFree(rng);
    }

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

        // Need to ensure aperture residual is not calculated
        psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
        psMetadataItem *item = psMetadataLookup(recipe, "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();
        }

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

    // Copy astrometry over
    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();
        }
    }

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

    psFree(outRO);

    return true;
}
