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

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

#include "ppSkycell.h"

#define BUFFER 16                       // Size of buffer for projections

static void regionMinMax(psRegion *base,// Base region; modified
                         const psRegion *compare // Comparison region
                         )
{
    base->x0 = PS_MIN(base->x0, compare->x0);
    base->x1 = PS_MAX(base->x1, compare->x1);
    base->y0 = PS_MIN(base->y0, compare->y0);
    base->y1 = PS_MAX(base->y1, compare->y1);
    return;
}

static psRegion *skycellRegion(pmAstromWCS *wcs, // World Coordinate System
                               int numCols, int numRows // Size of image
    )
{
    psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms
    psRegion *region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY); // Region for skycell

// Limit the region using the nominated point (X,Y)
#define REGION_RANGE(X,Y) \
    fromCoords->x = (X); \
    fromCoords->y = (Y); \
    psPlaneTransformApply(toCoords, wcs->trans, fromCoords); \
    toCoords->x /= wcs->cdelt1; \
    toCoords->y /= wcs->cdelt2; \
    toCoords->x += wcs->crpix1; \
    toCoords->y += wcs->crpix2; \
    region->x0 = PS_MIN(region->x0, toCoords->x); \
    region->x1 = PS_MAX(region->x1, toCoords->x); \
    region->y0 = PS_MIN(region->y0, toCoords->y); \
    region->y1 = PS_MAX(region->y1, toCoords->y);

    REGION_RANGE(0, 0);                 // Lower left
    psTrace("ppSkycell", 9, "Start: %.0f %.0f", toCoords->x, toCoords->y);
    REGION_RANGE(0, numRows);           // Upper left
    REGION_RANGE(numCols, 0);           // Lower right
    REGION_RANGE(numCols, numRows);     // Upper right
    psTrace("ppSkycell", 9, "Stop: %.0f %.0f\n", toCoords->x, toCoords->y);

    return region;
}

// Activate/deactivate a single element for a list
void fileActivationSingle(pmConfig *config, // Configuration
                          const char **files, // Files to turn on/off
                          bool state,   // Activation state
                          int num // Number of file in sequence
                          )
{
    assert(config);
    for (int i = 0; files[i] != NULL; i++) {
        pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file
    }
    return;
}

// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
static pmFPAview *filesIterateDown(pmConfig *config // Configuration
                                   )
{
    assert(config);

    pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
        return NULL;
    }
    view->chip = 0;
    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
        return NULL;
    }
    view->cell = 0;
    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
        return NULL;
    }
    view->readout = 0;
    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
        return NULL;
    }
    return view;
}

// Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells
static bool filesIterateUp(pmConfig *config // Configuration
                           )
{
    assert(config);

    pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
    view->chip = view->cell = view->readout = 0;
    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
        return false;
    }
    view->readout = -1;
    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
        return false;
    }
    view->cell = -1;
    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
        return false;
    }
    view->chip = -1;
    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
        return false;
    }
    psFree(view);
    return true;
}


