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

#include "ppstamp.h"
#include "pmAstrometryUtils.h"

#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS

typedef enum {
    PPSTAMP_OFF,
    PPSTAMP_PARTIALLY_ON,
    PPSTAMP_ON,
    PSTAMP_ERROR
} ppstampOverlap;

static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);

// convert the input chip's transforms to the output
static bool convertWCS(pmHDU *outHDU, pmFPA *outFPA, pmChip *outChip, pmChip *inChip, psRegion *roi)
{
    pmHDU *hdu = pmHDUFromChip(inChip);
    PS_ASSERT_PTR_NON_NULL(hdu, 1)
    PS_ASSERT_PTR_NON_NULL(hdu->header, 1)

    // Change the reference pixel to account for the change in origin between the stamp and
    // the original image
    outChip->toFPA   = psPlaneTransformSetCenter(NULL, inChip->toFPA, (int) roi->x0, (int) roi->y0);

    outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);

    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_NONLIN_TOL)) {
        psError(PS_ERR_UNKNOWN, false, "Failed to write WCS to output\n");
        return false;
    }

    return true;
}

static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options)
{
    pmChip    *outChip;
    pmFPAview *view = pmFPAviewAlloc(0);

    // our output file has a single chip
    view->chip = 0;
    outChip = pmFPAviewThisChip(view, output->fpa);
    psFree(view);

    // copy data from the input chip header to the output.
    // since some of the keywords might be duplicated we may not want to copy both

    // copy the fpa concepts
    outChip->parent->concepts = psMetadataCopy(outChip->parent->concepts, inChip->parent->concepts);

    pmHDU *inHDU  = pmHDUFromChip(inChip);
    pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);

    if (inHDU->header) {
        outHDU->header = psMetadataCopy(outHDU->header, inHDU->header);
    } else if (!outHDU->header) {
        outHDU->header = psMetadataAlloc();
    }


    // If input had WCS convert it for stamp
    if (input->fpa->toSky) {
        // steal the input fpa's transforms
        output->fpa->toTPA   = input->fpa->toTPA;
        output->fpa->fromTPA = input->fpa->fromTPA;
        output->fpa->toSky   = input->fpa->toSky;

        // drop the references
        input->fpa->toTPA   = 0;
        input->fpa->fromTPA = 0;
        input->fpa->toSky   = 0;

        if (!convertWCS(outHDU, output->fpa, outChip, inChip, &options->roi)) {
            return false;
        }
    } else {
        psWarning("No WCS present in input\n");
    }

    ppstampVersionMetadata(outHDU->header, options);

    return true;
}

