Index: trunk/psphot/src/psphotImageMedian.c
===================================================================
--- trunk/psphot/src/psphotImageMedian.c	(revision 6753)
+++ trunk/psphot/src/psphotImageMedian.c	(revision 6851)
@@ -1,11 +1,3 @@
 # 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
@@ -15,19 +7,24 @@
     bool status;
     psRegion region;
-    psImage *model = NULL;
+    int MAX_SAMPLE_PIXELS;
 
     psTimerStart ("psphot");
 
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
+
+    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
+    if (!status) MAX_SAMPLE_PIXELS = 1000;
+    psImageClippedStatsInit(MAX_SAMPLE_PIXELS);
+
+    // subtract this amount extra from the sky
+    float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS");
+    if (!status) SKY_BIAS = 0;
+
     // find the currently selected readout
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
-    pmFPAfile  *input   = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
-    pmReadout  *readout = pmFPAviewThisReadout (view, input->fpa);
+    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
 
     psImage *image = readout->image;
     psImage *mask  = readout->mask;
-
-    rnd = psRandomAlloc (PS_RANDOM_TAUS, 0);
-    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
-    if (!status) MAX_SAMPLE_PIXELS = 1000;
 
     // dimensions of input & output image
@@ -49,6 +46,14 @@
     int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY;
 
-    // select model pixels, from output background model file, or create
-    model = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
+    // select model pixels (from output background model file, or create internal file)
+    pmReadout *model = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKMDL");
+    if (model == NULL) {
+	// select model pixels, from output background model file, or create
+	model = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
+    } else {
+	// replace the supplied image data with an image of the desired size
+	psImageRecycle (model->image, nx, ny, PS_TYPE_F32);
+    }
+    psF32 **modelData = model->image->data.F32;
 
     // measure clipped median for subimages 
@@ -63,5 +68,6 @@
 	    psImage *submask = psImageSubset (mask, region);
 
-	    model->data.F32[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.50);
+	    // XXX the value of the upper and lower cuts probably should be studied...
+	    modelData[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.75) + SKY_BIAS;
 
 	    psFree (subset);
@@ -76,151 +82,36 @@
     psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot"));
 
-    // linear interpolation to full-scale
+    // XXX temporarily until we get the pmFPAfile output finished
+    psphotSaveImage (NULL, model->image, "model.fits");
 
     // select background pixels, from output background file, or create
-    psImage *background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", 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;
-		}
-	    }
-	}
+    pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
+    if (background == NULL) {
+	background = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
+    } else {
+	// replace the supplied image data with an image of the desired size
+	// XXX should not have to do this: it should be allocated correctly
+	// psImageRecycle (background->image, Nx, Ny, PS_TYPE_F32);
     }
 
-    // side pixels
-    int Xs = DX - dx;
-    int Xe = nx*DX - dx;
-    for (int Iy = 0; Iy < ny - 1; Iy++) {
+    psF32 **backData = background->image->data.F32;
 
-	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;
-	    }
-	}
-    }	
+    // linear interpolation to full-scale
+    psImageUnbin (background->image, model->image, DX, DY, dx, dy);
     psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot"));
 
-    psphotSaveImage (NULL, background, "back.fits");
+    // XXX temporarily until we get the pmFPAfile output finished...
+    psphotSaveImage (NULL, background->image, "back.fits");
 
     // back-sub image pixels, from output background file (don't create if not requested)
-    psImage *backSub = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKSUB", 0, 0, PS_TYPE_F32);
+    pmReadout *backSub = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKSUB");
 
-    // subtract the background model
+    // subtract the background model (save in backSub, if requested)
     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];
+		image->data.F32[j][i] -= backData[j][i];
 		if (backSub) {
-		    backSub->data.F32[j][i] = image->data.F32[j][i];
+		    backSub->image->data.F32[j][i] = image->data.F32[j][i];
 		}
 	    }
@@ -229,96 +120,7 @@
 
     psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
-    psFree (rnd);
+    psImageClippedStatsCleanup();
 
-    // it is safe to free all of these: either they are saved on config->files or are local
-    psFree (model);
-    psFree (background);
-    psFree (backSub);
+    // the pmReadout selected in this function are all view on entries in config->files
     return true;
 }
-
-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;
-    
-    // XXX correct for selection bias??
-
-    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;
-    }
-}
-