Index: trunk/psphot/src/image_stats.c
===================================================================
--- trunk/psphot/src/image_stats.c	(revision 4114)
+++ trunk/psphot/src/image_stats.c	(revision 4115)
@@ -13,22 +13,20 @@
     // pass on the stats
     { 
-	bool   status    = false;
+	bool     status  = false;
+	int      Npixels = image->numRows*image->numCols;
+	int      Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "NSUBSET"));
+	psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32);
 	psRandom *rnd    = psRandomAlloc (PS_RANDOM_TAUS, 0);
-	int    Nsubset   = psMetadataLookupF32 (&status, config, "NSUBSET");
-	double fSubset   = (double) Nsubset / (image->numRows*image->numCols);
-	psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32);
 
-	subset->n = 0;
-	for (int i = 0; i < image->numRows; i++) {
-	    for (int j = 0; j < image->numCols; j++) {
-		double frnd = psRandomUniform (rnd);
-		if ((fSubset < 1.0) && (fSubset < frnd)) continue;
-		subset->data.F32[subset->n] = image->data.F32[j][i];
-		subset->n ++;
-		if (subset->n == Nsubset) goto got_subset;
-	    }
+	// choose Nsubset points between 0 and Nx*Ny, convert to coords
+	Npixels = image->numRows*image->numCols;
+	for (int i = 0; i < Nsubset; i++) {
+	  double frnd = psRandomUniform (rnd);
+	  int pixel = Npixels * frnd;
+	  int ix = pixel / image->numCols;
+	  int iy = pixel % image->numCols;
+	  subset->data.F32[i] = image->data.F32[iy][ix];
 	}
 
-    got_subset:
 	// this should use ROBUST not SAMPLE median
 	// robust median is broken in pslib (0.5)