Changeset 31029

Timestamp:

Mar 23, 2011, 4:22:31 PM (15 years ago)

Author:

eugene

Message:

use Fisher-Yates to select unique pixels in the subsample for appropriate sample sizes

File:

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

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

@@ -13 +13 @@
 #include "psType.h"
 #include "psAssert.h"
+#include "psAbort.h"
 #include "psRandom.h"
 #include "psError.h"
@@ -23 +24 @@
 }
 
-// XXX allow the user to choose the stats method?
-// (SAMPLE_MEAN, CLIPPED_MEAN, ROBUST_MEDIAN, FITTED_MEAN)
 bool psImageBackground(psStats *stats, psVector **sample, const psImage *image, const psImage *mask, psImageMaskType maskValue, psRandom *rng)
 {
@@ -45 +44 @@
     long nx = image->numCols;
     long ny = image->numRows;
-
-    psImage *bad = psImageAlloc(nx, ny, PS_TYPE_U8); // Image with bad pixels
-    psImageInit(bad, 0);
-
-    int Npixels = 0;                    // Total number of pixels
-    for (int y = 0; y < ny; y++) {
-        for (int x = 0; x < nx; x++) {
-            if (!isfinite(image->data.F32[y][x]) ||
-                (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 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
-
-    psVector *values;                   // Vector containing subsample
+    long nPixels = nx * ny;
+
+    long nSubset = (stats->nSubsample == 0) ? (nx * ny) : PS_MIN(stats->nSubsample, (nx * ny));
+
+    /*** discussion of the sampling and algorithm choices:
+
+         We want stats->nSubsample pixels.  How many good pixels do we actually have?  
+
+         There are a few domains of interest:
+
+         If nGoodPixels < f1 * nSubset, we should fail [1]
+         If nSubset >= nGoodPixels, we should just loop over all pixels [2]
+         If (nSubset <  f2 * nGoodPixels) and (nGoodPixels >= f3 * nPixels), we should just use the random sample technique [3]
+         If (nSubset >= f2 * nGoodPixels) or  (nGoodPixel  <  f3 * nPixels), we should use the Fisher-Yates technique.
+
+         to use Fisher-Yates, we need to have a copy of the pixels so we can re-shuffle.  We
+         should not do that with the whole list unless we need to.
+
+         [1] we just have too few samples to be useful; f1 ~ 0.01?
+
+         [2] we need to select from all pixels to reach the desired sample size
+
+         [3] this avoids a full-image-sized alloc, but only makes sense if nGoodPixels is actually
+         large.  f2 ~ 0.01, f3 ~ 0.1 (4 tries per success on average)
+
+    ***/
+
+    // count the number of good pixels
+    long nGoodPixels = 0;
+    for (long iy = 0; iy < ny; iy++) {
+        for (long ix = 0; ix < nx; ix++) {
+            if (!isfinite(image->data.F32[iy][ix])) continue;
+            if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue;
+            nGoodPixels ++;
+        }
+    }
+
+    // XXX This value of 1% is somewhat arbitrary
+    if (nGoodPixels < 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)", nGoodPixels, nFailures);
+        }
+        nFailures ++;
+        psTrace ("psLib.imageops", 4, "case 1: nGoodPixels < 0.01*nSubset (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels);
+        return false;
+    }
+
+    // alloc or recycle the vector containing subsample pixels
+    psVector *values;
     if (sample) {
-        *sample = psVectorRecycle(*sample, Nsubset, PS_TYPE_F32);
+        *sample = psVectorRecycle(*sample, nSubset, PS_TYPE_F32);
         values = psMemIncrRefCounter(*sample);
         values->n = 0;
     } else {
-        values = psVectorAllocEmpty(Nsubset, PS_TYPE_F32);
+        values = psVectorAllocEmpty(nSubset, PS_TYPE_F32);
     }
 
@@ -75 +105 @@
     float max = -PS_MAX_F32;
 
-    // select a subset of the image pixels to measure the stats
-    long n = 0;                         // Number of actual pixels in subset
-    if (Nsubset >= Npixels) {
-        // 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];
-                // 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++;
-            }
-        }
+    if ((nSubset < 0.01*nGoodPixels) && (nGoodPixels >= 0.1*nPixels)) {
+        psTrace ("psLib.imageops", 4, "case 3: nSubset < 0.01*nGoodPixels && nGoodPixels > 0.1*nPixels (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels);
+        // for small subsets, just use simple random sampling (saves a full-image-sized alloc),
+        // unless we have lots of masked pixels
+        for (long i = 0; i < nSubset; i++) {
+            double frnd = psRandomUniform(rng);
+            long pixel = nPixels * frnd;
+            long ix = pixel % nx;
+            long iy = pixel / nx;
+
+            if (!isfinite(image->data.F32[iy][ix])) continue;
+            if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue;
+
+            float value = image->data.F32[iy][ix];
+            min = PS_MIN(value, min);
+            max = PS_MAX(value, max);
+            psVectorAppend (values, value);
+        }
+    } else if (nSubset >= nGoodPixels) {
+            psTrace ("psLib.imageops", 4, "case 2: nSubset >= nGoodPixels (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels);
+        // in this case, we have to select from all masked pixels just to get the desired
+        // sample size
+        for (long iy = 0; iy < ny; iy++) {
+            for (long ix = 0; ix < nx; ix++) {
+                if (!isfinite(image->data.F32[iy][ix])) continue;
+                if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue;
+                float value = image->data.F32[iy][ix];
+                min = PS_MIN(value, min);
+                max = PS_MAX(value, max);
+                psVectorAppend(values, value);
+            }
+        }
     } 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.......
-
-        // 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++;
-            }
-        }
-    }
-
-    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 ++;
-        psFree(values);
-        return false;
-    }
-
-    values->n = n;
+        psTrace ("psLib.imageops", 4, "case 4: Fisher-Yates (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels);
+        // use Knuth's version of Fisher-Yates to select a random subset of the pixels, but
+        // only drawing each value once
+            
+        // generate a vector of all pixels which may in theory be selected
+        psVector *pixelVector = psVectorAllocEmpty(nGoodPixels, PS_TYPE_F32);
+        for (long iy = 0; iy < ny; iy++) {
+            for (long ix = 0; ix < nx; ix++) {
+                if (!isfinite(image->data.F32[iy][ix])) continue;
+                if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue;
+                psVectorAppend(pixelVector, image->data.F32[iy][ix]);
+            }
+        }
+        psAssert (nGoodPixels == pixelVector->n, "we must have mis-counted above");
+        
+        // generate the unique random subset
+        for (int i = 0; i < nSubset; i++) {
+            double frnd = psRandomUniform(rng);
+            int pixel = pixelVector->n * frnd;
+                
+            psVectorAppend(values, pixelVector->data.F32[pixel]);
+            pixelVector->data.F32[pixel] = pixelVector->data.F32[pixelVector->n - 1];
+            pixelVector->n --;
+        }
+        psFree (pixelVector);
+    }
+    psAssert (values->n >= 0.01*nSubset, "didn't we test this above?");
 
     if (stats->options & PS_STAT_ROBUST_QUARTILE) {
@@ -181 +170 @@
         // XXX this hack is just for testing, drop when I am happy with the psStats version of the values
 
-        int imin = stats->min * n;
-        int imax = stats->max * n;
+        int imin = stats->min * values->n;
+        int imax = stats->max * values->n;
         int npts = imax - imin + 1;
 
@@ -196 +185 @@
         // Subtract the median when we add the numbers, so we don't get numerical problems
         float median = npts % 2 ? 0.5 * (values->data.F32[npts/2 - 1] + values->data.F32[npts/2]) :
-                       values->data.F32[npts/2];
+            values->data.F32[npts/2];
         double value = 0;
-        for (long i = imin; (i <= imax) && (i < n); i++) {
+        for (long i = imin; (i <= imax) && (i < values->n); i++) {
             value += values->data.F32[i] - median;
         }
@@ -234 +223 @@
     return true;
 }
+
+
+/***  old comments
+
+ // 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.......
+
+ // 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)
+
+ // 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))
+
+ ***/