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

#include <stdio.h>
#include <assert.h>
#include <string.h>
#include <pslib.h>

#include "pmHDU.h"
#include "pmFPA.h"
#include "pmHDUUtils.h"
#include "pmFPALevel.h"
#include "pmFPAview.h"
#include "pmFPACalibration.h"

#include "pmOverscan.h"
#include "pmBias.h"

// pmBiasSubtractFrame(): this routine will take as input a readout for the input image and a readout for the bias
// image.  The bias image is subtracted in place from the input image.
bool pmBiasSubtractFrame(pmReadout *in, // Input readout
                          const pmReadout *sub, // Readout to be subtracted from input
                          float scale   // Scale to apply before subtracting
                         )
{
    PS_ASSERT_PTR_NON_NULL(in, false);
    PS_ASSERT_PTR_NON_NULL(in->image, false);
    PS_ASSERT_IMAGE_TYPE(in->image, PS_TYPE_F32, false);
    PS_ASSERT_IMAGE_NON_EMPTY(in->image, false);
    PS_ASSERT_PTR_NON_NULL(sub, false);
    PS_ASSERT_PTR_NON_NULL(sub->image, false);
    PS_ASSERT_IMAGE_TYPE(sub->image, PS_TYPE_F32, false);
    PS_ASSERT_IMAGE_NON_EMPTY(sub->image, false);
    PS_ASSERT_IMAGES_SIZE_EQUAL(in->image, sub->image, false);

    psImage *inImage  = in->image;      // The input image
    psImage *inMask   = in->mask;       // The input mask
    psImage *subImage = sub->image;     // The image to be subtracted
    psImage *subMask  = sub->mask;      // The mask for the subtraction image

    // Check parities
    int xIpar = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.XPARITY");
    int xSpar = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.XPARITY");
    if (xIpar != xSpar) {
        psError(PS_ERR_UNKNOWN, true, "images for subtraction do not have the same "
                "CELL.XPARITY (%d vs %d).\n    pmSubtractBias must be upgraded to handle this situation\n",
                xIpar, xSpar);
        return false;
    }

    int yIpar = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.YPARITY");
    int ySpar = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.YPARITY");
    if (yIpar != ySpar) {
        psError(PS_ERR_UNKNOWN, true, "images for subtraction do not have the same "
                "CELL.YPARITY (%d vs %d).\n    pmSubtractBias must be upgraded to handle this situation\n",
                xIpar, xSpar);
        return false;
    }

    // Offsets of the cells
    int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0");
    int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0");
    int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0");
    int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0");

    if ((inImage->numCols + x0in - x0sub) > subImage->numCols) {
        psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction (%d vs %d).\n",
                inImage->numCols + x0in - x0sub, subImage->numCols);
        return false;
    }
    if ((inImage->numRows + y0in - y0sub) > subImage->numRows) {
        psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction (%d vs %d).\n",
                inImage->numRows + y0in - y0sub, subImage->numRows);
        return false;
    }

    if (scale == 1.0) {
        for (int i = 0; i < inImage->numRows; i++) {
            for (int j = 0; j < inImage->numCols; j++) {
                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub];
                if (inMask && subMask) {
                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
                }
            }
        }
    } else {
        for (int i = 0; i < inImage->numRows; i++) {
            for (int j = 0; j < inImage->numCols; j++) {
                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale;
                if (inMask && subMask) {
                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
                }
            }
        }
    }

    return true;
}

bool pmBiasSubtract(pmReadout *in, pmOverscanOptions *overscanOpts,
                    const pmReadout *bias, const pmReadout *dark, const pmFPAview *view)
{
    psTrace("psModules.detrend", 4,
            "---- pmBiasSubtract() begin ----\n");

    PS_ASSERT_PTR_NON_NULL(in, false);
    PS_ASSERT_IMAGE_NON_NULL(in->image, false);
    PS_ASSERT_IMAGE_TYPE(in->image, PS_TYPE_F32, false);
    if (bias) {
        PS_ASSERT_IMAGE_NON_NULL(bias->image, false);
        PS_ASSERT_IMAGE_TYPE(bias->image, PS_TYPE_F32, false);
    }
    if (dark) {
        PS_ASSERT_PTR_NON_NULL(view, false);
        PS_ASSERT_IMAGE_NON_NULL(dark->image, false);
        PS_ASSERT_IMAGE_TYPE(dark->image, PS_TYPE_F32, false);
    }

    pmHDU *hdu = pmHDUFromReadout(in);  // HDU of interest

    if (!pmOverscanSubtract (in, overscanOpts)) {
	return false;
    }

    // Bias frame subtraction
    if (bias) {
        psVector *md5 = psImageMD5(bias->image); // md5 hash
        psString md5string = psMD5toString(md5); // String
        psFree(md5);
        psStringPrepend(&md5string, "BIAS image MD5: ");
        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
                         md5string, "");
        psFree(md5string);

        if (!pmBiasSubtractFrame(in, bias, 1.0)) {
            return false;
        }
    }

    if (dark) {
        // Get the scaling
        float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME");
        float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME");
        if (isnan(inTime) || isnan(darkTime)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to determine dark scaling.");
            return false;
        }

	float darkNorm = 1.0;
        float inNorm = pmFPADarkNorm(in->parent->parent->parent, view, inTime);

	// if we have a normalized dark exposure, we simply multiply the master by inNorm.  if
	// we do not have a normalized exposure, we have to scale the master as well.  XXX do
	// we need to explicitly identify the master as normalized?

	if (darkTime != 1.0) {
	    darkNorm = pmFPADarkNorm(dark->parent->parent->parent, view, darkTime);
	}

	if (isnan(inNorm) || isnan(darkNorm)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to determine dark normalisations.");
            return false;
        }

        float scale = inNorm / darkNorm;// Scaling to apply to dark exposure

        psVector *md5 = psImageMD5(dark->image); // md5 hash
        psString md5string = psMD5toString(md5); // String
        psFree(md5);
        psStringPrepend(&md5string, "DARK image (scale %.3f) MD5: ", scale);
        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
                         md5string, "");
        psFree(md5string);

        if (!pmBiasSubtractFrame(in, dark, scale)) {
            return false;
        }
    }

    psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now, used for reporting
    psString timeString = psTimeToISO(time); // String with time
    psFree(time);
    psStringPrepend(&timeString, "Overscan/bias/dark processing completed at ");
    psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
                     timeString, "");
    psFree(timeString);


    return true;
}


