Changeset 27066 for trunk/psLib/src/imageops/psImageBackground.c

Timestamp:

Feb 24, 2010, 3:13:33 PM (16 years ago)

Author:

Paul Price

Message:

Better measurement of the background when the mask fraction is such that the number of good pixels is less than the subset size. Still not happy with the solution for when there's slightly more good pixels than the subset size, but the mask fraction is large. That probably requires a second pass.

File:

• trunk/psLib/src/imageops/psImageBackground.c (modified) (4 diffs)

trunk/psLib/src/imageops/psImageBackground.c

@@ -46 +46 @@
     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 +78 @@
     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 +137 @@
 
         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 +164 @@
                 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;