Index: trunk/psphot/src/psphotImageMedian.c
===================================================================
--- trunk/psphot/src/psphotImageMedian.c	(revision 6324)
+++ trunk/psphot/src/psphotImageMedian.c	(revision 6399)
@@ -1,65 +1,208 @@
 # include "psphot.h"
+double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax);
+void fsort (float *value, int N);
+
+// random number seed to select a fraction of the image pixels
+static psRandom *rnd;
+
+// use no more than MAX_SAMPLE_PIXELS pixels for each median box
+static int MAX_SAMPLE_PIXELS;
 
 // generate the median in NxN boxes, clipping heavily
-// build a 2D clipped spline to the data (toss out the outliers?
-
+// linear interpolation to generate full-scale model
 psImage *psphotImageMedian (pmReadout *readout, psMetadata *config, psStats *sky) 
 { 
-
+    bool status;
     psRegion region;
+
+    rnd = psRandomAlloc (PS_RANDOM_TAUS, 0);
+    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, config, "IMSTATS_NPIX");
+    if (!status) MAX_SAMPLE_PIXELS = 1000;
 
     // 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);
+    // dimensions of input & output image
+    int Nx = image->numCols;
+    int Ny = image->numRows;
+
+    // scaling factor
+    int DX = psMetadataLookupS32 (&status, config, "BACKGROUND_XBIN");
+    if (!status) {DX = 64;}
+    int DY = psMetadataLookupS32 (&status, config, "BACKGROUND_YBIN");
+    if (!status) {DY = 64;}
+
+    // overhang : we will balance this evenly
+    int dx = (Nx % DX) / 2;
+    int dy = (Ny % DY) / 2;
+
+    // dimensions of binned image
+    int nx = (Nx % DX) ? (int)(Nx / DX) + 1 : Nx / DX;
+    int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY;
+
+    // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+    psImage *model = psImageAlloc (nx, 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 iy = 0; iy < ny; iy++) {
+	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;
+	    region = psRegionSet (sx, sx + 2*DX, sy, sy + 2*DY);
+	    region = psRegionForImage (image, region);
+	    psImage *subset  = psImageSubset (image, region);
+	    psImage *submask = psImageSubset (mask, region);
+
+	    model->data.F32[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.50);
+
+	    psFree (subset);
+	    psFree (submask);
+
+	    // XXX psImageStats is still very inefficient and poorly coded...
+	    // stats = psImageStats (stats, subset, maskset, 0xff);
+	    // model->data.F32[iy][ix] = stats->clippedMean;
+	}
+    }
+    psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot"));
+
+    // linear interpolation to full-scale
+    
+    psImage *background = psImageAlloc (Nx, Ny, PS_TYPE_F32);
+
+    // XXX this code skips the initial pixels
+    for (int Iy = 0; Iy < ny-1; Iy ++) {
+	for (int Ix = 0; Ix < nx-1; Ix ++) {
+
+	    float V00 = model->data.F32[Iy+0][Ix+0];
+	    float V01 = model->data.F32[Iy+0][Ix+1];
+	    float V10 = model->data.F32[Iy+1][Ix+0];
+	    float V11 = model->data.F32[Iy+1][Ix+1];
+
+	    // a single binned pixel quad
+	    // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords
+	    // corresponding to (Ix,Iy), (Ix+1,Iy+1)
+	    int Xs = (Ix + 1)*DX - dx;
+	    int Ys = (Iy + 1)*DY - dy;
+	    int Xe = Xs + DX;
+	    int Ye = Ys + DY;
+
+	    for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
+		float Vxs = (V10 - V00)*(iy - Ys) / DY + V00;
+		float Vxe = (V11 - V01)*(iy - Ys) / DY + V01;
+		float dV = (Vxe - Vxs) / DX;
+		float V  = Vxs;
+		for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
+		    background->data.F32[iy][ix] = V;
+		    V += dV;
+		}
+	    }
+	}
+    }
+
+    // side pixels
+    int Xs = DX - dx;
+    int Xe = nx*DX - dx;
+    for (int Iy = 0; Iy < ny - 1; Iy++) {
+
+	int Ys = (Iy + 1)*DY - dy;
+	int Ye = Ys + DY;
+
+	// leading edge
+	float V0 = model->data.F32[Iy+0][0];
+	float V1 = model->data.F32[Iy+1][0];
+	float dV = (V1 - V0) / DY;
+	float V = V0;
+	for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
+	    for (int ix = 0; ix < Xs; ix++) {
+		background->data.F32[iy][ix] = V;
+	    }
+	    V += dV;
+	}
+
+	// trailing edge
+	V0 = model->data.F32[Iy+0][nx-1];
+	V1 = model->data.F32[Iy+1][nx-1];
+	dV = (V1 - V0) / DY;
+	V = V0;
+	for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
+	    for (int ix = Xe; ix < Nx; ix++) {
+		background->data.F32[iy][ix] = V;
+	    }
+	    V += dV;
+	}
+    }
+
+    // top and bottom pixels
+    int Ys = DY - dy;
+    int Ye = ny*DY - dy;
+    for (int Ix = 0; Ix < nx - 1; Ix++) {
+
+	int Xs = (Ix + 1)*DX - dx;
+	int Xe = Xs + DX;
+
+	// top edge
+	float V0 = model->data.F32[0][Ix+0];
+	float V1 = model->data.F32[0][Ix+1];
+	float dV = (V1 - V0) / DX;
+	for (int iy = 0; iy < Ys; iy++) {
+	    float V = V0;
+	    for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
+		background->data.F32[iy][ix] = V;
+		V += dV;
+	    }
+	}
+
+	// bottom edge
+	V0 = model->data.F32[ny-1][Ix+0];
+	V1 = model->data.F32[ny-1][Ix+1];
+	dV = (V1 - V0) / DX;
+	for (int iy = Ye; iy < Ny; iy++) {
+	    float V = V0;
+	    for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
+		background->data.F32[iy][ix] = V;
+		V += dV;
+	    }
+	}
+    }
+
+    // the four corners
+    {
+	float V;
+	// 0,0
+	V = model->data.F32[0][0];
+	for (int iy = 0; iy < DY - dy; iy++) {
+	    for (int ix = 0; ix < DX - dx; ix++) {
+		background->data.F32[iy][ix] = V;
+	    }
+	}
+	// Nx,0
+	V = model->data.F32[0][nx-1];
+	for (int iy = 0; iy < DY - dy; iy++) {
+	    for (int ix = nx*DX - dx; ix < Nx; ix++) {
+		background->data.F32[iy][ix] = V;
+	    }
+	}
+	// 0,Ny
+	V = model->data.F32[ny-1][0];
+	for (int iy = ny*DY - dy; iy < Ny; iy++) {
+	    for (int ix = 0; ix < DX - dx; ix++) {
+		background->data.F32[iy][ix] = V;
+	    }
+	}
+	// Nx,Ny
+	V = model->data.F32[nx-1][ny-1];
+	for (int iy = ny*DY - dy; iy < Ny; iy++) {
+	    for (int ix = nx*DX - dx; ix < Nx; ix++) {
+		background->data.F32[iy][ix] = V;
+	    }
+	}
+    }	
+    psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot"));
+
+    // subtract the background model
     for (int j = 0; j < image->numRows; j++) {
 	for (int i = 0; i < image->numCols; i++) {
@@ -69,9 +212,90 @@
 	}
     }
