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

# include "psphot.h"

// generate the median in NxN boxes, clipping heavily
// linear interpolation to generate full-scale model
bool psphotImageMedian (pmConfig *config, pmFPAview *view)
{
    bool status;
    pmFPA *inFPA;
    pmFPAfile *file;
    psRegion region;
    int MAX_SAMPLE_PIXELS;

    pmReadout *model = NULL;
    pmReadout *readout = NULL;
    pmReadout *background = NULL;
    pmReadout *backSub = NULL;

    psTimerStart ("psphot");

    // select the appropriate recipe information
    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");

    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
    if (!status) {
        MAX_SAMPLE_PIXELS = 1000;
    }
    unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED");
    if (!status) {
        seed = 0;
    }
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, seed);

    // subtract this amount extra from the sky
    float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS");
    if (!status) {
        SKY_BIAS = 0;
    }

    // find the currently selected readout
    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
    inFPA = file->fpa;
    readout = pmFPAviewThisReadout (view, inFPA);

    psImage *image = readout->image;
    psImage *mask  = readout->mask;

    // dimensions of input & output image
    int Nx = image->numCols;
    int Ny = image->numRows;

    // scaling factor
    int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
    int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");

    // overhang : we will balance this evenly
    int dx = (Nx % DX) / 2;
    int dy = (Ny % DY) / 2;

    // dimensions of binned image
    int nx = (Nx % DX) ? (int)(Nx / DX) + 1 : Nx / DX;
    int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY;

    // select model pixels (from output background model file, or create internal file)
    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL");
    if (file == NULL) {
	// we are not using PSPHOT.BACKMDL as an I/O file: define an internal version
        model = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
    } else {
	// we are using PSPHOT.BACKMDL as an I/O file: select readout or create
	model = pmFPAviewThisReadout (view, file->fpa);
	if (model == NULL) {
	    // readout does not yet exist: create from input
	    pmFPAfileCopyStructureView (file->fpa, inFPA, DX, DY, view);
	    model = pmFPAviewThisReadout (view, file->fpa);
	    if (nx != model->image->numCols) psAbort ("psphot", "inconsistent size");
	    if (ny != model->image->numRows) psAbort ("psphot", "inconsistent size");
	}
    }
    psF32 **modelData = model->image->data.F32;
    assert(model->analysis != NULL);
    psMetadataAdd(model->analysis, PS_LIST_TAIL, "XBIN", PS_DATA_S32 | PS_META_REPLACE,
		  "Background x-binsize", DX);
    psMetadataAdd(model->analysis, PS_LIST_TAIL, "YBIN", PS_DATA_S32 | PS_META_REPLACE,
		  "Background x-binsize", DY);
    psMetadataAdd(model->analysis, PS_LIST_TAIL, "x0", PS_DATA_S32 | PS_META_REPLACE,
		  "Background x-overhang", dx);
    psMetadataAdd(model->analysis, PS_LIST_TAIL, "y0", PS_DATA_S32 | PS_META_REPLACE,
		  "Background y-overhang", dy);

    // measure clipped median for subimages
    for (int iy = 0; iy < ny; iy++) {
        for (int ix = 0; ix < nx; ix++) {
            // sx, sy are in parent coords
            int sx = ix*DX - dx + image->col0;
            int sy = iy*DY - dy + image->row0;
            region = psRegionSet (sx, sx + 2*DX, sy, sy + 2*DY);
            region = psRegionForImage (image, region);
            psImage *subset  = psImageSubset (image, region);
	    if (!subset->numCols || !subset->numRows) {
		psFree (subset);
		continue;
	    }

            psImage *submask = psImageSubset (mask, region);

            // XXX the value of the upper and lower cuts probably should be studied...
            psStats *stats = psImageBackground(subset, submask, 0xff, 0.25, 0.75, MAX_SAMPLE_PIXELS, rng);
            modelData[iy][ix] = stats->robustMedian + SKY_BIAS;

            psFree (stats);
            psFree (subset);
            psFree (submask);

            // XXX psImageStats is still very inefficient and poorly coded...
            // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
            // stats = psImageStats (stats, subset, maskset, 0xff);
            // backMdl->data.F32[iy][ix] = stats->clippedMean;
        }
    }
    psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot"));

    // measure background stats and save for later output
    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_MIN | PS_STAT_MAX);
    stats = psImageStats (stats, model->image, NULL, 0);
    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MEAN", PS_DATA_F32 | PS_META_REPLACE, "sky model mean",          stats->sampleMean);
    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_STDV", PS_DATA_F32 | PS_META_REPLACE, "sky model stdev",         stats->sampleStdev);
    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MAX",  PS_DATA_F32 | PS_META_REPLACE, "sky model maximum value", stats->max);
    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MIN",  PS_DATA_F32 | PS_META_REPLACE, "sky model minimum value", stats->min);
    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_NX",   PS_DATA_S32 | PS_META_REPLACE, "sky model size (x)",      nx);
    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_NY",   PS_DATA_S32 | PS_META_REPLACE, "sky model size (y)",      ny);
    psFree (stats);

    // select background pixels, from output background file, or create
    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
    if (file) {
	// we are using PSPHOT.BACKMDL as an I/O file: select readout or create
	background = pmFPAviewThisReadout (view, file->fpa);
	if (background == NULL) {
	    // readout does not yet exist: create from input
	    pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view);
	    background = pmFPAviewThisReadout (view, file->fpa);
	    if (Nx != background->image->numCols) psAbort ("psphot", "inconsistent size");
	    if (Ny != background->image->numRows) psAbort ("psphot", "inconsistent size");
	}
    } else {
        background = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
    }
    psF32 **backData = background->image->data.F32;

    // linear interpolation to full-scale
    psImageUnbin (background->image, model->image, DX, DY, dx, dy);
    psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot"));

    // back-sub image pixels, from output background file (don't create if not requested)
    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKSUB");
    if (file) {
	// we are using PSPHOT.BACKMDL as an I/O file: select readout or create
	backSub = pmFPAviewThisReadout (view, file->fpa);
	if (backSub == NULL) {
	    // readout does not yet exist: create from input
	    pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view);
	    backSub = pmFPAviewThisReadout (view, file->fpa);
	}
    }

# ifdef TESTSAVE
    psphotSaveImage (NULL, image, "image.fits");
    psphotSaveImage (NULL, background->image, "back.fits");
    psphotSaveImage (NULL, mask, "mask.fits");
    psphotSaveImage (NULL, model->image, "backmdl.fits");
# endif

    // subtract the background model (save in backSub, if requested)
    for (int j = 0; j < image->numRows; j++) {
        for (int i = 0; i < image->numCols; i++) {
	    image->data.F32[j][i] -= backData[j][i];
	    if (backSub) {
		backSub->image->data.F32[j][i] = image->data.F32[j][i];
	    }
        }
    }

# ifdef TESTSAVE
    psphotSaveImage (NULL, image, "backsub.fits");
# endif

    psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
    psFree(rng);

    // the pmReadout selected in this function are all view on entries in config->files
    return true;
}
