Index: /trunk/psphot/src/psphotImageMedian.c
===================================================================
--- /trunk/psphot/src/psphotImageMedian.c	(revision 13012)
+++ /trunk/psphot/src/psphotImageMedian.c	(revision 13013)
@@ -2,107 +2,18 @@
 static int npass = 0;
 
-// 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;
-    static char *defaultStatsName = "FITTED_MEAN";
-
+// 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;
-    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 statsOption = psStatsOptionFromString (statsName);
-    if (!(statsOption & (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))) {
-	statsOption = PS_STAT_FITTED_MEAN;
-    }
-    psStats *stats = psStatsAlloc (statsOption);
-
-    // set range for old-version of sky statistic
-    if (statsOption & 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
-    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
-    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
-
-    // 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)
-    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL");
+
+    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, "PSPHOT.BACKMDL", binning->nXruff, binning->nYruff, PS_TYPE_F32);
+        model = pmFPAfileDefineInternal (config->files, name, binning->nXruff, binning->nYruff, PS_TYPE_F32);
     } else {
 	if (file->mode == PM_FPA_MODE_INTERNAL) {
@@ -123,5 +34,136 @@
 	}
     }
+
+    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
@@ -157,16 +199,18 @@
 		    modelData[iy][ix] = stats->robustMedian;
 		} else {	
-		    modelData[iy][ix] = psStatsGetValue(stats, statsOption);
+		    modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
 		}
+		modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
 	    } else {
 		psStatsOptions currentOptions = stats->options;
-		stats->options = PS_STAT_ROBUST_MEDIAN;
+		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] = NAN;
+		    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
@@ -183,15 +227,18 @@
     // patch over bad regions (use average of 8 possible neighbor pixels)
     // XXX consider testing all pixels against the 8 neighbors and replacing outliers...
-    float Count = 0;
-    float Value = 0;
+    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];
-		Count += 1;
+		ValueStdev += modelStdevData[iy][ix];
+		Count++;
 		continue;
 	    }
-	    float value = 0;
-	    float count = 0;
+
+	    double value = 0;
+	    double count = 0;
 	    for (int jy = iy - 1; jy <= iy + 1; jy++) {
 		if (jy <   0) continue;
@@ -206,9 +253,10 @@
 	    }
 	    if (count > 0) modelData[iy][ix] = value / count;
+	}
+    }
+    assert (Count > 0);
+    Value /= Count;
+    ValueStdev /= Count;
 	    
-	}
-    }
-    Value = Value / Count;
-
     // patch over remaining bad regions (use global average)
     for (int iy = 0; iy < model->image->numRows; iy++) {
@@ -216,4 +264,5 @@
 	    if (!isnan(modelData[iy][ix])) continue;
 	    modelData[iy][ix] = Value;
+	    modelStdevData[iy][ix] = ValueStdev;
 	}
     }
@@ -221,13 +270,26 @@
     psLogMsg ("psphot", PS_LOG_MINUTIA, "build median image: %f sec\n", psTimerMark ("psphot"));
 
-    // measure background stats and save for later output
-    psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_MIN | PS_STAT_MAX);
+    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_MEAN", PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKYSTDEV", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MAX",  PS_META_REPLACE, "sky model maximum value", statsBck->max);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MIN",  PS_META_REPLACE, "sky model minimum value", statsBck->min);
-    psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_NX",   PS_META_REPLACE, "sky model size (x)",      model->image->numCols);
-    psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_NY",   PS_META_REPLACE, "sky model size (y)",      model->image->numRows);
+    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);