// Extract pixels from an image. If part of the bounds of the region are
// partially outside of the input those pixels are set to zero.
psImage *extractStampOld(psImage *image, psRegion region)
{
    int width  = region.x1 - region.x0;
    int height = region.y1 - region.y0;

    if (width < 0) {
        return NULL;
    }
    if (height < 0) {
        return NULL;
    }

    int srcX  = region.x0;
    int srcY  = region.y0;
    int lastX = region.x1;
    int lastY = region.y1;
    if (lastX > (image->col0 + image->numCols)) {
        lastX = image->col0 + image->numCols;
    }
    if (lastY > (image->row0 + image->numRows)) {
        lastY = image->row0 + image->numRows;
    }

    int leftBlank = 0;
    int dstX  = 0;
    if (srcX < image->col0) {
        leftBlank = image->col0 - srcX;
        dstX += leftBlank;
        srcX = 0;
    }

    int copyWidth  = lastX - srcX;
    int rightBlank = width - copyWidth - leftBlank;
    if (copyWidth <= 0) {
        return NULL;
    }

    psImage *output = psImageAlloc(width, height, image->type.type);

    psElemType  type = image->type.type;
    psS32       elementSize =  PSELEMTYPE_SIZEOF(type);
    int         copySize = copyWidth * elementSize;

    for (int dstY=0 ; dstY < height; srcY++, dstY++) {
        psU8 *pdst = output->data.U8[dstY];

        if ((srcY < image->row0) || (srcY >= lastY)) {
            // zero the row
            memset(pdst, 0, width*elementSize);
        } else {
            psU8 *psrc = image->data.U8[srcY]  + (srcX * elementSize);

            if (leftBlank) {
                memset(pdst, 0, leftBlank*elementSize);
            }

            // copy copyWidth pixels from srcX,srcY to dstX, dstY
            pdst += dstX*elementSize;
            memcpy(pdst, psrc, copySize);

            if (rightBlank) {
                pdst += copyWidth * elementSize;
                memset(pdst, 0, rightBlank);
            }
        }
    }


    return output;
}
// Extract pixels from an image. If part of the bounds of the region are
// partially outside of the input those pixels are set to zero.
static psImage *extractStamp(psImage *image, psRegion region, double value)
{
    int width  = region.x1 - region.x0;
    int height = region.y1 - region.y0;

    if (width < 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "negative width\n");
        return NULL;
    }
    if (height < 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "negative height\n");
        return NULL;
    }

    int srcX  = region.x0;
    int srcY  = region.y0;
    int lastX = region.x1;
    int lastY = region.y1;
    if (lastX > (image->col0 + image->numCols)) {
        lastX = image->col0 + image->numCols;
    }
    if (lastY > (image->row0 + image->numRows)) {
        lastY = image->row0 + image->numRows;
    }

    int leftBlank = 0;
    int dstX  = 0;
    if (srcX < image->col0) {
        leftBlank = image->col0 - srcX;
        dstX += leftBlank;
        srcX = 0;
    }

    int copyWidth  = lastX - srcX;
//    int rightBlank = width - copyWidth - leftBlank;
    if (copyWidth <= 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "copyWidth is invalid: %d\n", copyWidth);
        return NULL;
    }

    psImage *output = psImageAlloc(width, height, image->type.type);
    if (!output) {
        psError(PS_ERR_UNKNOWN, false, "psImageAlloc failed\n");
        return NULL;
    }

    if (!psImageInit(output, value)) {
        psError(PS_ERR_UNKNOWN, false, "psInitImage failed\n");
        return NULL;
    }

    psElemType  type = image->type.type;
    psS32       elementSize =  PSELEMTYPE_SIZEOF(type);
    int         copySize = copyWidth * elementSize;

    int     dstY = 0;
    psS32   skip = image->row0 - srcY;
    if (skip > 0) {
        dstY = skip;
        height -= skip;
        srcY = image->row0;
    }
    for ( ; dstY < height; srcY++, dstY++) {
        if (srcY >= lastY) {
            break;
        }
        psU8 *pdst = output->data.U8[dstY];

        psU8 *psrc = image->data.U8[srcY]  + (srcX * elementSize);

        // copy copyWidth pixels from srcX,srcY to dstX, dstY
        pdst += dstX*elementSize;
        memcpy(pdst, psrc, copySize);
    }

    return output;
}

// Build the postage stamp output file

