# include "psphotInternal.h"
static int npass = 0;
static char *defaultStatsName = "FITTED_MEAN";

// Determine the binning characteristics for the background model
// Sets the ruff image size and skip based on recipe values for the binning
static psImageBinning *backgroundBinning(const psImage *image, // Image for which to generate a bg model
                                         const pmConfig *config // Configuration
                                         )
{
    bool status = true;

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

    // 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);

    return binning;
}

// Generate the background model
static bool backgroundModel(psImage *model,  // Model image
                            psImage *modelStdev, // Model stdev image
                            psMetadata *analysis, // Analysis metadata for outputs
                            pmReadout *readout, // Readout for which to generate a background model
                            psImageBinning *binning, // Binning parameters
                            const pmConfig *config // Configuration
                            )
{
    psTimerStart ("psphot.background");

    bool status = true;

    psImage *image = readout->image, *mask = readout->mask; // Image and mask for readout

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

    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    assert (maskVal);

    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);

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

    // we save the binning structure for use in psphotMagnitudes
    psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING",
                     PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);

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

    // XXXX we can thread this here by running blocks in parallel

    // measure clipped median for subimages
    psRegion ruffRegion = {0,0,0,0};
    psRegion fineRegion = {0,0,0,0};
    for (int iy = 0; iy < model->numRows; iy++) {
        for (int ix = 0; ix < model->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 + 1, iy, iy + 1);
            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) {
                (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 6);
                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
            } else {
                (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 0);
                (void) 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);
            } 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);
                }
                // drop errors caused by psImageBackground failures
                // XXX we probably should trap and exit on serious failures
                psErrorClear();
                stats->options = currentOptions;
            }
            psFree(sample);
            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, 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->numRows; iy++) {
        for (int ix = 0; ix < model->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->numRows) continue;
                for (int jx = ix - 1; jx <= ix + 1; jx++) {
                    if (!jx && !jy) continue;
                    if (jx   <   0) continue;
                    if (jx   >= model->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->numRows; iy++) {
        for (int ix = 0; ix < model->numCols; ix++) {
            if (!isnan(modelData[iy][ix])) continue;
            modelData[iy][ix] = Value;
            modelStdevData[iy][ix] = ValueStdev;
        }
    }

    psLogMsg ("psphot", PS_LOG_INFO, "built background image: %f sec\n", psTimerMark ("psphot.background"));

    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, 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->numCols);
    psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NY",
                      PS_META_REPLACE, "sky model size (y)",      model->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(stats);
    psFree(statsBck);
    psFree(statsDefaults);
    psFree(binning);
    psFree(rng);

    return model;
}



psImage *psphotBackgroundModel(pmReadout *ro, const pmConfig *config)
{
    PM_ASSERT_READOUT_NON_NULL(ro, NULL);
    PM_ASSERT_READOUT_IMAGE(ro, NULL);
    PS_ASSERT_PTR_NON_NULL(config, NULL);

    psImageBinning *binning = backgroundBinning(ro->image, config); // Image binning parameters
    psImage *model = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Background model
    psImage *modelStdev = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Standard deviation

    if (!backgroundModel(model, modelStdev, ro->analysis, ro, binning, config)) {
        psFree(model);
        psFree(modelStdev);
        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
        return NULL;
    }

    psFree(modelStdev);

    return model;
}


// generate the median in NxN boxes, clipping heavily
// linear interpolation to generate full-scale model
bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filename)
{
    bool status = true;

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

    psImageBinning *binning = backgroundBinning(readout->image, config); // Image binning parameters
    pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning);
    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning);

    if (!backgroundModel(model->image, modelStdev->image, model->analysis, readout, binning, config)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
        return false;
    }

    npass ++;
    return true;
}
