# include "psphotInternal.h"
static int npass = 0;

// generate the median in NxN boxes, clipping heavily
// linear interpolation to generate full-scale model
bool psphotImageMedian (pmConfig *config, const pmFPAview *view, psMaskType maskVal)
{
    bool status = true;
    static char *defaultStatsName = "FITTED_MEAN";

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

    psTimerStart ("psphot");

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

    // user supplied seed, if available
    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;
    }

    // supply the sky background statistics options
    char *statsName = psMetadataLookupStr (&status, recipe, "SKY_STAT");
    if (statsName == NULL) {
        statsName = defaultStatsName;
    }
    psStatsOptions statsOptionLocation = psStatsOptionFromString(statsName);
    if (!(statsOptionLocation & (PS_STAT_SAMPLE_MEAN |
                                 PS_STAT_SAMPLE_MEDIAN |
                                 PS_STAT_ROBUST_MEDIAN |
                                 PS_STAT_ROBUST_QUARTILE |
                                 PS_STAT_CLIPPED_MEAN |
                                 PS_STAT_FITTED_MEAN |
                                 PS_STAT_FITTED_MEAN_V2 |
                                 PS_STAT_FITTED_MEAN_V3 |
                                 PS_STAT_FITTED_MEAN_V4))) {
        statsOptionLocation = PS_STAT_FITTED_MEAN;
    }

    psStatsOptions statsOptionWidth = PS_STAT_NONE;
    if (statsOptionLocation & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN)) {
        statsOptionWidth = PS_STAT_SAMPLE_STDEV;
    } else if (statsOptionLocation & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE)) {
#if 1
        statsOptionWidth = PS_STAT_ROBUST_STDEV; // not set; => NaN
#else
        statsOptionWidth = PS_STAT_FITTED_STDEV;
#endif
    } else if (statsOptionLocation & PS_STAT_FITTED_MEAN) {
        statsOptionWidth = PS_STAT_FITTED_STDEV;
    } else if (statsOptionLocation & PS_STAT_CLIPPED_MEAN) {
        statsOptionWidth = PS_STAT_CLIPPED_STDEV;
    } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V2) {
        statsOptionWidth = PS_STAT_FITTED_STDEV_V2;
    } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V3) {
        statsOptionWidth = PS_STAT_FITTED_STDEV_V3;
    } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V4) {
        statsOptionWidth = PS_STAT_FITTED_STDEV_V4;
    } else {
        psAbort("Unable to estimate variance of selected statsOptionLocations 0x%x", statsOptionLocation);
    }

    const psStatsOptions statsOption = statsOptionLocation | statsOptionWidth;
    psStats *stats = psStatsAlloc (statsOption);

    // set range for old-version of sky statistic
    if (statsOptionLocation & PS_STAT_ROBUST_QUARTILE) {
        stats->min = 0.25;
        stats->max = 0.75;
    }

    // set user-option for number of pixels per region
    stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
    if (!status) {
        stats->nSubsample = 1000;
    }

    // optionally set the binsize
    stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
    if (status) {
        stats->options |= PS_STAT_USE_BINSIZE;
    }

    // optionally set the sigma clipping
    stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
    if (!status) {
        if ((stats->options & PS_STAT_FITTED_MEAN) || (stats->options & PS_STAT_FITTED_MEAN_V2)) {
            stats->clipSigma = 1.0;
        } else {
            stats->clipSigma = 3.0;
        }
    }

    // stats is not initialized by psStats???  use this to save the input options
    // XXX re-work this to use the new psStatsInit function
    psStats *statsDefaults = psStatsAlloc (statsOption);
    *statsDefaults = *stats;

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

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

    // I have the fine image size, I know the binning factor, determine the ruff image size
    psImageBinning *binning = psImageBinningAlloc();
    binning->nXfine = image->numCols;
    binning->nYfine = image->numRows;
    binning->nXbin  = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
    binning->nYbin  = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");

    psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER);
    psImageBinningSetSkip(binning, image);
    status = psMetadataAddPtr(recipe, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
    PS_ASSERT (status, false);

    // we save the binning structure for use in psphotMagnitudes
    pmReadout *model      = pmFPAGenerateReadout (config, view, "PSPHOT.BACKMDL",       inFPA, binning);
    pmReadout *modelStdev = pmFPAGenerateReadout (config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning);

    psF32 **modelData = model->image->data.F32;
    psF32 **modelStdevData = modelStdev->image->data.F32;

    // measure clipped median for subimages
    psRegion ruffRegion = {0,0,0,0};
    psRegion fineRegion = {0,0,0,0};
    for (int iy = 0; iy < model->image->numRows; iy++) {
        for (int ix = 0; ix < model->image->numCols; ix++) {

            // convert the ruff grid cell to the equivalent fine grid cell
            // XXX we need to watch out for row0,col0
            ruffRegion = psRegionSet (ix, ix + 2, iy, iy + 2);
            fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
            fineRegion = psRegionForImage (image, fineRegion);

            psImage *subset  = psImageSubset (image, fineRegion);
            if (!subset->numCols || !subset->numRows) {
                psFree (subset);
                continue;
            }
            psImage *submask = psImageSubset (mask, fineRegion);

            // reset the default values
            *stats = *statsDefaults;

            // Use the selected background statistic for the first pass
            // If it fails, fall back on the "ROBUST_MEDIAN" version
            // If both fail, set the pixel to NAN and (below) interpolate
            // XXX psImageBackground will probably be renamed psImageStats
            // XXX don't bother trying if there are no valid pixels...

	    psVector *sample = NULL;

	    // turn on stats tracing in desired cells
	    # if (0)
	    psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS");
	    assert (plots);
	    
	    int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X");
	    assert (status);
	    int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y");
	    assert (status);

	    bool gotX = (xPlot < 0) || (xPlot == ix);
	    bool gotY = (yPlot < 0) || (yPlot == iy);

	    if (gotX && gotY) {
		psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 6);
		psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
	    } else {
		psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 0);
		psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
	    }	    
	    # endif

            if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
                if (stats->options & PS_STAT_ROBUST_QUARTILE) {
                    modelData[iy][ix] = stats->robustMedian;
                } else {
                    modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
                }
                modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);

		// supply sample to plotting routing
		psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
		psFree (sample);
            } else {
                psStatsOptions currentOptions = stats->options;
                stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
                if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
                    psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
                               "(%dx%d, (row0,col0) = (%d,%d)",
                               subset->numRows, subset->numCols, subset->row0, subset->col0);
                    modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
                } else {
                    modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
                    modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);

		    // supply sample to plotting routing
		    psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
		    psFree (sample);
                }
                // drop errors caused by psImageBackground failures
                // XXX we probably should trap and exit on serious failures
                psErrorClear();
                stats->options = currentOptions;
            }
            modelData[iy][ix] += SKY_BIAS;
            psFree (subset);
            psFree (submask);
        }
    }

    if (psTraceGetLevel("psphot") > 5) {
        char name[256];
        sprintf (name, "backraw.%02d.fits", npass);
        psphotSaveImage (NULL, model->image, name);
    }

    // patch over bad regions (use average of 8 possible neighbor pixels)
    // XXX consider testing all pixels against the 8 neighbors and replacing outliers...
    double Count = 0;                   // number of good pixels
    double Value = 0;                   // sum of good pixel's value
    double ValueStdev = 0;              // sum of good pixel's standard deviations
    for (int iy = 0; iy < model->image->numRows; iy++) {
        for (int ix = 0; ix < model->image->numCols; ix++) {
            if (!isnan(modelData[iy][ix])) {
                Value += modelData[iy][ix];
                ValueStdev += modelStdevData[iy][ix];
                Count++;
                continue;
            }

            double value = 0;
            double count = 0;
            for (int jy = iy - 1; jy <= iy + 1; jy++) {
                if (jy <   0) continue;
                if (jy >= model->image->numRows) continue;
                for (int jx = ix - 1; jx <= ix + 1; jx++) {
                    if (!jx && !jy) continue;
                    if (jx   <   0) continue;
                    if (jx   >= model->image->numCols) continue;
                    value += modelData[jy][jx];
                    count += 1.0;
                }
            }
            if (count > 0) modelData[iy][ix] = value / count;
        }
    }
    if (Count == 0) {
        psError (PSPHOT_ERR_DATA, true, "failed to build background image");
	psFree(stats);
	psFree(statsDefaults);
	psFree(binning);
	psFree(rng);
        return false;
    }

    Value /= Count;
    ValueStdev /= Count;

    // patch over remaining bad regions (use global average)
    for (int iy = 0; iy < model->image->numRows; iy++) {
        for (int ix = 0; ix < model->image->numCols; ix++) {
            if (!isnan(modelData[iy][ix])) continue;
            modelData[iy][ix] = Value;
            modelStdevData[iy][ix] = ValueStdev;
        }
    }

    psLogMsg ("psphot", PS_LOG_MINUTIA, "build median image: %f sec\n", psTimerMark ("psphot"));

    psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value);
    psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
    psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev);

    // measure image and background stats and save for later output
    psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN |
                                      PS_STAT_SAMPLE_STDEV |
                                      PS_STAT_MIN |
                                      PS_STAT_MAX);
    psImageStats (statsBck, model->image, NULL, 0);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MN",
                      PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_SIG",
                      PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MAX",
                      PS_META_REPLACE, "sky model maximum value", statsBck->max);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MIN",
                      PS_META_REPLACE, "sky model minimum value", statsBck->min);
    psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NX",
                      PS_META_REPLACE, "sky model size (x)",      model->image->numCols);
    psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NY",
                      PS_META_REPLACE, "sky model size (y)",      model->image->numRows);
    psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f",
              statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
    psFree (statsBck);

    // select background pixels, from output background file, or create
    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
    if (file) {
        // we are using PSPHOT.BACKGND as an I/O file: select readout or create
        if (file->mode == PM_FPA_MODE_INTERNAL) {
            background = file->readout;
        } else {
            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 ((image->numCols != background->image->numCols) || (image->numRows != background->image->numRows)) {
                psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for background dimensions");
                return false;
            }
        }
    } else {
        background = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKGND", image->numCols, image->numRows, PS_TYPE_F32);
    }
    psF32 **backData = background->image->data.F32;

    // linear interpolation to full-scale
    if (!psImageUnbin (background->image, model->image, binning)) {
        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
        return false;
    }

    psLogMsg ("psphot", PS_LOG_MINUTIA, "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.BACKSUB 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);
        }
    }

    if (psTraceGetLevel("psphot") > 5) {
        char name[256];
        sprintf (name, "image.%02d.fits", npass);
        psphotSaveImage (NULL, image, name);
        sprintf (name, "back.%02d.fits", npass);
        psphotSaveImage (NULL, background->image, name);
        sprintf (name, "mask.%02d.fits", npass);
        psphotSaveImage (NULL, mask, name);
        sprintf (name, "backmdl.%02d.fits", npass);
        psphotSaveImage (NULL, model->image, name);
    }

    // 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];
            }
        }
    }

    if (psTraceGetLevel("psphot") > 5) {
        char name[256];
        sprintf (name, "backsub.%02d.fits", npass);
        psphotSaveImage (NULL, image, name);
    }
    npass ++;

    psLogMsg ("psphot", PS_LOG_INFO, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
    psFree(stats);
    psFree(statsDefaults);
    psFree(binning);
    psFree(rng);

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