Index: /trunk/psastro/src/psastroMosaicHeaders.c
===================================================================
--- /trunk/psastro/src/psastroMosaicHeaders.c	(revision 7384)
+++ /trunk/psastro/src/psastroMosaicHeaders.c	(revision 7385)
@@ -9,6 +9,6 @@
     psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO");
     if (!recipe) {
-	psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n");
-	exit(EXIT_FAILURE);
+        psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n");
+        exit(EXIT_FAILURE);
     }
 
@@ -16,6 +16,6 @@
     pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
     if (!input) {
-	psErrorStackPrint(stderr, "Can't find input data!\n");
-	exit(EXIT_FAILURE);
+        psErrorStackPrint(stderr, "Can't find input data!\n");
+        exit(EXIT_FAILURE);
     }
 
@@ -31,15 +31,15 @@
         psTrace (__func__, 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
         if (!chip->process || !chip->file_exists) { continue; }
-	
+
         // read WCS data from the corresponding header
         pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
 
-	pmAstromWriteBilevelChip (chip->toFPA, hdu->header, plateScale);
+        pmAstromWriteBilevelChip (chip->toFPA, hdu->header, plateScale);
     }
 
     psMetadata *mosaic = pmAstromWriteBilevelMosaic (fpa->projection, fpa->toTangentPlane, plateScale);
-    
+
     psFits *fits = psFitsOpen (mosastro, "w");
-    psFitsWriteHeaderNotImage (fits, mosaic);
+    psFitsWriteBlank(fits, mosaic);
     psFitsClose (fits);
 
Index: /trunk/psphot/src/Makefile.am
===================================================================
--- /trunk/psphot/src/Makefile.am	(revision 7384)
+++ /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 7384)
+++ /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