Index: /trunk/psphot/src/psphotImageMedian.c
===================================================================
--- /trunk/psphot/src/psphotImageMedian.c	(revision 10398)
+++ /trunk/psphot/src/psphotImageMedian.c	(revision 10399)
@@ -9,5 +9,4 @@
     pmFPAfile *file;
     psRegion region;
-    int MAX_SAMPLE_PIXELS;
     static char *defaultStatsName = "FITTED_MEAN";
 
@@ -22,8 +21,5 @@
     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
 
-    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
-    if (!status) {
-        MAX_SAMPLE_PIXELS = 1000;
-    }
+    // user supplied seed, if available
     unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED");
     if (!status) {
@@ -38,5 +34,5 @@
     }
 
-    // supply the sky background statistic name:
+    // supply the sky background statistics options
     char *statsName = psMetadataLookupStr (&status, recipe, "SKY_STAT");
     if (statsName == NULL) {
@@ -46,4 +42,23 @@
     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))) {
 	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;
     }
 
@@ -116,7 +131,5 @@
 
             // 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,
-					       statsOption, rng);
-	    if (stats != NULL) {
+            if (psImageBackground(stats, subset, submask, 0xff, rng)) {
 		if (stats->options & PS_STAT_ROBUST_QUARTILE) {
 		    modelData[iy][ix] = stats->robustMedian;
@@ -125,8 +138,7 @@
 		}
 	    } else {
-		psStatsOptions fallback_option = PS_STAT_ROBUST_MEDIAN;
-		stats = psImageBackground(subset, submask, 0xff, 0.25, 0.75, MAX_SAMPLE_PIXELS,
-					  fallback_option, rng);
-		if (stats == NULL) {
+		psStatsOptions currentOptions = stats->options;
+		stats->options |= PS_STAT_ROBUST_MEDIAN;
+		if (!psImageBackground(stats, subset, submask, 0xff, rng)) {
 		    psAbort(PS_FILE_LINE, "Failed to estimate background using ROBUST_MEDIAN for "
 			    "(%dx%d, (row0,col0) = (%d,%d)",
@@ -134,8 +146,8 @@
 		}
 		modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
+		stats->options = currentOptions;
 	    }
 	    modelData[iy][ix] += SKY_BIAS;
 
-            psFree (stats);
             psFree (subset);
             psFree (submask);
@@ -150,14 +162,14 @@
 
     // 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_STDEV", 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);
+    psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_MIN | PS_STAT_MAX);
+    psImageStats (statsBck, model->image, NULL, 0);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MEAN", PS_DATA_F32 | PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_STDEV", PS_DATA_F32 | PS_META_REPLACE, "sky model stdev",        statsBck->sampleStdev);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MAX",  PS_DATA_F32 | PS_META_REPLACE, "sky model maximum value", statsBck->max);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MIN",  PS_DATA_F32 | PS_META_REPLACE, "sky model minimum value", statsBck->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);
-    psLogMsg ("psphot", 3, "background sky : min %f mean %f max %f stdev %f", stats->min, stats->sampleMean, stats->max, stats->sampleStdev);
-    psFree (stats);
+    psLogMsg ("psphot", 3, "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
@@ -216,4 +228,5 @@
 
     psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
+    psFree(stats);
     psFree(rng);