Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 6314) +++ /trunk/psLib/src/math/psStats.c (revision 6315) @@ -16,6 +16,6 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.165 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-02 21:09:08 $ + * @version $Revision: 1.166 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-03 01:30:14 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -878,6 +878,4 @@ psF32 sumSquares = 0.0; // temporary variable psF32 sumDiffs = 0.0; // temporary variable - // psF32 sum1; - // psF32 sum2; psF32 errorDivisor = 0.0f; @@ -913,6 +911,5 @@ } else { for (i = 0; i < myVector->n; i++) { - if ((stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { + if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { diff = myVector->data.F32[i] - mean; sumSquares += (diff * diff); @@ -965,5 +962,10 @@ // data ranges are used correctly. if (errors != NULL) { - stats->sampleStdev = (1.0 / sqrtf(errorDivisor)); + if (0) { + stats->sampleStdev = (1.0 / sqrtf(errorDivisor)); + } else { + countFloat = (psF32)countInt; + stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); + } } else { countFloat = (psF32)countInt; @@ -998,4 +1000,5 @@ { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + psTrace(__func__, 4, "Trace level is %d\n", psTraceGetLevel(__func__)); psF32 clippedMean = 0.0; // self-explanatory psF32 clippedStdev = 0.0; // self-explanatory @@ -1071,4 +1074,5 @@ return(-2); } + psTrace(__func__, 6, "The initial sample median is %f\n", statsTmp->sampleMedian); // 2. Compute the sample standard deviation. @@ -1081,4 +1085,5 @@ return(-2); } + psTrace(__func__, 6, "The initial sample stdev is %f\n", statsTmp->sampleStdev); // 3. Use the sample median as the first estimator of the mean X. @@ -1090,9 +1095,9 @@ // 5. Repeat N (stats->clipIter) times: for (psS32 iter = 0; iter < stats->clipIter; iter++) { + psTrace(__func__, 6, "------------ Iteration %d ------------\n", iter); // a) Exclude all values x_i for which |x_i - x| > K * stdev if (errors != NULL) { for (psS32 j = 0; j < myVector->n; j++) { - if (fabs(myVector->data.F32[j] - clippedMean) > - (stats->clipSigma * errors->data.F32[j])) { + if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * errors->data.F32[j])) { tmpMask->data.U8[j] = 0xff; } @@ -1110,4 +1115,6 @@ p_psVectorSampleMean(myVector, errors, tmpMask, 0xff, statsTmp); p_psVectorSampleStdev(myVector, errors, tmpMask, 0xff, statsTmp); + psTrace(__func__, 6, "The new sample mean is %f\n", statsTmp->sampleMean); + psTrace(__func__, 6, "The new sample stdev is %f\n", statsTmp->sampleStdev); // If the new mean and stdev are NAN, we must exit the loop. @@ -1128,8 +1135,10 @@ if (stats->options & PS_STAT_CLIPPED_MEAN) { stats->clippedMean = clippedMean; + psTrace(__func__, 6, "The final clipped mean is %f\n", clippedMean); } // 8. The last calcuated value of stdev is the cliped stdev. if (stats->options & PS_STAT_CLIPPED_STDEV) { stats->clippedStdev = clippedStdev; + psTrace(__func__, 6, "The final clipped stdev is %f\n", clippedStdev); } @@ -1139,113 +1148,11 @@ } - - - -/****************************************************************************** -p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, getThisValue): This routine -takes as input a 1-D polynomial of arbitrary order and a range of x-values for -which it is defined: [rangeLow, rangeHigh]. It determines the x-value of -that polynomial such that f(x) == getThisValue. This function uses a -binary-search algorithm on the range and assumes that the polynomial is -monotonically increasing or decreasing within that range. - -XXX: Terminate when f(x)-getThisValue is within some error tolerance. - -XXX: Create a 2nd-order polynomial version and solve for X analytically. - *****************************************************************************/ -psF32 p_ps1DPolyMedian( - psPolynomial1D* myPoly, - psF32 rangeLow, - psF32 rangeHigh, - psF32 getThisValue) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psTrace(__func__, 4, "range (%f, %f) get (%f)\n", rangeLow, rangeHigh, getThisValue); - PS_ASSERT_POLY_NON_NULL(myPoly, NAN); - PS_ASSERT_FLOAT_LARGER_THAN(rangeHigh, rangeLow, NAN); - // We ensure that the requested f(y) value, which is getThisValue, is - // falls within the range of y-values of the polynomial "myPoly" in the - // specified x-range (rangeLow:rangeHigh). - psF32 fLo = psPolynomial1DEval(myPoly, rangeLow); - psF32 fHi = psPolynomial1DEval(myPoly, rangeHigh); - psTrace(__func__, 4, "function at endpoints are (%f, %f) get (%f)\n", fLo, fHi, getThisValue); - if (!((fLo <= getThisValue) && (fHi >= getThisValue))) { - psError(PS_ERR_UNKNOWN, true, "The requested y value (%f) does not fall within the specified range (%f to %f)\n", getThisValue, fLo, fHi); - psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); - return(NAN); - } - - psS32 numIterations = 0; - psF32 midpoint = 0.0; - psF32 oldMidpoint = 1.0; - psF32 f = 0.0; - - while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) { - midpoint = (rangeHigh + rangeLow) / 2.0; - if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) { - return (midpoint); - } - oldMidpoint = midpoint; - - f = psPolynomial1DEval(myPoly, midpoint); - printf("p_ps1DPolyMedian(): f(%f) is %f\n", midpoint, f); - if (fabs(f - getThisValue) <= FLT_EPSILON) { - return (midpoint); - } - - if (f > getThisValue) { - rangeHigh = midpoint; - } else { - rangeLow = midpoint; - } - numIterations++; - } - psTrace(__func__, 4, "---- %s(%f) end ----\n", __func__, midpoint); - return (midpoint); -} - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -#define PS_SQRT(x) sqrt(x) static psF32 QuadraticInverse( psF32 a, psF32 b, psF32 c, - psF32 y) + psF32 y, + psF32 xLo, + psF32 xHi) { psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a)); @@ -1254,13 +1161,20 @@ psF64 x2 = -b/(2.0*a) - tmp; - psF64 y1 = (a * x1 * x1) + (b * x1) + c; - psF64 y2 = (a * x2 * x2) + (b * x2) + c; - - printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c); - printf("QuadraticInverse: y is %f\n", y); - printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2); - printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2); - - return(x1); + if (0) { + psF64 y1 = (a * x1 * x1) + (b * x1) + c; + psF64 y2 = (a * x2 * x2) + (b * x2) + c; + printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c); + printf("QuadraticInverse: y is %f\n", y); + printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2); + printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2); + } + + if ((xLo <= x1) && (x1 <= xHi)) { + return(x1); + } else if ((xLo <= x2) && (x2 <= xHi)) { + return(x2); + } else { + return(0.5 * (xLo + xHi)); + } } @@ -1389,14 +1303,7 @@ (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2])); - // Call p_ps1DPolyMedian(), which does a binary search on the - // polynomial, looking for the value x such that f(x) = yVal psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); - printf("Calling p_ps1DPolyMedian() for the y-value (%.2f) that has an x in the range (%.2f - %.2f)\n", yVal, x->data.F64[0], x->data.F64[2]); - tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal); - printf("Cool, tmpFloat is %.2f\n", tmpFloat); - - QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], 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)) { @@ -1470,47 +1377,7 @@ deriv->data.F32[1] = tmp / (stdev * stdev * stdev); - // printf("f(x, mean, stdev) = f(%.2f %.2f %.2f) is %.2f\n", x, mean, stdev, psGaussian(x, mean, stdev, false)); - psTrace(__func__, 4, "---- %s() end ----\n", __func__); return(psGaussian(x, mean, stdev, false)); } - -//XXX: Use the psLib function? -psF32 LinInterpolate(psF32 x0, - psF32 x1, - psF32 y0, - psF32 y1, - psF32 y) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psTrace(__func__, 5, "LinInterpolate(%.2f %.2f %.2f %.2f %.2f)\n", x0, x1, y0, y1, y); - if (y0 == y1) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: y0 == y1. Cannot interpolate. Returning NaN.\n"); - return(NAN); - } - - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(x0 + y * (x1 - x0) / (y1 - y0)); -} - - - - - - - - - - - - - - - - - - - - /****************************************************************************** @@ -1532,5 +1399,5 @@ XXX: Review and ensure that all memory is free'ed at premature exits. *****************************************************************************/ -#define INITIAL_NUM_BINS 500.0 +#define INITIAL_NUM_BINS 1000.0 psS32 p_psVectorRobustStats(const psVector* myVector, const psVector* errors, @@ -1543,6 +1410,4 @@ PS_VECTOR_PRINT_F32(myVector); } - // psS32 currentId = psMemGetId(); - // psS32 memLeaks = 0; psHistogram *robustHistogram = NULL; psHistogram *cumulativeRobustHistogram = NULL; @@ -1589,5 +1454,4 @@ 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); } psTrace(__func__, 6, "Initial robust bin size is %.2f\n", binSize); @@ -1680,5 +1544,4 @@ // XXX: Check return codes. // - printf("ADD: step 3: Interpolate to the exact 50-percent position: this is the robust histogram median.\n"); stats->robustMedian = fitQuadraticSearchForYThenReturnX( *(psVector* *)&cumulativeRobustHistogram->bounds, @@ -1714,5 +1577,6 @@ } } - psTrace(__func__, 6, "The 15.8655-percent and 84.1345-percent 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)); @@ -1750,24 +1614,51 @@ binHi, totalDataPoints * 0.841345f); - psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32); - } - + psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", + binLoF32, binHiF32); + } + + // + // This code basically interpolates to find the positions exactly. + // if (1) { - psTrace(__func__, 6, "binLo is %d. Nums at that bin and the next are (%.2f, %.2f)\n", binLo, cumulativeRobustHistogram->nums->data.F32[binLo], cumulativeRobustHistogram->nums->data.F32[binLo+1]); - psTrace(__func__, 6, "binHi is %d. Nums at that bin and the next are (%.2f, %.2f)\n", binHi, cumulativeRobustHistogram->nums->data.F32[binHi], cumulativeRobustHistogram->nums->data.F32[binHi+1]); - // XXX: Add checks to ensure that these bins are not the first (seg fault). - binLoF32 = LinInterpolate( - cumulativeRobustHistogram->bounds->data.F32[binLo], - cumulativeRobustHistogram->bounds->data.F32[binLo+1], - cumulativeRobustHistogram->nums->data.F32[binLo-1], - cumulativeRobustHistogram->nums->data.F32[binLo], - totalDataPoints * 0.158655f); - binHiF32 = LinInterpolate( - cumulativeRobustHistogram->bounds->data.F32[binHi], - cumulativeRobustHistogram->bounds->data.F32[binHi+1], - cumulativeRobustHistogram->nums->data.F32[binHi-1], - cumulativeRobustHistogram->nums->data.F32[binHi], - totalDataPoints * 0.841345f); - // XXX: Check for NANs + psTrace(__func__, 6, "binLo is %d. Nums at that bin and the next are (%.2f, %.2f)\n", + binLo, cumulativeRobustHistogram->nums->data.F32[binLo], + cumulativeRobustHistogram->nums->data.F32[binLo+1]); + psTrace(__func__, 6, "binHi is %d. Nums at that bin and the next are (%.2f, %.2f)\n", + binHi, cumulativeRobustHistogram->nums->data.F32[binHi], + cumulativeRobustHistogram->nums->data.F32[binHi+1]); + + psF32 deltaNums; + psF32 deltaBounds; + psF32 prevPixels; + psF32 percentNums; + psF32 base; + deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binLo+1] - cumulativeRobustHistogram->bounds->data.F32[binLo]; + if (binLo == 0) { + deltaNums = cumulativeRobustHistogram->nums->data.F32[0]; + prevPixels = 0; + } else { + deltaNums = cumulativeRobustHistogram->nums->data.F32[binLo] - cumulativeRobustHistogram->nums->data.F32[binLo-1]; + prevPixels = cumulativeRobustHistogram->nums->data.F32[binLo-1]; + } + percentNums = (totalDataPoints * 0.158655f) - prevPixels; + base = cumulativeRobustHistogram->bounds->data.F32[binLo]; + binLoF32 = base + (deltaBounds / deltaNums) * percentNums; + psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", + base, deltaBounds, deltaNums, prevPixels, percentNums); + + deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binHi+1] - cumulativeRobustHistogram->bounds->data.F32[binHi]; + if (binHi == 0) { + deltaNums = cumulativeRobustHistogram->nums->data.F32[0]; + prevPixels = 0; + } else { + deltaNums = cumulativeRobustHistogram->nums->data.F32[binHi] - cumulativeRobustHistogram->nums->data.F32[binHi-1]; + prevPixels = cumulativeRobustHistogram->nums->data.F32[binHi-1]; + } + percentNums = (totalDataPoints * 0.841345f) - prevPixels; + base = cumulativeRobustHistogram->bounds->data.F32[binHi]; + binHiF32 = base + (deltaBounds / deltaNums) * percentNums; + psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", + base, deltaBounds, deltaNums, prevPixels, percentNums); psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32); } @@ -1789,13 +1680,16 @@ if (sigma < (2 * binSize)) { psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize); - psF32 medianLo = robustHistogram->bounds->data.F32[PS_MAX(0, (binMedian - 25))]; - psF32 medianHi = robustHistogram->bounds->data.F32[PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25))]; + psS32 maskLo = PS_MAX(0, (binMedian - 25)); + psS32 maskHi = PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25)); + psF32 medianLo = robustHistogram->bounds->data.F32[maskLo]; + psF32 medianHi = robustHistogram->bounds->data.F32[maskHi]; 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, "The median is at bin number %d. We mask bins outside the bin range (%d:%d)\n", + binMedian, maskLo, maskHi); 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)) { tmpMaskVec->data.U8[i] = 1; - psTrace(__func__, 8, "Masking element %d is %f\n", i, myVector->data.F32[i]); + psTrace(__func__, 6, "Masking element %d is %f\n", i, myVector->data.F32[i]); } } @@ -1847,5 +1741,4 @@ // XXX: Check for errors. // - printf("ADD: step 8: Interpolate for the lower quartile positions.\n"); psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX( *(psVector* *)&cumulativeRobustHistogram->bounds, @@ -1853,5 +1746,4 @@ binLo25, totalDataPoints * 0.25f); - printf("ADD: step 8: Interpolate for the upper quartile positions.\n"); psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX( *(psVector* *)&cumulativeRobustHistogram->bounds, @@ -1891,4 +1783,7 @@ psF32 newBinSize = sigma / dN; + // + // Determine the min/max of the vector (which prior outliers masked out) + // rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); @@ -1903,4 +1798,7 @@ } + // + // Calculate the number of bins. + // numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize); psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max); @@ -1986,5 +1884,4 @@ } psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); - psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); psArray *x = psArrayAlloc((1 + (binMax - binMin))); psS32 j = 0; @@ -1994,5 +1891,4 @@ x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32); ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i); - xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i); j++; } @@ -2005,5 +1901,5 @@ psF32 minY = FLT_MAX; psF32 maxY = -FLT_MAX; - for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { + for (psS32 i = 0 ; i < y->n ; i++) { if (y->data.F32[i] > maxY) { maxY = y->data.F32[i]; @@ -2017,5 +1913,5 @@ // XXX: Use the normalize routines for this. // - for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { + for (psS32 i = 0 ; i < y->n ; i++) { y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY); } @@ -2031,7 +1927,9 @@ PS_VECTOR_PRINT_F32(y); } - psFree(xTmp); + + // + // Fit a Gaussian to the data. + // rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D); - if (rcBool != true) { psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); @@ -2566,4 +2464,6 @@ Returns The stats structure. + +XXX: Should we free stats if the asserts fail? *****************************************************************************/ psStats* psVectorStats(psStats* stats, @@ -2575,5 +2475,5 @@ psTrace(__func__, 3,"---- %s() begin ----\n", __func__); PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(in, stats); + PS_ASSERT_VECTOR_NON_NULL(in, NULL); if (mask != NULL) { PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, stats);