static bool makeStamp(pmConfig *config, ppstampOptions *options, pmFPAfile *input,
                pmChip *inChip, pmFPAview *view)
{
    int status = false;

    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTAMP.OUTPUT");
    if (!output) {
        psError(PS_ERR_UNKNOWN, false, "Can't find output data\n");
        return false;
    }
    char *fpaName = psMetadataLookupStr(NULL, input->fpa->concepts, "FPA.OBS"); // Name of FPA
    pmFPAview *outview = pmFPAviewAlloc(0);
    pmFPAAddSourceFromView(output->fpa, fpaName, outview, output->format);

    outview->chip = 0;
    outview->cell = 0;
    pmCell *target =  pmFPAviewThisCell(outview, output->fpa); // Target cell
    psFree(outview);
    outview = NULL;

    //   psMetadataPrint(stderr, inChip->concepts, 0);

    // These default to zero. would that be ok?
    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_META_REPLACE, "Binning in x", 1);
    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE, "Binning in y", 1);

    pmFPAfile *srcFile;
    pmFPAview *srcView = pmFPAviewAlloc(0);
    if (psArrayLength(inChip->cells) > 1) {
        // we extract our postage stamp from a mosaic
        pmFPAfile *mosaic = ppstampBuildMosaic(config, input, view);
        if (mosaic == NULL) {
            return false;
        }
        srcFile = mosaic;
        srcView->chip = 0;
    } else {
        srcFile = input;
        *srcView = *view;
    }

    // At this point we know we have one cell
    srcView->cell = 0;

    pmReadout *readout;
    while ((readout = pmFPAviewNextReadout (srcView, srcFile->fpa, 1)) != NULL) {
        if (!readout->data_exists) {
            psError(PS_ERR_UNKNOWN, false, "no data in input readout!\n");
            continue;
        }
        pmReadout *outReadout = pmReadoutAlloc(target);
        if (!outReadout) {
            psError(PS_ERR_UNKNOWN, false, "failed to allocate output readout\n");
            status = false;
            break;
        }

        psRegion extractRegion = options->roi;

        // Close your eyes while I hack around bug 986
        if (ppstampMegacamWorkaround) {
            // the coordinates of the mosaic are shifted 32 pixels from the chip
            // TODO does this always apply? For example I doubt that applies to
            // skycells.
            extractRegion.x0 -= 32;
            extractRegion.x1 -= 32;
        }

        outReadout->image = extractStamp(readout->image, extractRegion,  0);
        if (!outReadout->image) {
            psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp image\n");
            status = false;
            break;
        }
        if (readout->mask) {
            outReadout->mask = extractStamp(readout->mask, extractRegion,  0);
            if (!outReadout->mask) {
                psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp mask image\n");
                status = false;
                break;
            }
        }
        if (readout->weight) {
            outReadout->weight = extractStamp(readout->weight, extractRegion,  0);
            if (!outReadout->weight) {
                psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp weight image\n");
                status = false;
                break;
            }
        }

        outReadout->data_exists = true;
        outReadout->parent->data_exists = true;
        outReadout->parent->parent->data_exists = true;
        status = true;

        psFree(outReadout); // drop reference
    }

    psFree(srcView);

    if (status) {
        status = copyMetadata(output, input, inChip, options);
    }
    return status;
}



static bool regionContainsPoint(psRegion *r, psPlane *pt)
{
    if (pt->x < r->x0)
        return false;
    if (pt->x >= r->x1)
        return false;
    if (pt->y < r->y0)
        return false;
    if (pt->y >= r->y1)
        return false;

    return true;
}

// true if the inner region is equal to or completely contained in
// the outer region
static bool regionContainsRegion(psRegion *outer, psRegion *inner)
{
    if ((outer->x0 <= inner->x0) &&
        (outer->y0 <= inner->y0) &&
        (outer->x1 >= inner->x1) &&
        (outer->y1 >= inner->y1)) {
        return true;
    } else {
        return false;
    }
}


static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
{
    // convert from sky to TP to FP
    psProject(pt->TP, pt->sky, fpa->toSky);
    psPlaneTransformApply(pt->FP, fpa->fromTPA, pt->TP);
    // convert from FP to chip
    psPlaneTransformApply(pt->chip, chip->fromFPA, pt->FP);
}

static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
{
    // chip to FP
    psPlaneTransformApply(pt->FP, chip->toFPA, pt->chip);
    // FP to TP to sky
    psPlaneTransformApply(pt->TP, fpa->toTPA, pt->FP);
    psDeproject(pt->sky, pt->TP, fpa->toSky);

}

static void compareToBox(psRegion *region, psPlane *pt)
{
    if (pt->x < region->x0)
        region->x0 = pt->x;
    if (pt->x > region->x1)
        region->x1 = pt->x;
    if (pt->y < region->y0)
        region->y0 = pt->y;
    if (pt->y > region->y1)
        region->y1 = pt->y;
}

