Index: branches/eam_branches/ipp-20110213/psLib/src/imageops/psImageBackground.c
===================================================================
--- branches/eam_branches/ipp-20110213/psLib/src/imageops/psImageBackground.c	(revision 30991)
+++ branches/eam_branches/ipp-20110213/psLib/src/imageops/psImageBackground.c	(revision 31014)
@@ -100,34 +100,65 @@
         // 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.......
-        // 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; 
-        // Note that it's the number of good samples found, "n", that will generally terminate the loop
-        for (long i = 0; n < Nsubset && i < Npixels; i++) {
-            int pixel = frndSortedPixelOrder->data.S32[i];
-            int ix = pixel % nx;
-            int iy = pixel / nx;
-
-            if (bad->data.U8[iy][ix]) {
-                continue;
-            }
-
-            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);
-            values->data.F32[n] = value;
-            n++;
+
+        // 2011/03/22 - MWV:  On prompting from Gene, changed the sampling so that it doesn't build and then sort 
+        //     a huge array when Nsubset << Npixel.  
+        //   O(Npixel) log(Npixel) is unnecessarily painful when Npixel=38 million (one skycell) 
+        //     but we only needed Nsubset pixels and Nsubset << Npixel.
+
+        //   It also perhaps partially defeats any gain of taking only subsamples if we're sorting the array.
+        //   I should do some performance tests on this anyway.
+
+        // If our chance of collision is low and the gain from not sorting the array is likely to be high, 
+        // then just pick randomnly
+        // Note that we're guaranteeing Nsubset pixels returned here (up to the limit of Npixels)
+        if (Nsubset < 0.1 * Npixels) {
+            for (long i = 0; n < Nsubset && i < Npixels; i++) {
+                float frnd = psRandomUniform(rng);
+                int pixel = Npixels * frnd;
+                int ix = pixel % nx;
+                int iy = pixel / nx;
+
+                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++;
+            }
+        } else {
+            // 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
+            // Instead we take the following approach:
+            //   1) Calculate a random ordering of the pixels
+            //   2) Go through this ordering up to Nsubset to select pixels
+            // This is O(n log(n))
+            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; 
+            // Note that it's the number of good samples found, "n", that will generally terminate the loop
+            for (long i = 0; n < Nsubset && i < Npixels; i++) {
+                int pixel = frndSortedPixelOrder->data.S32[i];
+                int ix = pixel % nx;
+                int iy = pixel / nx;
+
+                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++;
+            }
         }
     }