-    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"));
+    psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
+    psFree (rnd);
     return (background);
 }
+
+double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax) {
+
+    double value;
+    int nx = image->numCols;
+    int ny = image->numRows;
+    
+    if (nx*ny <= 0) return 0.0;
+
+    int Nsubset = PS_MIN (MAX_SAMPLE_PIXELS, nx*ny);
+    int Npixels = nx*ny;
+
+    psVector *values = psVectorAlloc (Nsubset, PS_TYPE_F32);
+
+    float min = values->data.F32[0];
+    float max = values->data.F32[0];
+
+    int n = 0;
+    for (int i = 0; i < Nsubset; i++) {
+	double frnd = psRandomUniform (rnd);
+	int pixel = Npixels * frnd;
+	int ix = pixel % nx;
+	int iy = pixel / nx;
+
+	if (mask->data.U8[iy][ix] & maskValue) continue;
+
+	value = image->data.F32[iy][ix];
+	min = PS_MIN (value, min);
+	max = PS_MIN (value, max);
+	values->data.F32[n] = value;
+	n++;
+    }
+    values->n = n;
+
+    int imin = fmin * n;
+    int imax = fmax * n;
+    int npts = imax - imin + 1;
+
+    fsort (values->data.F32, n);
+
+    value = 0;
+    for (int i = imin; (i <= imax) && (i < n); i++) {
+	value += values->data.F32[i];
+    }
+    value = value / npts;
+    
+    psFree (values);
+    return value;
+}
+
+void fsort (float *value, int N) {
+
+    int l,j,ir,i;
+    float temp;
+  
+    if (N < 2) return;
+    l = N >> 1;
+    ir = N - 1;
+    for (;;) {
+	if (l > 0) {
+	    temp = value[--l];
+	}
+	else {
+	    temp = value[ir];
+	    value[ir] = value[0];
+	    if (--ir == 0) {
+		value[0] = temp;
+		return;
+	    }
+	}
+	i = l;
+	j = (l << 1) + 1;
+	while (j <= ir) {
+	    if (j < ir && value[j] < value[j+1]) ++j;
+	    if (temp < value[j]) {
+		value[i]=value[j];
+		j += (i=j) + 1;
+	    }
+	    else j = ir + 1;
+	}
+	value[i] = temp;
+    }
+}
+