Changeset 6270

Timestamp:

Jan 31, 2006, 4:07:24 PM (20 years ago)

Author:

gusciora

Message:

Fixed bug 674.

File:

• trunk/psLib/src/math/psStats.c (modified) (13 diffs)

trunk/psLib/src/math/psStats.c

@@ -22 +22 @@
  *
  *
- *  @version $Revision: 1.163 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-26 23:49:11 $
+ *  @version $Revision: 1.164 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-01 02:07:24 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1644 +1644 @@
 }
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 /******************************************************************************
 p_psVectorRobustStats(myVector, maskVector, maskVal, stats): This is the new
@@ -1662 +1754 @@
 XXX: Review and ensure that all memory is free'ed at premature exits.
 *****************************************************************************/
+#define INITIAL_NUM_BINS 1000.0
 psS32 p_psVectorRobustStats(const psVector* myVector,
                             const psVector* errors,
@@ -1697 +1790 @@
 
     while (iterate > 0) {
-        psTrace(__func__, 6, "Iterating on Bin size.\n");
+        psTrace(__func__, 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
         psF32 binSize = 0.0;
         if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) {
@@ -1717 +1810 @@
                 return(1);
             }
-            psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
-            binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f;
+            psTrace(__func__, 6, "Data min/max is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
+            binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / INITIAL_NUM_BINS;
             psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize);
         }
@@ -1756 +1849 @@
         numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize);
         psTrace(__func__, 6, "Numbins is %d\n", numBins);
+        psTrace(__func__, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
         robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
         cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
@@ -1792 +1886 @@
         } else {
             tmpScalar->data.F32 = totalDataPoints/2.0;
-            binMedian = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
             if (binMedian < 0) {
                 psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian);
@@ -1804 +1898 @@
             }
         }
-        psTrace(__func__, 6, "The median bin is %d\n", binMedian);
+        psTrace(__func__, 6, "The median bin is %d (%.2f to %.2f)\n", binMedian, cumulativeRobustHistogram->bounds->data.F32[binMedian], cumulativeRobustHistogram->bounds->data.F32[binMedian+1]);
 
         //
@@ -1839 +1933 @@
         } else {
             // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
-            binHi = 1+ p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
             if (binHi > cumulativeRobustHistogram->nums->n-1) {
                 binHi = cumulativeRobustHistogram->nums->n-1;
             }
         }
-        psTrace(__func__, 6, "The 15.8655% and 84.1345% data point bins are (%d, %d).\n", binLo, binHi);
+        psTrace(__func__, 6, "The 15.8655-percent and 84.1345-percent data point bins are (%d, %d).\n", binLo, binHi);
         psTrace(__func__, 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binLo));
         psTrace(__func__, 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binHi));
@@ -1920 +2014 @@
         if (sigma < (2 * binSize)) {
             psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
-            psF32 medianLo = robustHistogram->bounds->data.F32[binMedian - 25];
-            psF32 medianHi = robustHistogram->bounds->data.F32[binMedian + 25];
+            psF32 medianLo;
+            psF32 medianHi;
+            if ((binMedian - 25) > 0) {
+                medianLo = robustHistogram->bounds->data.F32[binMedian - 25];
+            } else {
+                medianLo = robustHistogram->bounds->data.F32[0];
+            }
+            if ((binMedian + 25) < robustHistogram->bounds->n) {
+                medianHi = robustHistogram->bounds->data.F32[binMedian + 25];
+            } else {
+                medianHi = robustHistogram->bounds->data.F32[robustHistogram->bounds->n - 1];
+            }
+
+            psTrace(__func__, 6, "Masking data more than 25 bins from the median (%.2f)\n", stats->robustMedian);
+            psTrace(__func__, 6, "The median is at bin number %d.  We mask bins outside the bin range (%d:%d)\n", binMedian, binMedian - 25, binMedian + 25);
             psTrace(__func__, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
             for (psS32 i = 0 ; i < myVector->n ; i++) {
-                if ((myVector->data.F32[i] < medianLo) ||
-                        (myVector->data.F32[i] > medianHi)) {
+                if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
                     tmpMaskVec->data.U8[i] = 1;
-                    psTrace(__func__, 6, "Masking element %d is %f\n", i, myVector->data.F32[i]);
+                    psTrace(__func__, 8, "Masking element %d is %f\n", i, myVector->data.F32[i]);
                 }
             }
@@ -1951 +2057 @@
         binLo25 = 0;
     } else {
-        binLo25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+        binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
     }
     tmpScalar->data.F32 = totalDataPoints * 0.75f;
@@ -1958 +2064 @@
         binHi25 = 0;
     } else {
-        binHi25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+        binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
     }
     if ((binLo25 < 0) || (binHi25 < 0)) {
@@ -1970 +2076 @@
         return(1);
     }
-    psTrace(__func__, 6, "The 25% and 75% data point bins are (%d, %d).\n", binLo25, binHi25);
+    psTrace(__func__, 6, "The 25-percent and 75-precent data point bins are (%d, %d).\n", binLo25, binHi25);
 
     //