Changeset 30981

Timestamp:

Mar 18, 2011, 5:29:30 PM (15 years ago)

Author:

mwv

Message:

Fixes from MWV and Paul for randomn sky sampling.

File:

• branches/eam_branches/ipp-20110213/psLib/src/imageops/psImageBackground.c (modified) (3 diffs)

branches/eam_branches/ipp-20110213/psLib/src/imageops/psImageBackground.c

@@ -86 +86 @@
 
                 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 +94 @@
         }
     } 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; 
+        // 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;
@@ -108 +125 @@
 
             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);