// For roi width and height given in sky coordinates find a bounding box in chip coordinates
// that encloses the requested range
static void findBoundingBox(ppstampOptions *options, pmFPA *fpa, pmChip *chip, pmAstromObj *center)
{
    pmAstromObj *pt = pmAstromObjAlloc();

    if (!options->roip.celestialCenter) {
        // Center was specified in chip coordinates, transform to the sky
        chipToSky(center, fpa, chip);
    }

    // calculate the four corners of the bounding box in sky coordinates, translate them to
    // chip coordinates and build the ROI by comparison.

    options->roi.x0 = INFINITY;
    options->roi.x1 = -INFINITY;
    options->roi.y0 = INFINITY;
    options->roi.y1 = -INFINITY;

    pt->sky->rErr = 0;
    pt->sky->dErr = 0;

    pt->sky->r = center->sky->r - options->roip.dRA/2;
    pt->sky->d = center->sky->d - options->roip.dDEC/2;
    skyToChip(pt, fpa, chip);
    compareToBox(&options->roi, pt->chip);

    pt->sky->r = center->sky->r + options->roip.dRA/2;
    pt->sky->d = center->sky->d - options->roip.dDEC/2;
    skyToChip(pt, fpa, chip);
    compareToBox(&options->roi, pt->chip);

    pt->sky->r = center->sky->r + options->roip.dRA/2;
    pt->sky->d = center->sky->d + options->roip.dDEC/2;
    skyToChip(pt, fpa, chip);
    compareToBox(&options->roi, pt->chip);

    pt->sky->r = center->sky->r - options->roip.dRA/2;
    pt->sky->d = center->sky->d + options->roip.dDEC/2;
    skyToChip(pt, fpa, chip);
    compareToBox(&options->roi, pt->chip);

    psFree(pt);

    // save the width and height in case we need them later
    options->roip.dX = options->roi.x1 - options->roi.x0;
    options->roip.dY = options->roi.y1 - options->roi.y0;
}

// findROI
// calculate the region of interest in chip coordinates and determine whether that region
// is completely contained on a single chip

static ppstampOverlap findROI(ppstampOptions *options, pmFPAview *view,
                    pmFPAfile *input, pmFPAfile *astrom, bool bilevelAstrometry,
                    pmChip *chip, pmAstromObj *center)
{
    psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
    psRegion    *chipBounds = ppstampChipRegion(chip);
    bool        onChip = false;
    ppstampOverlap   returnval = PPSTAMP_OFF;

//    fprintf(stderr, "ppstampChipBounds: %s\n", psRegionToString(*chipBounds));

    // set up the astrometry
    pmHDU *hdu = pmFPAviewThisHDU(view, astrom->fpa);
    PS_ASSERT_PTR_NON_NULL(hdu, 1)
    PS_ASSERT_PTR_NON_NULL(hdu->header, 1)

    // we can live without astrometry if the ROI is completely specified in pixel coordinates
    bool mustHaveAstrometry = false;
    if (options->roip.celestialCenter || options->roip.celestialRange) {
        mustHaveAstrometry = true;
    }
    if (bilevelAstrometry) {
        if (!pmAstromReadBilevelChip (chip, hdu->header) && mustHaveAstrometry) {
            psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel chip astrometry for input FPA.");
            return false;
        }
    } else {
        // we use a default FPA pixel scale of 1.0
        if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0) && mustHaveAstrometry) {
            psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA.");
            return false;
        }
    }

    if (options->roip.celestialCenter) {

        center->sky->r = options->roip.centerRA;
        center->sky->d = options->roip.centerDEC;
        center->sky->rErr = 0;
        center->sky->dErr = 0;

        skyToChip(center, input->fpa, chip);

        if (regionContainsPoint(chipBounds, center->chip)) {
            psLogMsg("ppstampMakeStamp", 2, "Found center (%f %f) on chip: %s\n",
                center->chip->x, center->chip->y, chipName);
            onChip = true;
        }
    } else {
        // center specified in pixels.
        // If the user specified a name of a chip name wait until we get to that one.
        // If no chip name was specified, select this one (the first one that had data)
        if ((options->chipName == NULL) || !strcasecmp(chipName, options->chipName)) {
            psLogMsg("ppstampMakeStamp", 2, "Center on chip: %s\n", chipName);
            center->chip->x = options->roip.centerX;
            center->chip->y = options->roip.centerY;
            center->chip->xErr = 0;
            center->chip->yErr = 0;
            onChip = true;
        }
    }

    if (onChip) {
        if (options->roip.celestialRange) {
            findBoundingBox(options, input->fpa, chip, center);
        } else {
            int width  = options->roip.dX;
            int height = options->roip.dY;

            // calculate the ROI in chip coordinates
            options->roi.x0 = center->chip->x - width / 2;
            options->roi.x1 = options->roi.x0 + width;
            options->roi.y0 = center->chip->y - height / 2;
            options->roi.y1 = options->roi.y0 + height;
        }


        if (regionContainsRegion(chipBounds, &options->roi)) {
            // returnval = findCell(view, options, input, center, ppCell);
            psLogMsg("ppstampMakeStamp", 2, "ROI contained on: %s\n", chipName);
            returnval = PPSTAMP_ON;
        } else {
            psLogMsg("ppstampMakeStamp", 2, "Partial Overlap chip %s\n", chipName);
            psLogMsg("ppstampMakeStamp", 2, "ROI:         %s\n", psRegionToString(options->roi));
            psLogMsg("ppstampMakeStamp", 2, "Chip Extent: %s\n", psRegionToString(*chipBounds));

            returnval = PPSTAMP_PARTIALLY_ON;
        }
    }
    psFree(chipBounds);

    return returnval;
}

