Changeset 36121 for trunk/psLib/src/math/psStats.c

Timestamp:

Sep 17, 2013, 11:13:22 AM (13 years ago)

Author:

eugene

Message:

fix boundaries on cumulative histogram; iterative clip uses 25 fixed bins, not 25 cumulative hist bins; tweak linear fit to use cumulative vector mid points

File:

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

trunk/psLib/src/math/psStats.c

@@ -879 +879 @@
         cumulative = psHistogramAlloc(min, max, numBins);
         cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
-
+        cumulative->bounds->data.F32[0] = histogram->bounds->data.F32[0];
+
+        // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins).
+        long Nc = 1;  // track the current bin of cumulative
+        // the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the 
+        // upper end of the current histogram bin
+        for (long i = 1; i < histogram->nums->n; i++) {
+            if (histogram->nums->data.F32[i] == 0.0) continue;
+            cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i];
+            cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i];
+            Nc ++;
+        }
+        long Nlast = Nc - 1;  // last valid cumulative bin 
+        for (long i = Nc; i < histogram->nums->n; i++) { // Ensure the unused entries are filled.
+            cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast];
+            cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0;
+        }
+
+# if (0)
         // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins).
         long delta = 0;
@@ -898 +916 @@
           cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0;
         }
+# endif
         
         if (psTraceGetLevel("psLib.math") >= 8) {
@@ -1051 +1070 @@
         if (sigma < (3.0 * binSize)) {
             psTrace(TRACE, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
-            long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
-            long maskHi = PS_MIN(histogram->bounds->n - 1, (binMedian + 25)); // High index for masking
-            psF32 medianLo = histogram->bounds->data.F32[maskLo]; // Value at low index
-            psF32 medianHi = histogram->bounds->data.F32[maskHi]; // Value at high index
+
+            // these limits are supposed to be 25 x the raw bin size, NOT 25 of the cumulative histogram bins
+            psF32 medianLo = stats->robustMedian - 25*binSize;
+            psF32 medianHi = stats->robustMedian + 25*binSize;
+
+            // long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
+            // long maskHi = PS_MIN(cumulative->bounds->n - 1, (binMedian + 25)); // High index for masking
+            // psF32 medianLo = cumulative->bounds->data.F32[maskLo]; // Value at low index
+            // psF32 medianHi = cumulative->bounds->data.F32[maskHi]; // Value at high index
             psTrace(TRACE, 6, "Masking data more than 25 bins from the median\n");
-            psTrace(TRACE, 6,
-                    "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n",
-                    binMedian, maskLo, maskHi);
+            // psTrace(TRACE, 6, "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n", binMedian, maskLo, maskHi);
             psTrace(TRACE, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
+            int Nmasked = 0;
             for (long i = 0 ; i < myVector->n ; i++) {
                 if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
-                    // XXXX is this correct?  is MASK_MARK safe?
+                    if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK) continue;
                     mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK;
                     psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]);
+                    Nmasked ++;
                 }
             }
+
+            if (Nmasked == 0) {
+                // no significant change to the sigma & binsize -- we are done here 
+                iterate = -1;
+                continue;
+            }
 
             // Free the histograms; they will be recreated on the next iteration, with new bounds
@@ -1108 +1138 @@
     PS_BIN_FOR_VALUE (binLo25, cumulative->nums, totalDataPoints * 0.25f, 0);
     PS_BIN_FOR_VALUE (binHi25, cumulative->nums, totalDataPoints * 0.75f, 0);
-    psTrace(TRACE, 6, "The 25-percent and 75-precent data point bins are (%ld, %ld).\n", binLo25, binHi25);
+    psTrace(TRACE, 6, "The 25-percent and 75-percent data point bins are (%ld, %ld).\n", binLo25, binHi25);
 
     // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
@@ -2530 +2560 @@
                  yVec->data.F32[binNum - 0] +
                  yVec->data.F32[binNum + 1] + yVec->data.F32[binNum + 2]);
-    double Sy = (xVec->data.F32[binNum - 2] + xVec->data.F32[binNum - 1] +
-                 xVec->data.F32[binNum - 0] +
-                 xVec->data.F32[binNum + 1] + xVec->data.F32[binNum + 2]);
+
+    double YM = 0.5*(xVec->data.F32[binNum - 2] + xVec->data.F32[binNum - 1]);
+    double Ym = 0.5*(xVec->data.F32[binNum - 1] + xVec->data.F32[binNum - 0]);
+    double Yo = 0.5*(xVec->data.F32[binNum - 0] + xVec->data.F32[binNum + 1]);
+    double Yp = 0.5*(xVec->data.F32[binNum + 1] + xVec->data.F32[binNum + 2]);
+    double YP = 0.5*(xVec->data.F32[binNum + 2] + xVec->data.F32[binNum + 3]);
+
+    double Sy = YM + Ym + Yo + Yp + YP;
+
     double Sxx = (PS_SQR(yVec->data.F32[binNum - 2]) + PS_SQR(yVec->data.F32[binNum - 1]) +
                   PS_SQR(yVec->data.F32[binNum - 0]) +
                   PS_SQR(yVec->data.F32[binNum + 1]) + PS_SQR(yVec->data.F32[binNum + 2]));
-    double Sxy = (xVec->data.F32[binNum - 2] * yVec->data.F32[binNum - 2] +
-                  xVec->data.F32[binNum - 1] * yVec->data.F32[binNum - 1] +
-                  xVec->data.F32[binNum - 0] * yVec->data.F32[binNum - 0] +
-                  xVec->data.F32[binNum + 1] * yVec->data.F32[binNum + 1] +
-                  xVec->data.F32[binNum + 2] * yVec->data.F32[binNum + 2]);
-    double value = ((Sy * Sxx - Sx * Sxy) + (5 * Sxy - Sx * Sy) * yVal)/(5 * Sxx - Sx * Sx);
+
+    double Sxy = (yVec->data.F32[binNum - 2] * YM +
+                  yVec->data.F32[binNum - 1] * Ym +
+                  yVec->data.F32[binNum - 0] * Yo +
+                  yVec->data.F32[binNum + 1] * Yp +
+                  yVec->data.F32[binNum + 2] * YP);
+
+    double Det = 5*Sxx - Sx*Sx;
+    if (Det == 0.0) return NAN;
+
+    double C0 = (Sy*Sxx - Sx*Sxy) / Det;
+    double C1 = (Sxy*5 - Sx*Sy) / Det;
+
+    double value = C0 + yVal*C1;
 
     return value;