Index: /branches/eam_branches/ipp-20101205/psLib/src/imageops/psImageBackground.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psLib/src/imageops/psImageBackground.c	(revision 30955)
+++ /branches/eam_branches/ipp-20101205/psLib/src/imageops/psImageBackground.c	(revision 30956)
@@ -86,4 +86,5 @@
 
                 float value = image->data.F32[iy][ix];
+                // 2011/03/16 - MWV:  Keep a running overall min and max
                 min = PS_MIN(value, min);
                 max = PS_MAX(value, max);
@@ -93,11 +94,32 @@
         }
     } else {
+        // 2011/03/16 - MWV: Hmmm... this overwites the previously defined Npixels that only counted finite-valued pixels.
+        Npixels = nx * ny; 
+
         // 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;
+        // 2011/03/16 - MWV:  Fixed double-sampling of sky pixels.
+        //   The pick-a-random-pixel-at-a-time approach  doesn't work well when Npixels is close to Nsubset
+        //   If Nsubset is 0.8*Npixels, then we will get lots of pixels double-counted
+        //   This is correct on average, but isn't the optimal way to estimate the sky level
+        // I suggest instead the following approach:
+        //   1) Calculate a random ordering of the pixels
+        //   2) Go through this ordering up to Nsubset to select pixels
+        psVector *frndPixelOrder = psVectorAlloc(Npixels, PS_TYPE_F32);
+        for (long i = 0; i < Npixels; i++) {  
+            frndPixelOrder->data.F32[i] = psRandomUniform(rng);  
+        }
+        // Now sort the array so that we end up with a list of the pixels in random order
+        psVector *frndSortedPixelOrder = psVectorSortIndex(NULL, frndPixelOrder);
+        // Now loop in our new sorted order
+        //  Paul Price suggests fixing this up so that it gets all Nsubset pixels, skipping over the masked ones.
+        //  This would be implemented with a while loop.  I'm leaving this comment-out till I test the simpler version that just makes this one conceptual change.
+        //  Ah, we can just use n
+        int i=0;
+        while ((n < Nsubset) && (i < Npixels)) {
+            i++;
+//        for (long i = 0; i < Nsubset; i++) {
+            int pixel = frndSortedPixelOrder->data.S32[i];
             int ix = pixel % nx;
             int iy = pixel / nx;
@@ -108,4 +130,5 @@
 
             float value = image->data.F32[iy][ix];
+            // 2011/03/16 - MWV:  Keep a running overall min and max
             min = PS_MIN(value, min);
             max = PS_MAX(value, max);