Index: trunk/psLib/src/imageops/psImageBackground.c
===================================================================
--- trunk/psLib/src/imageops/psImageBackground.c	(revision 27056)
+++ trunk/psLib/src/imageops/psImageBackground.c	(revision 27066)
@@ -46,5 +46,18 @@
     long ny = image->numRows;
 
-    const int Npixels = nx*ny;                // Total number of pixels
+    psImage *bad = psImageAlloc(nx, ny, PS_TYPE_U8); // Image with bad pixels
+    psImageInit(bad, 0);
+
+    int Npixels;                        // Total number of pixels
+    for (int y = 0; y < ny; y++) {
+        for (int x = 0; x < nx; x++) {
+            if (!isfinite(image->data.F32[iy][ix]) ||
+                (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue)) {
+                bad->data.U8[y][x] = 0xFF;
+            }
+            Npixels++;
+        }
+    }
+
     const int Nsubset = (stats->nSubsample == 0) ? Npixels : PS_MIN(stats->nSubsample, Npixels); // Number of pixels in subset
 
@@ -65,41 +78,48 @@
     long n = 0;                         // Number of actual pixels in subset
     if (Nsubset >= Npixels) {
-	// if we have an image smaller than Nsubset, just loop over the image pixels
-	for (int iy = 0; iy < ny; iy++) {
-	    for (int ix = 0; ix < nx; ix++) {
-		if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue)) {
-		    continue;
-		}
+        // if we have an image smaller than Nsubset, just loop over the (good) image pixels
+        for (int iy = 0; iy < ny; iy++) {
+            for (int ix = 0; ix < nx; ix++) {
+                if (bad->data.U8[iy][ix]) {
+                    continue;
+                }
 
-		float value = image->data.F32[iy][ix];
-		min = PS_MIN(value, min);
-		max = PS_MAX(value, max);
-		values->data.F32[n] = value;
-		n++;
-	    }
-	}
+                float value = image->data.F32[iy][ix];
+                min = PS_MIN(value, min);
+                max = PS_MAX(value, max);
+                values->data.F32[n] = value;
+                n++;
+            }
+        }
     } else {
-	for (long i = 0; i < Nsubset; i++) {
-	    double frnd = psRandomUniform(rng);
-	    int pixel = Npixels * frnd;
-	    int ix = pixel % nx;
-	    int iy = pixel / nx;
+        // Subsample all pixels
+        // This is not optimal, since there may be a large masked fraction that leaves us with few good pixels.
+        // In that case, you should have set Nsubset.......
+        Npixels = nx * ny;
+        for (long i = 0; i < Nsubset; i++) {
+            double frnd = psRandomUniform(rng);
+            int pixel = Npixels * frnd;
+            int ix = pixel % nx;
+            int iy = pixel / nx;
 
-	    if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue)) {
-		continue;
-	    }
+            if (bad->data.U8[iy][ix]) {
+                continue;
+            }
 
-	    float value = image->data.F32[iy][ix];
-	    min = PS_MIN(value, min);
-	    max = PS_MAX(value, max);
-	    values->data.F32[n] = value;
-	    n++;
-	}
+            float value = image->data.F32[iy][ix];
+            min = PS_MIN(value, min);
+            max = PS_MAX(value, max);
+            values->data.F32[n] = value;
+            n++;
+        }
     }
+
+    psFree(bad);
+
     if (n < 0.01*Nsubset) {
-	if ((nFailures < 3) || (nFailures % 100 == 0)) {
-	    psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld) (%d failures)", n, nFailures);
-	}
-	nFailures ++;
+        if ((nFailures < 3) || (nFailures % 100 == 0)) {
+            psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld) (%d failures)", n, nFailures);
+        }
+        nFailures ++;
         psFree(values);
         return false;
@@ -117,8 +137,8 @@
 
         if (!psVectorSort(values, values)) {
-	    if ((nFailures < 3) || (nFailures % 100 == 0)) {
-		psError(PS_ERR_UNKNOWN, false, "Unable to sort values.(%d failures)\n", nFailures);
-	    }
-	    nFailures ++;
+            if ((nFailures < 3) || (nFailures % 100 == 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to sort values.(%d failures)\n", nFailures);
+            }
+            nFailures ++;
             psFree(values);
             return false;
@@ -144,10 +164,10 @@
                 fclose (f);
             }
-	    if ((nFailures < 3) || (nFailures % 100 == 0)) {
-		psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background "
-			"(%dx%d, (row0,col0) = (%d,%d) (%d failures)\n",
-			image->numRows, image->numCols, image->row0, image->col0, nFailures);
-	    }		
-	    nFailures ++;
+            if ((nFailures < 3) || (nFailures % 100 == 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background "
+                        "(%dx%d, (row0,col0) = (%d,%d) (%d failures)\n",
+                        image->numRows, image->numCols, image->row0, image->col0, nFailures);
+            }
+            nFailures ++;
             psFree(values);
             return false;