Changeset 26868

Timestamp:

Feb 10, 2010, 4:08:49 PM (16 years ago)

Author:

eugene

Message:

updates from trunk (includes fix to robust stats)

Location:

branches/eam_branches/20091201/psLib/src

Files:

• fits/psFitsTable.c (modified) (1 diff)
• imageops/psImageBackground.c (modified) (4 diffs)
• imageops/psImageBackground.h (modified) (1 diff)
• math/psStats.c (modified) (5 diffs)

branches/eam_branches/20091201/psLib/src/fits/psFitsTable.c

@@ -162 +162 @@
 
         switch (typecode) {
+           // TBYTE and TSHORT fall though to read into psS32
           case TBYTE:
           case TSHORT:
-          case TLONGLONG:
             READ_TABLE_ROW_CASE(TLONG, long, S32, S32);
+            READ_TABLE_ROW_CASE(TLONGLONG, long, S64, S64);
             READ_TABLE_ROW_CASE(TFLOAT, float, F32, F32);
             READ_TABLE_ROW_CASE(TDOUBLE, double, F64, F64);

branches/eam_branches/20091201/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?
@@ -91 +98 @@
     }
     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 +117 @@
 
         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 +144 @@
                 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;

branches/eam_branches/20091201/psLib/src/imageops/psImageBackground.h

@@ -22 +22 @@
 #include <psRandom.h>
 
+void psImageBackgroundInit();
+
 // Get the background for an image
 bool psImageBackground(psStats *stats, // desired measurement and options

branches/eam_branches/20091201/psLib/src/math/psStats.c

@@ -749 +749 @@
     // Iterate to get the best bin size; an iteration limit is enforced at the bottom of the loop.
     for (int iterate = 1; iterate > 0; iterate++) {
-        psTrace(TRACE, 6,
-                "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n",
-                iterate);
+        psTrace(TRACE, 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
+
+        if (iterate >= PS_ROBUST_MAX_ITERATIONS) {
+          // This occurs when a large number of the values are identical --- a bin size cannot be found
+          // that will spread out the distribution.  Therefore, set what we can, and fall over
+          // gracefully.
+          COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
+          goto escape;
+        }
 
         // Get the minimum and maximum values
@@ -791 +797 @@
         psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize);
 
-        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can not specify the bin
-        // size precisely since the argument to psHistogramAlloc() is the number of bins, not the binSize.  If
-        // we get here, we know that binSize != 0.0.
-        long numBins = (max - min) / binSize; // Number of bins
+        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can
+        // not specify the bin size precisely since the argument to psHistogramAlloc() is
+        // the number of bins, not the binSize.  If we get here, we know that binSize !=
+        // 0.0.  We can also have a floating-point round-off error such that the last bin
+        // of the histogram does not correspond exactly with the value of 'max'.  Let's be
+        // a bit generous and extend the histogram by two bins in either direction
+        long numBins = 4 + (max - min) / binSize; // Number of bins
         psTrace(TRACE, 6, "Numbins is %ld\n", numBins);
         psTrace(TRACE, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
         // Generate the histogram
-        histogram = psHistogramAlloc(min, max, numBins);
+        histogram = psHistogramAlloc(min - 2.0*binSize, max + 2.0*binSize, numBins);
         // XXXXX we need to consider this step if errors -> variance
         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
@@ -807 +816 @@
             psFree(statsMinMax);
             psFree(mask);
-
             return false;
         }
@@ -814 +822 @@
             PS_VECTOR_PRINT_F32(histogram->nums);
         }
+
+        // perversity check: if most of the values land in a single bin, then we probably
+        // have a perverse case (eg, small number of points at extremely large / small
+        // values; nearly bi-modal distribution).  if so, keep only points within 5? 10?
+        // bins of that excess bin:
+        int nMaxBin = 0;
+        int iMaxBin = 0;
+        for (long i = 1; i < histogram->nums->n; i++) {
+            if (histogram->nums->data.F32[i] > nMaxBin) {
+                nMaxBin = histogram->nums->data.F32[i];
+                iMaxBin = i;
+            }
+        }
+        if (nMaxBin > numValid / 2) {
+            float minKeep = histogram->bounds->data.F32[iMaxBin] - 10*binSize;
+            float maxKeep = histogram->bounds->data.F32[iMaxBin + 1] + 10*binSize;
+            int nInvalid = 0;
+            for (long i = 0; i < myVector->n; i++) {
+                // skip the already-masked values
+                if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskVal) continue;
+                bool invalid = false;
+                invalid |= (myVector->data.F32[i] <= minKeep);
+                invalid |= (myVector->data.F32[i] >= maxKeep);
+                invalid |= (!isfinite(myVector->data.F32[i]));
+                if (!invalid) continue;
+                mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal;
+                nInvalid ++;
+            }
+
+            if (nInvalid) {
+              psTrace(TRACE, 6, "data is concentrated in a single bin, masking %d extreme outliers and retrying\n", nInvalid);
+              psFree(histogram);
+              psFree(cumulative);
+              histogram = NULL;
+              cumulative = NULL;
+              continue;
+            }
+            // if we did not mask anything, give up.
+        }
 
         // ADD step 1: Convert the specific histogram to a cumulative histogram
@@ -1007 +1054 @@
     stats->robustN50 = N50;
     psTrace(TRACE, 6, "The robustN50 is %ld.\n", N50);
+    psTrace(TRACE, 6, "The robust median and stdev are %f, %f\n", stats->robustMedian, stats->robustStdev);
 
     // Clean up