Index: trunk/psLib/src/imageops/psImageBackground.c
===================================================================
--- trunk/psLib/src/imageops/psImageBackground.c	(revision 10166)
+++ trunk/psLib/src/imageops/psImageBackground.c	(revision 10169)
@@ -11,7 +11,8 @@
 #include "psError.h"
 
-
+// XXX allow the user to choose the stats method?
+// (SAMPLE_MEAN, CLIPPED_MEAN, ROBUST_MEDIAN, FITTED_MEAN)
 psStats *psImageBackground(const psImage *image, const psImage *mask, psMaskType maskValue,
-                           double fmin, double fmax, long nMax, psRandom *rng)
+                           double fmin, double fmax, long nMax, psStatsOptions option, psRandom *rng)
 {
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
@@ -32,6 +33,4 @@
     long ny = image->numRows;
 
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE); // Statistics, for return
-
     int Nsubset = PS_MIN(nMax, nx*ny);  // Number of pixels in nubset
     int Npixels = nx*ny;                // Total number of pixels
@@ -43,4 +42,5 @@
     float max = values->data.F32[0];
 
+    // select a subset of the image pixels to measure the stats
     long n = 0;                         // Number of actual pixels in subset
     for (long i = 0; i < Nsubset; i++) {
@@ -62,30 +62,41 @@
     values->n = n;
 
-    int imin = fmin * n;
-    int imax = fmax * n;
-    int npts = imax - imin + 1;
+    psStats *stats = psStatsAlloc(option); // Statistics, for return
 
-    if (!psVectorSort(values, values)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to sort values.\n");
-        psFree(stats);
-        psFree(values);
-        return NULL;
+    if (stats->options & PS_STAT_ROBUST_QUARTILE) {
+        // use simple stats code (old verions)
+        // XXX this hack is just for testing, drop when I am happy with the psStats version of the values
+
+        int imin = fmin * n;
+        int imax = fmax * n;
+        int npts = imax - imin + 1;
+
+        if (!psVectorSort(values, values)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to sort values.\n");
+            psFree(stats);
+            psFree(values);
+            return NULL;
+        }
+
+        // Subtract the median when we add the numbers, so we don't get numerical problems
+        float median = npts % 2 ? 0.5 * (values->data.F32[npts/2 - 1] + values->data.F32[npts/2]) :
+                       values->data.F32[npts/2];
+        double value = 0;
+        for (long i = imin; (i <= imax) && (i < n); i++) {
+            value += values->data.F32[i] - median;
+        }
+        value = value / npts + median;
+        stats->robustMedian = value;
+        stats->robustUQ = values->data.F32[imax];
+        stats->robustLQ = values->data.F32[imin];
+    } else {
+        // XXX leave this as a psphot user option (passed in as part of stats?)
+        if (stats->options & PS_STAT_FITTED_MEAN) {
+            stats->clipSigma = 1.0;
+        }
+        psVectorStats (stats, values, NULL, NULL, 0);
     }
-
-    // Subtract the median when we add the numbers, so we don't get numerical problems
-    float median = npts % 2 ? 0.5 * (values->data.F32[npts/2 - 1] + values->data.F32[npts/2]) :
-                   values->data.F32[npts/2];
-    double value = 0;
-    for (long i = imin; (i <= imax) && (i < n); i++) {
-        value += values->data.F32[i] - median;
-    }
-    value = value / npts + median;
-
-    stats->robustMedian = value;
-    stats->robustUQ = values->data.F32[imax];
-    stats->robustLQ = values->data.F32[imin];
 
     psFree(values);
     return stats;
-    // XXX correct for selection bias??
 }
Index: trunk/psLib/src/imageops/psImageBackground.h
===================================================================
--- trunk/psLib/src/imageops/psImageBackground.h	(revision 10166)
+++ trunk/psLib/src/imageops/psImageBackground.h	(revision 10169)
@@ -13,4 +13,5 @@
                            double fmax, // Fraction to return in the upper quartile field (0.75 for LQ)
                            long nMax,   // Maximum number of pixels to subsample
+                           psStatsOptions option,
                            psRandom *rng // Random number generator (for pixel selection)
                           );