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

// we have 4 possibilities: (INTERNAL or I/O file) and (exists or not) 
// select model pixels (from output background model file, or create internal file)
static pmReadout *get_model_readout(const char *name, // name of internal/external file
				    const pmConfig *config, // configuration information
				    pmFPAview *view,
				    pmFPA *inFPA,
				    const psImageBinning *binning) {
    pmReadout *model = NULL;

    bool status = true;
    pmFPAfile *file = psMetadataLookupPtr(&status, config->files, name);
    if (file == NULL) {
	// we are not using PSPHOT.BACKMDL as an I/O file: define an internal version
        model = pmFPAfileDefineInternal (config->files, name, binning->nXruff, binning->nYruff, PS_TYPE_F32);
    } else {
	if (file->mode == PM_FPA_MODE_INTERNAL) {
	    // we are not using PSPHOT.BACKMDL as an I/O file: already defined above
	    model = file->readout;
	} 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
		// XXX we have an inconsistency in this calculation here and in pmFPACopy
		// XXX use the psImageBinning functions to set the output image size
		pmFPAfileCopyStructureView (file->fpa, inFPA, binning->nXbin, binning->nYbin, view);
		model = pmFPAviewThisReadout (view, file->fpa);
		PS_ASSERT (binning->nXruff == model->image->numCols, false);
		PS_ASSERT (binning->nYruff == model->image->numRows, false);
	    }
	}
    }

    return model;
}

// generate the median in NxN boxes, clipping heavily
// linear interpolation to generate full-scale model
bool psphotImageMedian (pmConfig *config, pmFPAview *view)
{
    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))) {
	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 {
	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
    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 = get_model_readout("PSPHOT.BACKMDL", config, view, inFPA, binning);
    pmReadout *modelStdev = get_model_readout("PSPHOT.BACKMDL.STDEV", config, view, 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...

	    if (psImageBackground(stats, subset, submask, 0xff, 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);
	    } else {
		psStatsOptions currentOptions = stats->options;
		stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
		if (!psImageBackground(stats, subset, submask, 0xff, 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);
		}
		// 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);
        }
    }

    // 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;
	}
    }
    assert (Count > 0);
    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, "SKY_MODEL_MEAN",
		      PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_STDEV",
		      PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_MAX",
		      PS_META_REPLACE, "sky model maximum value", statsBck->max);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_MIN",
		      PS_META_REPLACE, "sky model minimum value", statsBck->min);
    psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_MODEL_NX",
		      PS_META_REPLACE, "sky model size (x)",      model->image->numCols);
    psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_MODEL_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, "failed to build background image");
	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;
}
