Index: trunk/psphot/src/psphotImageMedian.c
===================================================================
--- trunk/psphot/src/psphotImageMedian.c	(revision 6324)
+++ trunk/psphot/src/psphotImageMedian.c	(revision 6324)
@@ -0,0 +1,77 @@
+# include "psphot.h"
+
+// generate the median in NxN boxes, clipping heavily
+// build a 2D clipped spline to the data (toss out the outliers?
+
+psImage *psphotImageMedian (pmReadout *readout, psMetadata *config, psStats *sky) 
+{ 
+
+    psRegion region;
+
+    // fit background to subset of image points within 3 sigma of sky
+    psTimerStart ("psphot");
+    psRandom *rnd     = psRandomAlloc (PS_RANDOM_TAUS, 0);
+
+    psImage *image = readout->image;
+    psImage *mask  = readout->mask;
+
+    bool      status  = false;
+    int       Nx = image->numCols + image->col0;
+    int       Ny = image->numRows + image->row0;
+    int       xbin = psMetadataLookupS32 (&status, config, "BACKGROUND_XBIN");
+    if (!status) {xbin = 64;}
+    int       ybin = psMetadataLookupS32 (&status, config, "BACKGROUND_YBIN");
+    if (!status) {ybin = 64;}
+
+    // we slide the box in half-steps from -0.5 off the left edge to +0.5 off the right edge 
+    int nx = 2 + 2 * Nx / xbin;
+    int ny = 2 + 2 * Ny / ybin;
+
+    psStats *stats = psStatsAlloc (PS_STATS_CLIPPED_MEAN, PS_STATS_CLIPPED_STDEV);
+    psImage *medimage = psImageAlloc (nx, ny, PS_TYPE_F32);
+
+    psVector *xval = psVectorAlloc (nx, PS_TYPE_F32);
+    psVector *yval = psVectorAlloc (ny, PS_TYPE_F32);
+
+    // measure clipped median for subimages 
+    for (int iy = -1; iy < ny - 1; iy++) {
+	for (int ix = -1; ix < nx - 1; ix++) {
+	    int sx = 0.5*ix*xbin;
+	    int sy = 0.5*iy*ybin;
+	    region = psRegionSet (sx, sx + xbin, sy, sy + ybin);
+	    region = psRegionForImage (readout->image, region);
+	    psImage *subset = psImageSubset (readout->image, region);
+	    stats = psImageStats (stats, subset);
+	    
+	    // XXX filter these on the clippedStdev?
+	    medimage->data.F32[iy][ix] = stats->clippedMean;
+
+	}
+    }
+    for (int ix = -1; ix < nx - 1; ix++) {
+	xval->data.F32[ix] = 0.5*ix*xbin + 0.5*xbin;
+    }
+    for (int iy = -1; iy < ny - 1; iy++) {
+	yval->data.F32[iy] = 0.5*iy*ybin + 0.5*ybin;
+    }
+
+    // x-dir / y-dir
+    psImage *devimage = psImageFitSpline (medimage, true);
+    psImage *background = psImageEvalSpline (medimage, devimage, xval, yval, Nx, Ny);
+
+    // this is a very inefficient way to evaluate the function..
+    // lame, temporary step to zero masked pixels
+    for (int j = 0; j < image->numRows; j++) {
+	for (int i = 0; i < image->numCols; i++) {
+	    if (!mask->data.U8[j][i]) {
+		image->data.F32[j][i] -= background->data.F32[j][i];
+	    }
+	}
+    }
+    psFree (medimage);
+    psFree (devimage);
+    psFree (rnd);  // XXX should this be made available at a higher level?
+
+    psLogMsg ("psphot", 3, "fit background model: %f sec\n", psTimerMark ("psphot"));
+    return (background);
+}