Index: /trunk/psphot/src/psphotImageMedian.c
===================================================================
--- /trunk/psphot/src/psphotImageMedian.c	(revision 11310)
+++ /trunk/psphot/src/psphotImageMedian.c	(revision 11311)
@@ -1,3 +1,4 @@
 # include "psphot.h"
+# define TESTSAVE
 
 // generate the median in NxN boxes, clipping heavily
@@ -80,10 +81,23 @@
 
     // overhang : we will balance this evenly
-    int dx = (Nx % DX) / 2;
-    int dy = (Ny % DY) / 2;
+    int xExtra = (Nx % DX) / 2;
+    int yExtra = (Ny % DY) / 2;
+    int xOffset = (xExtra > 0) ? DX - xExtra : 0;
+    int yOffset = (yExtra > 0) ? DY - yExtra : 0;
+    xOffset -= image->col0;
+    yOffset -= image->row0;
 
     // dimensions of binned image
-    int nx = (Nx % DX) ? (int)(Nx / DX) + 1 : Nx / DX;
-    int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY;
+    int nx, ny;
+    if (Nx % DX == 0) {
+	nx = Nx / DX;
+    } else {
+	nx = (xExtra) ? (Nx / DX) + 2 : (Nx / DX) + 1;
+    }
+    if (Ny % DY == 0) {
+	ny = Ny / DY;
+    } else {
+	ny = (yExtra) ? (Ny / DY) + 2 : (Ny / DY) + 1;
+    }
 
     // select model pixels (from output background model file, or create internal file)
@@ -108,12 +122,13 @@
 
     assert(model->analysis != NULL);
-    psMetadataAdd(model->analysis, PS_LIST_TAIL, "XBIN", PS_DATA_S32 | PS_META_REPLACE,
+    psMetadataAdd(model->analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.XBIN", PS_DATA_S32 | PS_META_REPLACE,
 		  "Background x-binsize", DX);
-    psMetadataAdd(model->analysis, PS_LIST_TAIL, "YBIN", PS_DATA_S32 | PS_META_REPLACE,
+    psMetadataAdd(model->analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.YBIN", PS_DATA_S32 | PS_META_REPLACE,
 		  "Background x-binsize", DY);
-    psMetadataAdd(model->analysis, PS_LIST_TAIL, "x0", PS_DATA_S32 | PS_META_REPLACE,
-		  "Background x-overhang", dx);
-    psMetadataAdd(model->analysis, PS_LIST_TAIL, "y0", PS_DATA_S32 | PS_META_REPLACE,
-		  "Background y-overhang", dy);
+    psMetadataAdd(model->analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.XOFF", PS_DATA_S32 | PS_META_REPLACE,
+		  "Background x-overhang", xOffset);
+    psMetadataAdd(model->analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.YOFF", PS_DATA_S32 | PS_META_REPLACE,
+		  "Background y-overhang", yOffset);
+
 
     // measure clipped median for subimages
@@ -121,6 +136,6 @@
         for (int ix = 0; ix < nx; ix++) {
             // sx, sy are in parent coords
-            int sx = ix*DX - dx + image->col0;
-            int sy = iy*DY - dy + image->row0;
+            int sx = ix*DX - xOffset;
+            int sy = iy*DY - yOffset;
             region = psRegionSet (sx, sx + 2*DX, sy, sy + 2*DY);
             region = psRegionForImage (image, region);
@@ -133,5 +148,8 @@
             psImage *submask = psImageSubset (mask, region);
 
-            // XXX the value of the upper and lower cuts probably should be studied...
+            // 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
             if (psImageBackground(stats, subset, submask, 0xff, rng)) {
 		if (stats->options & PS_STAT_ROBUST_QUARTILE) {
@@ -151,16 +169,12 @@
 		    modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
 		}
-		psErrorClear(); // drop errors caused by this condition
+		// 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);
-
-            // XXX psImageStats is still very inefficient and poorly coded...
-            // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-            // stats = psImageStats (stats, subset, maskset, 0xff);
-            // backMdl->data.F32[iy][ix] = stats->clippedMean;
         }
     }
@@ -243,5 +257,9 @@
 
     // linear interpolation to full-scale
-    psImageUnbin (background->image, model->image, DX, DY, dx, dy);
+    if (!psImageUnbin (background->image, model->image, DX, DY, xOffset, yOffset)) {
+	psError (PSPHOT_ERR_PROG, true, "failed to build background iamge");
+	return false;
+    }
+	
     psLogMsg ("psphot", PS_LOG_MINUTIA, "build resampled image: %f sec\n", psTimerMark ("psphot"));
 
