Changeset 27840 for branches/simtest_nebulous_branches/psLib/src/imageops/psImageBackground.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psLib/src/imageops/psImageBackground.c (modified) (5 diffs)

branches/simtest_nebulous_branches/psLib/src/imageops/psImageBackground.c

@@ -15 +15 @@
 #include "psRandom.h"
 #include "psError.h"
+
+static int nFailures = 0;
+
+void psImageBackgroundInit() {
+    nFailures = 0;
+    return;
+}
 
 // XXX allow the user to choose the stats method?
@@ -39 +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 = 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
 
@@ -58 +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) {
-        psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL,
-                 "Unable to measure image background: too few data points (%ld)", n);
+        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;
@@ -108 +137 @@
 
         if (!psVectorSort(values, values)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to sort values.\n");
+            if ((nFailures < 3) || (nFailures % 100 == 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to sort values.(%d failures)\n", nFailures);
+            }
+            nFailures ++;
             psFree(values);
             return false;
@@ -132 +164 @@
                 fclose (f);
             }
-            psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background "
-                    "(%dx%d, (row0,col0) = (%d,%d)",
-                    image->numRows, image->numCols, image->row0, image->col0);
+            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;