# include "pswarp.h"
#include <ppStats.h>

// Loop over the inputs, warp them to the output skycell and then write out the output.
bool pswarpLoop(pmConfig *config)
{

    psMaskType maskPoor = psMetadataLookupU8(NULL, config->arguments, "MASK.POOR"); // Mask for "poor" data
    psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Mask for bad data

    // select the input data sources
    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT");
    if (!input) {
        psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
        return false;
    }

    // use the external astrometry source if supplied
    pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
    if (!astrom) {
        astrom = input;
    }

    // select the output readout
    pmFPAview *view = pmFPAviewAlloc(0);
    view->chip = 0;
    view->cell = 0;
    view->readout = 0;
    pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT");
    if (!output) {
        psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n");
        return false;
    }
    psFree (view);

    bool mdok;                          // Status of MD lookup
    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 (mdok && statsName && strlen(statsName) > 0) {
        psString resolved = pmConfigConvertFilename(statsName, config, true);
        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;
        }
        psFree(resolved);
        stats = psMetadataAlloc();
    }

    // Read the input astrometry
    {
        pmFPAfileActivate(config->files, false, "PSWARP.INPUT");
        pmFPAfileActivate(config->files, false, "PSWARP.MASK");
        pmFPAfileActivate(config->files, false, "PSWARP.WEIGHT");
        pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
        pmFPAfileActivate(config->files, false, "PSWARP.OUTPUT");
        pmFPAfileActivate(config->files, false, "PSWARP.OUTPUT.MASK");
        pmFPAfileActivate(config->files, false, "PSWARP.OUTPUT.WEIGHT");
        pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);
        pmChip *chip;
        pmFPAview *view = pmFPAviewAlloc(0);
        while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
            psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
            if (!chip->process || !chip->file_exists) { continue; }
            pmFPAfileIOChecks(config, view, PM_FPA_BEFORE);
            pmCell *cell;
            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
                psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
                if (!cell->process || !cell->file_exists) { continue; }
                pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);
            }
        }
        psFree(view);

        pmFPAfileActivate(config->files, true, "PSWARP.INPUT");
        pmFPAfileActivate(config->files, true, "PSWARP.MASK");
        pmFPAfileActivate(config->files, true, "PSWARP.WEIGHT");
        pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
        pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT");
        pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.MASK");
        pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.WEIGHT");
    }

    pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL");
    view = pmFPAviewAlloc(0);

    // find the FPA phu
    bool bilevelAstrometry = false;
    pmHDU *phu = pmFPAviewThisPHU (view, astrom->fpa);
    if (phu) {
        char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
        if (ctype) {
            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
        }
    }
    if (bilevelAstrometry) {
        if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA.");
            psFree(view);
            psFree(stats);
            return false;
        }
    }

    psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging

    // files associated with the science image
    pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);

    pmChip *chip;
    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
        psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
        if (!chip->process || !chip->file_exists) { continue; }
        pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);

        // read WCS data from the corresponding header
        pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
        if (bilevelAstrometry) {
            if (!pmAstromReadBilevelChip (chip, hdu->header)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel chip astrometry for input FPA.");
                psFree(view);
                psFree(stats);
                return false;
            }
        } else {
            // we use a default FPA pixel scale of 1.0
            if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA.");
                psFree(view);
                psFree(stats);
                return false;
            }
        }

        pmCell *cell;
        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
            psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
            if (!cell->process || !cell->file_exists) { continue; }
            pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);

            psListAdd(cells, PS_LIST_TAIL, cell);

            // process each of the readouts
            pmReadout *readout;
            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
                pmFPAfileIOChecks(config, view, PM_FPA_BEFORE);
                if (!readout->data_exists) {
                    continue;
                }

                pswarpTransformReadout_Opt(output, readout, config);

                pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
            }
            pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
        }
        pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
    }

    pmCell *outCell = output->parent;   // Output cell
    pmChip *outChip = outCell->parent;  // Output chip
    pmFPA *outFPA = outChip->parent;    // Output FPA

    if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts.");
        psFree(stats);
        psFree(view);
        return false;
    }

    if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA concepts from input to output.");
        psFree(stats);
        psFree(view);
        return false;
    }

    output->data_exists = true;
    outCell->data_exists = true;
    outChip->data_exists = true;

#if 0
    float gain = 0.0, readnoise = 0.0, saturation = INFINITY, exposure = 0.0, darktime = 0.0;
    double mjd = 0.0;
#endif

#if 0
    // XXX Hack so that INSTRUME header keyword and the like are not updated properly
    // If INSTRUME (etc) is updated, then the FITS file with the warped image will be often
    // recognised as coming from the old camera, not the new, warped camera.
    // The proper solution to this, I believe, would be to define a skycell derivative of the
    // old camera, so that the skycell can be recognised as coming from the old camera (and would
    // also therefore inherit all its recipes, which would be good), rather than having a completely
    // different camera.  That will take some work; the below solution should be good enough for now.
    {
        psMetadataItem *item;           // Concepts item to update

        item = psMetadataLookup(outFPA->concepts, "FPA.INSTRUMENT");
        psFree(item->data.V);
        item->data.V = NULL;

        item = psMetadataLookup(outFPA->concepts, "FPA.DETECTOR");
        psFree(item->data.V);
        item->data.V = NULL;
    }
#endif


    // Perform statistics on the output image
    if (stats) {
        if (!ppStats(stats, output->parent->parent->parent, view, maskBad | maskPoor, config)) {
            psWarning("Unable to perform statistics on warped image.");
        }
    }

    pmHDU *hdu = output->parent->parent->parent->hdu; // HDU for the output warped image
    if (!hdu->header) {
        hdu->header = psMetadataAlloc();
    }

    // Add MD5 information for readout
    const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
    const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
    psString headerName = NULL; // Header name for MD5
    psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
    psVector *md5 = psImageMD5(output->image); // md5 hash
    psString md5string = psMD5toString(md5); // String
    psFree(md5);
    psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
                     "Image MD5", md5string);
    psFree(md5string);
    psFree(headerName);

    pmFPAfileIOChecks (config, view, PM_FPA_AFTER);
    psFree (view);

    // Write out summary statistics
    if (stats) {
        const char *statsMDC = psMetadataConfigFormat(stats);
        if (!statsMDC) {
            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
            psFree(stats);
            fclose(statsFile);
            return false;
        }
        fprintf(statsFile, "%s", statsMDC);
        psFree((void*)statsMDC);
        fclose(statsFile);

        psFree(stats);
    }

    return true;
}
