Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 4898) +++ trunk/psLib/src/math/psStats.c (revision 4958) @@ -14,6 +14,6 @@ * stats->binsize * - * @version $Revision: 1.139 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-08-30 01:14:13 $ + * @version $Revision: 1.140 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-07 21:35:50 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -46,5 +46,6 @@ /* DEFINE STATEMENTS */ /*****************************************************************************/ -#define PS_GAUSS_WIDTH 5 // The width of the Gaussian or boxcar smoothing. +#define PS_GAUSS_WIDTH 5 // The width of the Gaussian smoothing. +// This corresponds to N in the ADD. #define PS_CLIPPED_NUM_ITER_LB 1 #define PS_CLIPPED_NUM_ITER_UB 10 @@ -292,4 +293,6 @@ max of the input vector. If there was a problem with the max calculation, this routine sets stats->max to NAN. + +XXX: Do we need to factor errors into it? *****************************************************************************/ psS32 p_psVectorMax(const psVector* myVector, @@ -614,5 +617,6 @@ /****************************************************************************** p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input -robustHistogram with a Gaussian of width sigma. +robustHistogram with a Gaussian of width sigma. It returns a psVector of the +smoothed data. XXX: Only PS_TYPE_F32 is supported. @@ -621,77 +625,227 @@ call that. Is that possible? *****************************************************************************/ -psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, +psVector *p_psVectorSmoothHistGaussian(psHistogram *histogram, psF32 sigma) { - PS_ASSERT_PTR_NON_NULL(robustHistogram, NULL); - PS_ASSERT_PTR_NON_NULL(robustHistogram->bounds, NULL); - - psS32 i = 0; // Loop index variable - psS32 j = 0; // Loop index variable - psF32 iMid; - psF32 jMid; - psS32 numBins = robustHistogram->nums->n; - psS32 numBounds = robustHistogram->bounds->n; - psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32); + PS_ASSERT_PTR_NON_NULL(histogram, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); + + psS32 numBins = histogram->nums->n; + psS32 numBounds = histogram->bounds->n; + psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); + psF32 firstBound = histogram->bounds->data.F32[0]; + psF32 lastBound = histogram->bounds->data.F32[numBounds-1]; + psScalar x; + x.type.type = PS_TYPE_F32; psS32 jMin = 0; psS32 jMax = 0; - psF32 firstBound = robustHistogram->bounds->data.F32[0]; - psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1]; + + if (histogram->uniform == false) { + // + // We get here if the histogram is non-uniform. + // + + for (psS32 i = 0; i < numBins; i++) { + // Determine the midpoint of bin i. + psS32 iMid = PS_BIN_MIDPOINT(histogram, i); + + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + // + x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); + if (jMin < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } + } else if (x.data.F32 <= firstBound) { + jMin = 0; + } else if (x.data.F32 >= lastBound) { + jMin = histogram->bounds->n - 1; + } + + x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); + if (jMax < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } + } else if (x.data.F32 <= firstBound) { + jMax = 0; + } else if (x.data.F32 >= lastBound) { + jMax = histogram->bounds->n - 1; + } + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + for (psS32 j = jMin ; j <= jMax ; j++) { + psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } + } + } else { + // + // We get here if the histogram is uniform. + // + + for (psS32 i = 0; i < numBins; i++) { + psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0]; + psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize); + + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + // + psS32 jMin = i - gaussWidth; + if (jMin < 0 ) { + jMin = 0; + } + psS32 jMax = i + gaussWidth; + if (jMax > (histogram->bounds->n - 1)) { + jMax = (histogram->bounds->n - 1); + } + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + psS32 iMid = PS_BIN_MIDPOINT(histogram, i); + for (psS32 j = jMin ; j <= jMax ; j++) { + psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } + } + } + + return(smooth); +} +/****************************************************************************** +p_psVectorSmoothHistGaussianNEW(): This routine smoothes the data in the input +robustHistogram with a Gaussian of width sigma. It returns a psVector of the +smoothed data. + +XXX: Only PS_TYPE_F32 is supported. + +XXX: Write a general routine which smoothes a psVector. This routine should +call that. Is that possible? + *****************************************************************************/ +psVector *p_psVectorSmoothHistGaussianNEW(psHistogram *histogram, + psF32 sigma) +{ + PS_ASSERT_PTR_NON_NULL(histogram, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); + + psS32 numBins = histogram->nums->n; + psS32 numBounds = histogram->bounds->n; + psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); + psF32 firstBound = histogram->bounds->data.F32[0]; + psF32 lastBound = histogram->bounds->data.F32[numBounds-1]; psScalar x; - x.type.type = PS_TYPE_F32; - for (i = 0; i < numBins; i++) { - // Determine the midpoint of bin i. - iMid = (robustHistogram->bounds->data.F32[i] + - robustHistogram->bounds->data.F32[i+1]) / 2.0; - - - // We determine the bin numbers corresponding to a range of data - // values surrounding iMid. The ranges is of size - // s*PS_GAUSS_WIDTH*sigma - - // YYY: The p_psVectorBinDisect() routine does much of the work of - // the following conditionals, however, it also reports a warning - // message. I don't want the warning message so I reproduce the - // conditionals here. Maybe p_psVectorBinDisect() should not produce - // warnings? - - x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); - if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { - jMin = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x); - if (jMin < 0) { - psError(PS_ERR_UNEXPECTED_NULL, - false, - PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); - return(NULL); - } - } else if (x.data.F32 <= firstBound) { - jMin = 0; - } else if (x.data.F32 >= lastBound) { - jMin = robustHistogram->bounds->n - 1; - } - - x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); - if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { - jMax = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x); - if (jMax < 0) { - psError(PS_ERR_UNEXPECTED_NULL, - false, - PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); - return(NULL); - } - } else if (x.data.F32 <= firstBound) { - jMax = 0; - } else if (x.data.F32 >= lastBound) { - jMax = robustHistogram->bounds->n - 1; - } - - smooth->data.F32[i] = 0.0; - for (j = jMin ; j <= jMax ; j++) { - jMid = (robustHistogram->bounds->data.F32[j] + - robustHistogram->bounds->data.F32[j+1]) / 2.0; - smooth->data.F32[i] += - robustHistogram->nums->data.F32[j] * - psGaussian(jMid, iMid, sigma, true); + psS32 jMin = 0; + psS32 jMax = 0; + + if (histogram->uniform == false) { + // + // We get here if the histogram is non-uniform. + // + + for (psS32 i = 0; i < numBins; i++) { + // Determine the midpoint of bin i. + psS32 iMid = PS_BIN_MIDPOINT(histogram, i); + + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + // + x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); + if (jMin < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } + } else if (x.data.F32 <= firstBound) { + jMin = 0; + } else if (x.data.F32 >= lastBound) { + jMin = histogram->bounds->n - 1; + } + + x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); + if (jMax < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } + } else if (x.data.F32 <= firstBound) { + jMax = 0; + } else if (x.data.F32 >= lastBound) { + jMax = histogram->bounds->n - 1; + } + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + for (psS32 j = jMin ; j <= jMax ; j++) { + psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } + } + } else { + // + // We get here if the histogram is uniform. + // + + for (psS32 i = 0; i < numBins; i++) { + psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0]; + psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize); + + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + // + psS32 jMin = i - gaussWidth; + if (jMin < 0 ) { + jMin = 0; + } + psS32 jMax = i + gaussWidth; + if (jMax > (histogram->bounds->n - 1)) { + jMax = (histogram->bounds->n - 1); + } + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + psS32 iMid = PS_BIN_MIDPOINT(histogram, i); + for (psS32 j = jMin ; j <= jMax ; j++) { + psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } } } @@ -1142,5 +1296,5 @@ These macros and functions define the following functions: -< p_psNormalizeVectorRange(myData, low, high) + p_psNormalizeVectorRange(myData, low, high) That assumes that the low/high arguments are PS_TYPE_F64; the vector myData @@ -1251,11 +1405,10 @@ /****************************************************************************** -p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes -as input a 1-D polynomial of arbitrary order (though we are using 2nd-order -polynomials here) 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) == midpoint. This functions uses a binary-search algorithm on the -range and assumes that the polynomial is monotonically increasing or -decreasing within that range. +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. @@ -1269,16 +1422,10 @@ { PS_ASSERT_POLY_NON_NULL(myPoly, NAN); - PS_FLOAT_COMPARE(rangeLow, rangeHigh, 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 - ); + psF32 fLo = psPolynomial1DEval(myPoly, rangeLow); + psF32 fHi = psPolynomial1DEval(myPoly, rangeHigh); if (!((fLo <= getThisValue) && (fHi >= getThisValue))) { psError(PS_ERR_UNKNOWN, @@ -1300,8 +1447,5 @@ oldMidpoint = midpoint; - f = psPolynomial1DEval( - myPoly, - midpoint - ); + f = psPolynomial1DEval(myPoly, midpoint); if (fabs(f - getThisValue) <= FLT_EPSILON) { return (midpoint); @@ -1331,4 +1475,8 @@ XXX: the vectors do not have to be the same length. Must insert the proper tests to ensure that binNum is within acceptable ranges for both vectors. + +XXX: This currently assumes that the three points are monotonically increasing +or decreasing: so, it works for the cumulative histogram vectors, but not for +arbitrary vectors. We should probably test that condition. *****************************************************************************/ psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec, @@ -1345,8 +1493,4 @@ PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (yVec->n - 1), NAN); - // PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64); - // PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64); - // PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64); - // PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD); psVector *x = psVectorAlloc(3, PS_TYPE_F64); psVector *y = psVectorAlloc(3, PS_TYPE_F64); @@ -1364,4 +1508,30 @@ y->data.F64[1] = yVec->data.F32[binNum]; y->data.F64[2] = yVec->data.F32[binNum + 1]; + + // + // Ensure that the y values are monotonic. + // + // XXX: This routine should probably be rewritten in a more general fashion + // so that the folloiwng checks are not necessary. + // + if (y->data.F64[0] < y->data.F64[1]) { + if (!(y->data.F64[1] <= y->data.F64[2])) { + psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); + return(NAN); + } + } else { + if (!(y->data.F64[1] >= y->data.F64[2])) { + psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); + return(NAN); + } + } // Ensure that yVal is within the range of the bins we are using. @@ -1637,5 +1807,5 @@ // code no longer produces sensible results. // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably - // remove some of theabove code that calculated the initial estimate + // remove some of the above code that calculated the initial estimate // for the mean and sigma. @@ -1645,6 +1815,6 @@ int index = i - modeBinNum + dL; // XXX: Should this be the natural log? - y->data.F32[index] = robustHistogramVector->data.F32[i]; - // y->data.F32[index] = logf(robustHistogramVector->data.F32[i]); + // y->data.F32[index] = robustHistogramVector->data.F32[i]; + y->data.F32[index] = logf(robustHistogramVector->data.F32[i]); x->data.F32[index] = (psF32) index; } @@ -1742,6 +1912,363 @@ psFree(robustHistogramVector); psFree(cumulativeRobustSums); + return(0); } + + +/***************************************************************************** +XXX: Is there a psLib function for this? + *****************************************************************************/ +psVector *PsVectorDup(psVector *in) +{ + psVector *out = psVectorAlloc(in->n, in->type.type); + + if (in->type.type == PS_TYPE_F32) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F32[i] = in->data.F32[i]; + } + } else if (in->type.type == PS_TYPE_F64) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F64[i] = in->data.F64[i]; + } + } else { + printf("XXX: Generate an error here.\n"); + return(NULL); + } + return(out); +} + +/****************************************************************************** +XXX: This function need to be written. Actually, it simply needs to be +retrieved from the CVS repository, since it was written earlier, then +discarded. At present, it was deleted from the CVS repository, so we might +have to retrieve it from tape. +*****************************************************************************/ +psVector *Fit1DGaussian(psVector *x, psVector*y) +{ + printf("XXX: Generate an error here.\n"); + printf("XXX: Error: This function was previously part of psStats.c, was removed, was purged from CVS, and now needs to be retrieved from tape.\n"); + return(NULL); +} + +/****************************************************************************** + +p_psVectorRobustStatsNew(myVector, maskVector, maskVal, stats): This is the new +version of the robust stats routine. + +XXX: MUST DO: If the errors in the input values are known, then the same +approach is used, except that the histograms become probability density +functions (PDFs). In this case, the input values are spread out, so that they +do not simply contribute a single unit to the histogram, but rather contribute +a fraction of a value, equivalent to the weight. In the interests of speed, a +boxcar PDF may be used to represent each input value (as opposed to a +Gaussian), where the boxcar width is equal to 2p2 ln 2 times the error and +each input value contributes constant area. Then the robust median and +standard deviation are estimated in the same manner as above. + +XXX: Check for errors in psLib routines that we call. +*****************************************************************************/ +psS32 p_psVectorRobustStatsNew(const psVector* myVector, + const psVector* errors, + const psVector* maskVector, + psU32 maskVal, + psStats* stats) +{ + psHistogram *robustHistogram = NULL; + psHistogram *cumulativeRobustHistogram = NULL; + psS32 numBins = 0; + psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32); + tmpScalar->type.type = PS_TYPE_F32; + psS32 totalDataPoints = 0; + psS32 rc = 0; + psVector *tmpMaskVec = PsVectorDup((psVector *) maskVector); + + while (1) { + // + // Determine the bin size of the robust histogram. This is done + // by computing the total range of data values and dividing by 1000.0. + // + psStats* tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); + rc = p_psVectorMin(myVector, tmpMaskVec, maskVal, tmpStatsMinMax); + rc|= p_psVectorMax(myVector, tmpMaskVec, maskVal, tmpStatsMinMax); + if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(tmpStatsMinMax); + psFree(tmpMaskVec); + psFree(tmpScalar); + return(1); + } + psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f; + + // + // If all data points have the same value, then we set the appropiate + // members of stats and return. + // + if (fabs(tmpStatsMinMax->max - tmpStatsMinMax->min) <= FLT_EPSILON) { + if (stats->options & PS_STAT_ROBUST_MEDIAN) { + stats->robustMedian = tmpStatsMinMax->min; + } + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + stats->robustUQ = tmpStatsMinMax->min; + stats->robustLQ = tmpStatsMinMax->min; + } + // XXX: Set these to the number of unmasked data points? + stats->robustNfit = 0.0; + stats->robustN50 = 0.0; + psFree(tmpStatsMinMax); + psFree(tmpMaskVec); + psFree(tmpScalar); + + return(0); + } + + // + // ADD: Step 0. + // Construct the histogram with the specified bin size. + // + // NOTE: we can not specify the bin size precisely since the argument + // to psHistogramAlloc() is the number of bins, not the binSize. + // If we get here, we know that binSize != 0.0. + // + numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize); + robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); + cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); + + // Populate the histogram array. + psVectorHistogram(robustHistogram, myVector, errors, tmpMaskVec, maskVal); + + // + // ADD: Step 1. + // Construct the cumulative histogram from the specific histogram + // + cumulativeRobustHistogram->nums->data.F32[0] = robustHistogram->nums->data.F32[0]; + for (psS32 i = 1 ; i < robustHistogram->nums->n ; i++) { + cumulativeRobustHistogram->nums->data.F32[i] = cumulativeRobustHistogram->nums->data.F32[i-1] + + robustHistogram->nums->data.F32[i]; + } + + // + // ADD: Step 2. + // Find the bin which contains the 50% data point. + // + totalDataPoints = cumulativeRobustHistogram->nums->data.F32[numBins - 1]; + tmpScalar->data.F32 = totalDataPoints/2.0; + psS32 binMedian = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if (binMedian != 0) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50% data point.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + return(1); + } + + // + // ADD: Step 3. + // Interpolate to the exact 50% position: this is the robust histogram median. + // XXX: Check for errors here! + // + stats->robustMedian = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binMedian, + totalDataPoints/2.0); + + // + // ADD: Step 4. + // Find the bins which contains the 15.8655% and 84.1345% data points. + // + tmpScalar->data.F32 = totalDataPoints * 0.158655f; + psS32 binLo = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + tmpScalar->data.F32 = totalDataPoints * 0.841345f; + psS32 binHi = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if ((binLo != 0) || (binHi != 0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the15.8655% and 84.1345% data point\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + return(1); + } + + // + // ADD: Step 4b. + // Interpolate Sigma to find these two positions exactly: these are the 1sigma positions. + // + psF32 binLoF32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binLo, + totalDataPoints * 0.158655f); + psF32 binHiF32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binHi, + totalDataPoints * 0.841345f); + + // + // ADD: Step 5. + // Determine SIGMA as 1/2 of the distance between these positions. + // + psF32 sigma = (binHiF32 - binLoF32) / 2.0; + + // + // 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)) { + psF32 medianLo = robustHistogram->bounds->data.F32[binMedian - 25]; + psF32 medianHi = robustHistogram->bounds->data.F32[binMedian + 25]; + for (psS32 i = 0 ; i < myVector->n ; i++) { + if ((myVector->data.F32[i] < medianLo) || + (myVector->data.F32[i] > medianHi)) { + tmpMaskVec->data.U8[i] = 1; + } + } + } else { + // + // ADD: Step 7. + // Find the bins which contains the 25% and 75% data points. + // + tmpScalar->data.F32 = totalDataPoints * 0.25f; + psS32 binLo25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + tmpScalar->data.F32 = totalDataPoints * 0.75f; + psS32 binHi25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if ((binLo25 != 0) || (binHi25 != 0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25% and 75% data points\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + return(1); + } + + // + // ADD: Step 8. + // Interpolate to find these two positions exactly: these are the upper + // and lower quartile positions. + // XXX: Check for errors. + // + psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binLo25, + totalDataPoints * 0.25f); + psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binHi25, + totalDataPoints * 0.75f); + + stats->robustLQ = binLo25F32; + stats->robustUQ = binHi25F32; + // XXX: No idea how to calculate stats->stdev + + // PS_BIN_MIDPOINT(robustHistogram, modeBinNum); + + // XXX: I think sumNfit == sumN50 here. + stats->robustNfit = -1; + stats->robustN50 = -1; + + // + // Perform the Robust Histogram Statistics algorithm above + // + + // + // Smooth the resulting histogram with a Gaussian with SIGMA_s = 1 + // bin. + // + // XXX: SIGMA_s is defined nowhere in the ADD. + // + psF32 SIGMA_S = 1.0; + p_psVectorSmoothHistGaussian(robustHistogram, SIGMA_S); + + // + // Find the bin with the peak value in the range 2 SIGMA of the + // robust histogram median. + // + // XXX: SIGMA is defined nowhere in the ADD. + // + psF32 SIGMA = 2.0; + psS32 binMin = binMedian - (SIGMA * PS_GAUSS_WIDTH); + if (binMin < 0) { + binMin = 0; + } + psS32 binMax = binMedian + (2 * PS_GAUSS_WIDTH); + if (binMin > (robustHistogram->nums->n - 1)) { + binMin = (robustHistogram->nums->n - 1); + } + psS32 binNum = binNum; + psF32 binMaxNums = robustHistogram->nums->data.F32[binNum]; + for (psS32 i = binNum+1 ; i <= binMax ; i++) { + if (robustHistogram->nums->data.F32[i] > binMaxNums) { + binNum = i; + binMaxNums = robustHistogram->nums->data.F32[i]; + } + } + + // + // Fit a Gaussian to the bins in the range 2 SIGMA of the robust + // histogram median. + // + // XXX: SIGMA is defined nowhere in the ADD. + // + psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psVector *x = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psS32 j = 0; + for (psS32 i = binNum ; i <= binMax ; i++) { + y->data.F32[j] = robustHistogram->nums->data.F32[i]; + x->data.F32[j] = PS_BIN_MIDPOINT(robustHistogram, i); + } + // + // XXX: This function need to be written. Actually, it simply + // needs to be retrieved from the CVS repository, since it was + // written earlier, then discarded. At present, it was deleted + // from the CVS repository, so we might have to retrieve it from + // tape. + // + psVector *params = Fit1DGaussian(x, y); + + // + // The robust mean mean_r is derived directly from the fitted + // Gaussian mean. + // + stats->robustMean = params->data.F32[0]; + + // + // The robust standard deviation, SIGMA_r is determined by + // subtracting the smoothing scale in quadrature: SIGMA_r^2 = SIGMA^2 + // - SIGMA_s^2 + // + // XXX: SIGMA and SIGMA_s are defined nowhere in the ADD. We must figure + // out what they are. + // + stats->robustStdev = sqrt(PS_SQR(SIGMA) - PS_SQR(SIGMA_S)); + + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(params); + + return(0); + } + + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + } + return(1); +} + + + + + + + + /*****************************************************************************/