Changeset 17612

Timestamp:

May 9, 2008, 11:34:42 AM (18 years ago)

Author:

eugene

Message:

fixed error with cumulative histogram: values correspond to UPPER BOUNDS for input histogram

File:

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

trunk/psLib/src/math/psStats.c

@@ -13 +13 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.222 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-04-17 23:43:02 $
+ *  @version $Revision: 1.223 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-05-09 21:34:42 $
  *
  *  Copyright 2006 IfA, University of Hawaii
@@ -152 +152 @@
 // static prototypes:
 static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords);
-static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum,
-                                               psF32 yVal);
+// static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
+static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
 
 /******************************************************************************
@@ -817 +817 @@
 
         // ADD step 1: Convert the specific histogram to a cumulative histogram
+        // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
         cumulative = psHistogramAlloc(min, max, numBins);
         cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
         for (long i = 1; i < histogram->nums->n; i++) {
             cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
+            cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
         }
         if (psTraceGetLevel("psLib.math") >= 8) {
@@ -831 +833 @@
         psTrace(TRACE, 6, "Total data points is %ld\n", totalDataPoints);
 
+        // find bin which is the lower bound of median (value[binMedian] < median < value[binMedian+1]
         PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0);
+
         psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian,
                 cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);
 
         // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median.
-        stats->robustMedian = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
-                                                                binMedian, totalDataPoints/2.0);
+        stats->robustMedian = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
         if (isnan(stats->robustMedian)) {
             psError(PS_ERR_UNKNOWN, false,
@@ -911 +914 @@
         sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
 
-        psTrace(TRACE, 6, "The 2x sigma is %f.\n", sigma1);
+        psTrace(TRACE, 6, "The 1x sigma is %f.\n", sigma1);
         psTrace(TRACE, 6, "The 2x sigma is %f.\n", sigma2);
         psTrace(TRACE, 6, "The 4x sigma is %f.\n", sigma4);
@@ -920 +923 @@
         // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more
         // than 25 bins from the median, recalculate the bin size, and perform the algorithm again.
-        if (sigma < (2 * binSize)) {
+        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
@@ -963 +966 @@
     // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
     // positions.
-    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
-                                                         binLo25, totalDataPoints * 0.25f);
-    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
-                                                         binHi25, totalDataPoints * 0.75f);
+    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
+    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
     if (isnan(binLo25F32) || isnan(binHi25F32)) {
         psError(PS_ERR_UNKNOWN, false,
@@ -2588 +2589 @@
 }
 
-
+# if (0)
 /******************************************************************************
 fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general
@@ -2696 +2697 @@
         if (binNum == 0) {
             // We have two points only at the beginning of the vectors x and y.
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
+                              xVec->data.F32[binNum + 1]);
+        } else if (binNum == (xVec->n - 1)) {
+            // The special case where we have two points only at the end of
+            // the vectors x and y.
+            // XXX: Is this right?
+            tmpFloat = xVec->data.F32[binNum];
+        } else if (binNum == (xVec->n - 2)) {
+            // XXX: Is this right?
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]);
+        }
+    }
+
+    psTrace(TRACE, 6, "FIT: return %f\n", tmpFloat);
+    psFree(x);
+    psFree(y);
+
+    psTrace(TRACE, 5, "---- %s(%f) end ----\n", __func__, tmpFloat);
+    return tmpFloat;
+}
+# endif
+
+/******************************************************************************
+fitQuadraticSearchForYThenReturnXusingValues(*xVec, *yVec, binNum, yVal): A general routine
+which fits a quadratic to three points and returns the x-value corresponding to the input
+y-value.  This routine takes psVectors of x/y pairs as input, and fits a quadratic to the 3
+points surrounding element binNum in the vectors.  This version uses the values of x[i] for the
+x coordinates (not the midpoints).  This is appropriate for a cumulative histogram.  It then
+determines for what value x does that quadratic f(x) = yVal (the input parameter).
+
+*****************************************************************************/
+static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec,
+                                                          psVector *yVec,
+                                                          psS32 binNum,
+                                                          psF32 yVal
+    )
+{
+    psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
+    psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(xVec);
+        PS_VECTOR_PRINT_F32(yVec);
+    }
+
+    PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
+    PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
+    PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
+    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN);
+    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
+
+    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
+    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
+    psF32 tmpFloat = 0.0f;
+
+    if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) {
+        // The general case.  We have all three points.
+        x->data.F64[0] = xVec->data.F32[binNum - 1];
+        x->data.F64[1] = xVec->data.F32[binNum];
+        x->data.F64[2] = xVec->data.F32[binNum+1];
+        y->data.F64[0] = yVec->data.F32[binNum - 1];
+        y->data.F64[1] = yVec->data.F32[binNum];
+        y->data.F64[2] = yVec->data.F32[binNum + 1];
+        psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1],
+                xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
+        psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
+        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
+
+        //
+        // Ensure that the y value lies within range of the y values.
+        //
+        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
+               ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    _("Specified yVal, %g, is not within y-range, %g to %g."),
+                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
+        }
+
+        //
+        // Ensure that the y values are monotonic.
+        //
+        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
+            ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
+            psError(PS_ERR_UNKNOWN, true,
+                    "This routine must be called with monotonically increasing or decreasing data points.\n");
+            psFree(x);
+            psFree(y);
+            psTrace(TRACE, 5, "---- %s() end ----\n", __func__);
+            return NAN;
+        }
+
+        // Determine the coefficients of the polynomial.
+        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
+        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
+            psError(PS_ERR_UNEXPECTED_NULL, false,
+                    _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
+                      "Returning NAN."));
+            psFree(x);
+            psFree(y);
+            psTrace(TRACE, 5, "---- %s(NAN) end ----\n", __func__);
+            return NAN;
+        }
+        psTrace(TRACE, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
+        psTrace(TRACE, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
+        psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
+        psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n",
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
+
+        psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
+        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal,
+                                    x->data.F64[0], x->data.F64[2]);
+        psFree(myPoly);
+
+        if (isnan(tmpFloat)) {
+            psError(PS_ERR_UNEXPECTED_NULL,
+                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
+            psFree(x);
+            psFree(y);
+            psTrace(TRACE, 5, "---- %s(NAN) end ----\n", __func__);
+            return(NAN);
+        }
+    } else {
+        // These are special cases where the bin is at the beginning or end of the vector.
+        if (binNum == 0) {
+            // We have two points only at the beginning of the vectors x and y.
+            // XXX this does not seem to be doing the linear interpolation / extrapolation
             tmpFloat = 0.5 * (xVec->data.F32[binNum] +
                               xVec->data.F32[binNum + 1]);