Index: trunk/psphot/src/Makefile.am
===================================================================
--- trunk/psphot/src/Makefile.am	(revision 7350)
+++ trunk/psphot/src/Makefile.am	(revision 7385)
@@ -7,4 +7,6 @@
 psphot_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(psphot_CFLAGS)
 psphot_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS)
+
+PSLIB_DIR = `pslib-config --prefix`
 
 # psphot_CPPFLAGS = $(libpsphot_la_CPPFLAGS)
Index: trunk/psphot/src/psphotImageMedian.c
===================================================================
--- trunk/psphot/src/psphotImageMedian.c	(revision 7350)
+++ trunk/psphot/src/psphotImageMedian.c	(revision 7385)
@@ -3,6 +3,6 @@
 // generate the median in NxN boxes, clipping heavily
 // linear interpolation to generate full-scale model
-bool psphotImageMedian (pmConfig *config, pmFPAview *view) 
-{ 
+bool psphotImageMedian (pmConfig *config, pmFPAview *view)
+{
     bool status;
     psRegion region;
@@ -15,12 +15,18 @@
 
     MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
-    if (!status) MAX_SAMPLE_PIXELS = 1000;
+    if (!status) {
+        MAX_SAMPLE_PIXELS = 1000;
+    }
     unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED");
-    if (!status) seed = 0;
-    psImageClippedStatsInit(MAX_SAMPLE_PIXELS, seed);
+    if (!status) {
+        seed = 0;
+    }
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, seed);
 
     // subtract this amount extra from the sky
     float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS");
-    if (!status) SKY_BIAS = 0;
+    if (!status) {
+        SKY_BIAS = 0;
+    }
 
     // find the currently selected readout
@@ -49,36 +55,36 @@
     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);
+        // 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);
+        // 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 
+    // measure clipped median for subimages
     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);
+        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);
 
-	    // XXX the value of the upper and lower cuts probably should be studied...
-	    psStats *stats = psImageClippedStats (subset, submask, 0xff, 0.25, 0.75);
-	    modelData[iy][ix] = stats->robustMedian + SKY_BIAS;
+            // XXX the value of the upper and lower cuts probably should be studied...
+            psStats *stats = psImageBackground(subset, submask, 0xff, 0.25, 0.75, MAX_SAMPLE_PIXELS, rng);
+            modelData[iy][ix] = stats->robustMedian + SKY_BIAS;
 
-	    psFree (stats);
-	    psFree (subset);
-	    psFree (submask);
+            psFree (stats);
+            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;
-	}
+            // 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;
+        }
     }
     psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot"));
@@ -98,6 +104,6 @@
     pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
     if (background == NULL) {
-	background = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
-    } 
+        background = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
+    }
     psF32 **backData = background->image->data.F32;
 
@@ -115,16 +121,16 @@
     // 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] -= backData[j][i];
-		if (backSub) {
-		    backSub->image->data.F32[j][i] = image->data.F32[j][i];
-		}
-	    }
-	}
+        for (int i = 0; i < image->numCols; i++) {
+            if (!mask->data.U8[j][i]) {
+                image->data.F32[j][i] -= backData[j][i];
+                if (backSub) {
+                    backSub->image->data.F32[j][i] = image->data.F32[j][i];
+                }
+            }
+        }
     }
 
     psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
-    psImageClippedStatsCleanup();
+    psFree(rng);
 
     // the pmReadout selected in this function are all view on entries in config->files