Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 17447) +++ trunk/psLib/src/math/psStats.c (revision 17612) @@ -13,6 +13,6 @@ * 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,6 +152,6 @@ // 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,8 +817,10 @@ // 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,11 +833,12 @@ 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,5 +914,5 @@ 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,5 +923,5 @@ // 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,8 +966,6 @@ // 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,5 +2589,5 @@ } - +# if (0) /****************************************************************************** fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general @@ -2696,4 +2697,132 @@ 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]);