bool ppSkycellLoop(ppSkycellData *data // Run-time data
    )
{
    psVector *crval1 = psVectorAllocEmpty(BUFFER, PS_TYPE_F64); // CRVAL1 values
    psVector *crval2 = psVectorAllocEmpty(BUFFER, PS_TYPE_F64); // CRVAL2 values
    psVector *cdelt1 = psVectorAllocEmpty(BUFFER, PS_TYPE_F64); // CDELT1 values
    psVector *cdelt2 = psVectorAllocEmpty(BUFFER, PS_TYPE_F64); // CDELT2 values
    psArray *projRegions = psArrayAllocEmpty(BUFFER); // Region for projection
    int numProj = 0;                    // Number of projections

    psVector *target = psVectorAlloc(data->numInputs, PS_TYPE_S32); // Target for each input
    psArray *imageRegions = psArrayAlloc(data->numInputs); // Region for image

    for (int i = 0; i < data->numInputs; i++) {
        pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine
        // Header in the FPA should have been read as a part of defining the file...
#if 0
        bool pmFPAReadHeaderSet(pmFPA *fpa,     // FPA to read into
                                psFits *fits,   // FITS file from which to read
                                pmConfig *config // Configuration
            );
#endif
        pmHDU *hdu = file->fpa->hdu;    // Header of interest

        int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns
        int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows

        pmAstromWCS *wcs = pmAstromWCSfromHeader(hdu->header); // World Coordinate System
        if (!wcs) {
            psError(psErrorCodeLast(), false, "Unable to read WCS for image %d", i);
            return false;
        }

        psRegion *region = imageRegions->data[i] = skycellRegion(wcs, numCols, numRows); // Region of image
        psTrace("ppSkycell", 5, "Image region %d is: [%.0f:%.0f,%.0f:%.0f]\n",
                i, region->x0, region->x1, region->y0, region->y1);

        bool found = false;             // Found a projection?
        for (int j = 0; j < numProj && !found; j++) {
            if (wcs->crval1 == crval1->data.F64[j] && wcs->crval2 == crval2->data.F64[j] &&
                wcs->cdelt1 == cdelt1->data.F64[j] && wcs->cdelt2 == cdelt1->data.F64[j]) {
                regionMinMax(projRegions->data[j], region);
                target->data.S32[i] = j;
                found = true;
                psTrace("ppSkycell", 3, "Image %d uses projection %d\n", i, j);
            }
        }

        if (!found) {
            psVectorAppend(crval1, wcs->crval1);
            psVectorAppend(crval2, wcs->crval2);
            psVectorAppend(cdelt1, wcs->cdelt1);
            psVectorAppend(cdelt2, wcs->cdelt2);
            psRegion *projRegion = psRegionAlloc(region->x0, region->x1, region->y0, region->y1);
            psArrayAdd(projRegions, projRegions->n, projRegion);
            psFree(projRegion);
            target->data.S32[i] = numProj;
            psTrace("ppSkycell", 3, "Image %d uses new projection\n", i);

            numProj++;
        }
    }

    pmFPAfileActivate(data->config->files, false, NULL);

    for (int i = 0; i < numProj; i++) {
        psRegion *projRegion = projRegions->data[i]; // Region for skycell projection
        psTrace("ppSkycell", 2, "Projection %d: [%.0f:%.0f,%.0f:%.0f]\n",
                i, projRegion->x0, projRegion->x1, projRegion->y0, projRegion->y1);
        int xSize = projRegion->x1 - projRegion->x0 + 1; // Size of unbinned image
        int ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image
        // Size of binned image 1
        int numCols1 = xSize / (float)data->bin1 + 1.5, numRows1 = ySize / (float)data->bin1 + 1.5;
        // Size of binned image 2
        int numCols2 = numCols1 / (float)data->bin2 + 1.5, numRows2 = numRows1 / (float)data->bin2 + 1.5;

        psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image
        psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image
        psImageInit(image1, 0);
        psImageInit(image2, 0);

        psImage *mask1 = NULL, *mask2 = NULL; // Binned masks
        if (data->masksName) {
            mask1 = psImageAlloc(numCols1, numRows1, PS_TYPE_IMAGE_MASK);
            mask2 = psImageAlloc(numCols2, numRows2, PS_TYPE_IMAGE_MASK);
            psImageInit(mask1, 0xFF);
            psImageInit(mask2, 0xFF);
        }

        for (int j = 0; j < data->numInputs; j++) {
            if (target->data.S32[j] != i) {
                continue;
            }

            pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.IMAGE", j);
            if (data->masksName) {
                pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.MASK", j);
            }

            pmFPAview *view = filesIterateDown(data->config); // View to readout
            if (!view) {
                psError(psErrorCodeLast(), false, "Unable to iterate down.");
                // XXX Cleanup
                return false;
            }

            pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j);
            pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input
            psFree(view);

            // Flip images; no idea why this has to be done, but apparently it does
            {
                psImage *rot = psImageRotate(NULL, inRO->image, M_PI, NAN, PS_INTERPOLATE_BILINEAR);
                psFree(inRO->image);
                inRO->image = rot;
            }
            if (inRO->mask) {
                psImage *rot = psImageRotate(NULL, inRO->mask, M_PI, 0, PS_INTERPOLATE_FLAT);
                psFree(inRO->mask);
                inRO->mask = rot;
            }

            pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts
            if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) {
                psError(psErrorCodeLast(), false, "Unable to rebin image");
                // XXX Cleanup
                return false;
            }
            if (!pmReadoutRebin(bin2RO, bin1RO, data->maskVal, data->bin2, data->bin2)) {
                psError(psErrorCodeLast(), false, "Unable to rebin image");
                // XXX Cleanup
                return false;
            }

            psRegion *imageRegion = imageRegions->data[j]; // Region for image
            // Offsets for image on skycell
            int xOffset1 = (imageRegion->x0 - projRegion->x0) / (float)data->bin1;
            int yOffset1 = (imageRegion->y0 - projRegion->y0) / (float)data->bin1;
            int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2;

            // XXX Completely neglecting rotations
            // The skycells are divided up neatly with them all having the same orientation
            psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "=");
            psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "=");
            if (data->masksName) {
                psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "=");
                psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "=");
            }

            psFree(bin1RO);
            psFree(bin2RO);
            filesIterateUp(data->config);
            psFree(file->fpa);
            file->fpa = NULL;
            if (data->masksName) {
                pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.MASK", j);
                psFree(file->fpa);
                file->fpa = NULL;
            }
            pmFPAfileActivate(data->config->files, false, NULL);
        }

        pmFPAfileActivate(data->config->files, true, "PPSKYCELL.JPEG1");
        pmFPAfileActivate(data->config->files, true, "PPSKYCELL.JPEG2");
        pmFPAview *view = filesIterateDown(data->config); // View to readout

        pmCell *cell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.JPEG1"); // Rebinned cell 1
        pmCell *cell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.JPEG2"); // Rebinned cell 2
        psFree(view);
        pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts

        ro1->image = image1;
        ro2->image = image2;
        ro1->mask = mask1;
        ro2->mask = mask2;

        ro1->data_exists = cell1->data_exists = cell1->parent->data_exists = true;
        ro2->data_exists = cell2->data_exists = cell2->parent->data_exists = true;

        pmFPAfile *file1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.JPEG1", 0);
        file1->save = true;
        file1->index = i;
        pmFPAfile *file2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.JPEG2", 0);
        file2->save = true;
        file2->index = i;
#if 0
        {
            psString filename = NULL;   // Filename for image
            psStringAppend(&filename, "skycell_%d_1.fits", i);
            psFits *fits = psFitsOpen(filename, "w");
            psFree(filename);
            psFitsWriteImage(fits, NULL, image1, 0, NULL);
            psFitsClose(fits);
        }

        {
            psString filename = NULL;   // Filename for image
            psStringAppend(&filename, "skycell_%d_2.fits", i);
            psFits *fits = psFitsOpen(filename, "w");
            psFree(filename);
            psFitsWriteImage(fits, NULL, image2, 0, NULL);
            psFitsClose(fits);
        }
#endif

        filesIterateUp(data->config);
    }

    // XXX Cleanup
    return true;
}