bool ppstampMakeStamp (pmConfig *config, ppstampOptions *options)
{
    bool        status = false;
    bool        returnval = false;;
    bool        foundOverlap = false;

    pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PPSTAMP.INPUT");
    if (!status) {
        psError(PS_ERR_UNKNOWN, true, "Can't find input file!\n");
        return false;
    }

    pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
    if (!astrom) {
        astrom = input;
    } else if (astrom->camera != input->camera) {
        psError(PS_ERR_UNKNOWN, true, "Input camera and astrometry camera do not match");
        return false;
    }

    pmFPAview *view = pmFPAviewAlloc(0);// View for level of interest

    // files associated with the science image
    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
        psError(PS_ERR_UNKNOWN, false, "Failed to load input.");
        psFree (view);
        return false;
    }
    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);
            return false;
        }
    }

    // Loop over the chips and find the one that contains the center
    pmAstromObj *center = pmAstromObjAlloc();
    pmChip *chip;
    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
        bool allDone = false;

        if (!chip->process || !chip->file_exists) {
            continue;
        }

        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
            psError(PS_ERR_UNKNOWN, false, "failed to load chip");
            status = false;
            break;
        }

        ppstampOverlap overlap = findROI(options, view, input, astrom, bilevelAstrometry, chip, center);

        switch (overlap) {
        case PPSTAMP_OFF:
            // keep looking
            break;
        case PPSTAMP_ON:
        case PPSTAMP_PARTIALLY_ON:
            returnval = makeStamp(config, options, input, chip, view);
            allDone = true;
            foundOverlap = true;
            break;
        case PSTAMP_ERROR:
            returnval = false;
            allDone = true;
            break;
        default:
            psError(PS_ERR_PROGRAMMING, false, "findROI returned unexpected value %d\n", overlap);
            break;
        }

        pmFPAfileIOChecks(config, view, PM_FPA_AFTER);

        if (allDone) {
            view->chip = -1;
            break;
        }
    }
    pmFPAfileIOChecks(config, view, PM_FPA_AFTER);

    psFree(center);
    psFree(view);

    if (!foundOverlap) {
        fprintf(stderr, "ROI not found in input\n");
    }

    return returnval;
}


