#include "ppbgrestore.h"
#include "psphot.h"

#define ESCAPE(MESSAGE) { \
  psError(PS_ERR_UNKNOWN, false, MESSAGE); \
  psFree(view); \
  return false; \
}


bool ppbgrestoreLoop(pmConfig *config, ppbgrestoreOptions *options)
{
    bool status;                        // Status of MD lookup
    pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PPBGRESTORE.INPUT");
    if (!status) {
        psErrorStackPrint(stderr, "Can't find input data!\n");
        ppbgrestoreCleanup(config, options);
        exit(PS_EXIT_PROG_ERROR);
    }
    pmFPAfile *backmdl = psMetadataLookupPtr(&status, config->files, "PPBGRESTORE.BACKMDL");
    if (!status) {
        psErrorStackPrint(stderr, "Can't find background model data!\n");
        ppbgrestoreCleanup(config, options);
        exit(PS_EXIT_PROG_ERROR);
    }
    pmFPAfile *output= psMetadataLookupPtr(&status, config->files, "PPBGRESTORE.OUTPUT");
    if (!status) {
        psErrorStackPrint(stderr, "Can't find output data!\n");
        ppbgrestoreCleanup(config, options);
        exit(PS_EXIT_PROG_ERROR);
    }
    pmFPAfile *background= psMetadataLookupPtr(&status, config->files, "PPBGRESTORE.BACKGROUND");
    if (!status) {
        psErrorStackPrint(stderr, "Can't find output data!\n");
        ppbgrestoreCleanup(config, options);
        exit(PS_EXIT_PROG_ERROR);
    }
    bool save_background = psMetadataLookupBool(&status, config->arguments, "SAVE_BACKGROUND");

    pmFPAview *view = pmFPAviewAlloc(0);// View for level of interest
//    pmHDU *lastHDU = NULL;              // Last HDU that was updated

    // files associated with the science image
    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
        ESCAPE("load failure for FPA");
    }

    pmChip *chip;                       // Chip from FPA
    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
        psLogMsg ("ppbgrestoreLoop", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
        if (!chip->process || !chip->file_exists) {
            continue;
        }
        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
            ESCAPE("load failure for Chip");
        }

        // much of the following is adapted from ppImageSubtractBackground

        if (chip->cells->n == 0) {
            psError(PS_ERR_UNKNOWN, true, "Chip has no cells");
            return false;
        }

        if (chip->cells->n > 1) {
            psError(PS_ERR_UNKNOWN, true, "Chip has too many cells: %ld", chip->cells->n);
            return false;
        }
        pmCell *cell = chip->cells->data[0]; // Cell of interest
        if (!cell || !cell->process || !cell->file_exists) {
            // Nothing to process
            return false;
        }
        if (cell->readouts->n == 0) {
            psError(PS_ERR_UNKNOWN, true, "Cell has no readouts");
            return false;
        }
        if (cell->readouts->n > 1) {
            psWarning("Cell has %ld readouts; only the first will be processed", cell->readouts->n);
        }
        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
        if (!ro || !ro->data_exists) {
            psError(PS_ERR_UNKNOWN, true, "readout has no data");
            return false;
        }
        psImage *image = ro->image;
//        psImage *mask = NULL;
//        psImageMaskType maskVal = 0;

        pmFPAview roView = *view;
        roView.cell = 0;
        roView.readout = 0;
        pmReadout *modelRO = pmFPAviewThisReadout(&roView, backmdl->fpa);
        if (!modelRO || !modelRO->data_exists) {
            psError(PS_ERR_UNKNOWN, true, "model readout has no data");
            return false;
        }
        output->save = true;
        background->save = save_background;

        psImageBinning *binning = psphotBackgroundBinning(image, config);
        if (!binning) {
            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find background binning");
            return false;
        }
        pmReadout *backgroundRO;
#define USE_CELL
#ifdef USE_CELL
        pmCell *backgroundCell = pmFPAviewThisCell(&roView, background->fpa);
        backgroundRO = pmReadoutAlloc(backgroundCell);
#else
        backgroundRO = pmFPAviewThisReadout(&roView, background->fpa);
//        backgroundRO = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKGND", image->numCols, image->numRows, PS_TYPE_F32);
#endif
        backgroundRO->image = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
	backgroundRO->data_exists = true;
	backgroundRO->parent->data_exists = true;
	backgroundRO->parent->parent->data_exists = true;

        // linear interpolation to full-scale
        if (!psImageUnbin (backgroundRO->image, modelRO->image, binning)) {
            psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
            psFree(binning);
            return false;
        }
        psFree(binning);
        psF32 **backData = backgroundRO->image->data.F32;

        // Undo the background subtraction
        int numCols = image->numCols, numRows = image->numRows; // Size of image
        long nancount = 0;
        long naninputs = 0;
        long finitecount = 0;
        for (int y = 0; y < numRows; y++) {
            for (int x = 0; x < numCols; x++) {
                float imageval =  image->data.F32[y][x] ;
                float bgvalue = backData[y][x];
                if (!isfinite(bgvalue)) {
                    // if background is NAN leave the image unchanged
                    // According to ppImageSubtractBackground pixels where the background is not finite should be NAN
                    // experimentally
                    // for breakpoint
                    nancount++;
                    if (isfinite(imageval)) {
                         finitecount++;
                         //   psWarning("unexpected finite image value %f found at %d %d background is infinite", imageval, x, y);
                    }
                    //      Doing this is causes a lot of new NAN pixels to appear.
                    //      NANS are appearing in the model for some reason.
                    //      image->data.F32[y][x] = NAN;
                } else if (isfinite(imageval)) {
                    // XXX: TODO:following ppImage we should not modify pixels that are masked.
                    image->data.F32[y][x] = imageval + bgvalue;
                } else {
		    naninputs++;
		}
            }
        }
        if (naninputs) {
            printf("%ld infinite pixels found in input image\n", naninputs);
        }
        if (nancount) {
            psWarning("%ld infinite pixels found in background %ld pixels were finite\n", nancount, finitecount);
        }
//        pmFPAfileDropInternal(config->files, "PSPHOT.BACKGND");

        // Close chip
        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
            ESCAPE("save failure for Chip");
        }
    }

    if (psTraceGetLevel("ppbgrestore") >= 3) {
//        ppbgrestoreFileCheck(config);
    }

    // Output and Close FPA
    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
        ESCAPE("save failure for FPA");
    }
    psFree(view);

    return true;
}
