#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

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
    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
    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.
    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

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

    // Add kernel descrption to header
    pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, outRO->analysis,
                                                        "SUBTRACTION.KERNEL"); // The subtraction kernels
    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0,
                     "Subtraction kernel", kernels->description);


    psImage *kernelImage = psMetadataLookupPtr(NULL, outRO->analysis,
                                               "SUBTRACTION.KERNEL.IMAGE"); // Image of the kernels
    psFits *fits = psFitsOpen("kernel.fits", "w");
    psFitsWriteImage(fits, NULL, kernelImage, 0, NULL);
    psFitsClose(fits);

    // Do the subtraction
    if (reverse) {
        psBinaryOp(outRO->image, outRO->image, "-", inRO->image);
    } else {
        psBinaryOp(outRO->image, inRO->image, "-", outRO->image);
    }
    psBinaryOp(outRO->mask, outRO->mask, "|", inRO->mask);
    if (outRO->weight) {
        psBinaryOp(outRO->weight, outRO->weight, "+", inRO->weight);
    }

#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 (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) {
        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
        pmFPACopy(photFile->fpa, outRO->parent->parent->parent);

        if (!psphotReadout(config, view)) {
            psError(psErrorCodeLast(), false, "Unable to perform photometry on subtracted image.\n");
            return false;
        }

        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)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from reference FPA.");
            psFree(outRO);
            return false;
        }
        if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA.");
            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);

    psFree(outRO);

    return true;
}
