Index: trunk/psLib/src/dataManip/psStats.c =================================================================== --- trunk/psLib/src/dataManip/psStats.c (revision 2074) +++ trunk/psLib/src/dataManip/psStats.c (revision 2197) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.67 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-13 02:40:13 $ + * @version $Revision: 1.68 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:43 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -40,5 +40,5 @@ /*****************************************************************************/ // will use robust statistical methods. -#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. +#define GAUSS_WIDTH 5 // The width of the Gaussian or boxcar smoothing. #define CLIPPED_NUM_ITER_LB 1 #define CLIPPED_NUM_ITER_UB 10 @@ -57,4 +57,5 @@ /* TYPE DEFINITIONS */ /*****************************************************************************/ +psVector* p_psConvertToF32(psVector* in); /*****************************************************************************/ @@ -423,22 +424,4 @@ float rangeMax = 0.0; // Exclude date above this - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - // XXX: This code is no longer used. We now calculate the median exactly - // regardless of the vector size. - /* - * // Calculate the robust quartiles. stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - * p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - * - * // Store the robust quartiles into the sample quartile members. stats->sampleMedian = - * stats2->robustMedian; - * - * // Free temporary data buffers. psFree(stats2); - * - * // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. stats->options = stats->options | - * PS_STAT_ROBUST_FOR_SAMPLE; - * - * return; } */ - // Determine how many data points fit inside this min/max range // and are not masked, IF the maskVector is not NULL> @@ -507,67 +490,66 @@ /****************************************************************************** - This routine smoothes the data in the input robustHistogram with a - Gaussian of width sigma. - - XXX: Must consult with IfA on the proper values for the Gaussian - smoothing. Currently, the Gaussian coefficients are such that only - the center coefficient is non-zero. This is because "e" is being - raised to a very large negative number for all points except the - center point. - - XXX: use a static variable for gaussianCoefs[] and compute them once. - *****************************************************************************/ -psVector* p_psVectorsmoothHistGaussian(psHistogram* robustHistogram, +p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input +robustHistogram with a Gaussian of width sigma. + *****************************************************************************/ +#define GAUSS_WIDTH_OTHER 5.0 +psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, float sigma) { int i = 0; // Loop index variable int j = 0; // Loop index variable - float denom = 0.0; // Temporary variable - float expo = 0.0; // Temporary variable - float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients - psVector* smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32); - - for (i = 0; i < (1 + (2 * GAUSS_WIDTH)); i++) { - if (fabs(sigma) >= FLT_EPSILON) { - // If sigma does not equal zero, then we use Gaussian smoothing. - #ifdef DARWIN - denom = (float)sqrt(2.0 * M_PI * sigma * sigma); - #else - - denom = sqrtf(2.0 * M_PI * sigma * sigma); - #endif - - expo = -(float)((i - GAUSS_WIDTH) * (i - GAUSS_WIDTH)); - expo /= (2.0 * sigma * sigma); - gaussianCoefs[i] = exp(expo / denom); - - //printf("p_psVectorsmoothHistGaussian(): gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); - // NOTE: Gaussian smoothing just isn't working with low sigma - // values. The problem is that the Gaussian coefficients are - // all zero, except for the middle coefficient, which is exactly - // one. Therefore, I'm using boxcar smoothing. - gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float)GAUSS_WIDTH)); - // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); - } else { - /* If sigma equals zero (all pixels have the same value) the above code will divide by zero. - * Therefore, we don't need to smooth the data. */ - for (i = 0; i < robustHistogram->nums->n; i++) { - smooth->data.F32[i] = (float)robustHistogram->nums->data.S32[i]; - } - return (smooth); - } - } - - // Perform the actual smoothing. - for (i = 0; i < robustHistogram->nums->n; i++) { + float iMid; + float jMid; + int numBins = robustHistogram->nums->n; + int numBounds = robustHistogram->bounds->n; + psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32); + int jMin = 0; + int jMax = 0; + float firstBound = robustHistogram->bounds->data.F32[0]; + float lastBound = robustHistogram->bounds->data.F32[numBounds-1]; + psScalar x; + + x.type.type = PS_TYPE_F32; + for (i = 0; i < numBins; i++) { + iMid = (robustHistogram->bounds->data.F32[i] + + robustHistogram->bounds->data.F32[i+1]) / 2.0; + + // 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 - (GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMin = p_psVectorBinDisect(robustHistogram->bounds, &x); + } else if (x.data.F32 <= firstBound) { + jMin = 0; + } else if (x.data.F32 >= lastBound) { + jMin = robustHistogram->bounds->n - 1; + + } + + x.data.F32 = iMid + (GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMax = p_psVectorBinDisect(robustHistogram->bounds, &x); + } 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 = -GAUSS_WIDTH; j <= +GAUSS_WIDTH; j++) { - if (((j + i) >= 0) && ((j + i) < smooth->n)) { - smooth->data.F32[i] += (gaussianCoefs[j + GAUSS_WIDTH] * - (float)robustHistogram->nums->data.S32[j + i]); - } - } - } - return (smooth); + for (j = jMin ; j <= jMax ; j++) { + jMid = (robustHistogram->bounds->data.F32[j] + + robustHistogram->bounds->data.F32[j+1]) / 2.0; + smooth->data.F32[i] += + ((float) robustHistogram->nums->data.U32[j]) * + psGaussian(jMid, iMid, sigma, true); + } + } + + return(smooth); } @@ -913,5 +895,6 @@ specified data point. *****************************************************************************/ -float p_psGaussian(const psVector* restrict myData, const psVector* restrict myParams) +float p_psGaussian(const psVector* restrict myData, + const psVector* restrict myParams) { float x = myData->data.F32[0]; @@ -930,5 +913,6 @@ calculates the specified partial derivative of the above Gaussian function. *****************************************************************************/ -float p_psGaussianDeriv(const psVector* restrict myData, const psVector* restrict myParams, int whichParam) +float p_psGaussianDeriv(const psVector* restrict myData, + const psVector* restrict myParams, int whichParam) { float x = myData->data.F32[0]; @@ -1036,45 +1020,75 @@ /****************************************************************************** -p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given -the specified histogram, with the specified pre-calculated cumulativeSums, -this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum", -and then finds the point x such that f(x) == fitFloat. +fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): This routine +takes psVectors of x/y pairs as input, and fits a quadratic to the 3 points +surrounding element binNum in the vectors (the midpoint between element i +and i+1 is used for x[i]). It then determines for what value x does that +quadratic f(x) = yVal (the input parameter). -XXX: This function is currently not being used. - *****************************************************************************/ -float p_psFitQuadratic(psHistogram* histogram, - psVector* cumulativeSums, - int binNum, - float fitFloat) -{ +// XXX: Use static variables. +*****************************************************************************/ +float fitQuadraticSearchForYThenReturnX(psVector *xVec, + psVector *yVec, + int binNum, + float yVal) +{ + // static psVector* x = NULL; + // static psVector* y = NULL; + // static psVector* yErr = NULL; psVector* x = psVectorAlloc(3, PS_TYPE_F64); psVector* y = psVectorAlloc(3, PS_TYPE_F64); psVector* yErr = psVectorAlloc(3, PS_TYPE_F64); + + /* + if (x == NULL) { + x = psVectorAlloc(3, PS_TYPE_F64); + p_psMemSetPersistent(x, true); + } + if (y == NULL) { + y = psVectorAlloc(3, PS_TYPE_F64); + p_psMemSetPersistent(y, true); + } + if (yErr == NULL) { + yErr = psVectorAlloc(3, PS_TYPE_F64); + p_psMemSetPersistent(yErr, true); + } + */ psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - - if ((binNum > 0) && (binNum < (histogram->nums->n + 1))) { - x->data.F64[0] = (double)0.5 * - (histogram->bounds->data.F32[binNum - 1] + histogram->bounds->data.F32[binNum]); - x->data.F64[1] = (double)0.5 * - (histogram->bounds->data.F32[binNum] + histogram->bounds->data.F32[binNum + 1]); - x->data.F64[2] = (double)0.5 * - (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum + 2]); - - y->data.F64[0] = cumulativeSums->data.F32[binNum - 1]; - y->data.F64[1] = cumulativeSums->data.F32[binNum]; - y->data.F64[2] = cumulativeSums->data.F32[binNum + 1]; - - if (!((y->data.F64[0] <= fitFloat) && (fitFloat <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); - } - + float tmpFloat; + + if ((binNum > 0) && (binNum < (yVec->n - 2))) { + x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); + x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); + x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); + 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]; + + // Ensure that yVal is within the range of the bins we are using. + if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) { + printf("((%f), %f, %f)\n", yVal, y->data.F64[0], y->data.F64[2]); + psError(__func__, "yVal not within y-range\n"); + } yErr->data.F64[0] = 1.0; yErr->data.F64[1] = 1.0; yErr->data.F64[2] = 1.0; + // Determine the coefficients of the polynomial. myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - return (p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat)); + // Call p_ps1DPolyMedian(), which does a binary search on the + // polynomial, looking for the value x such that f(x) = yVal + tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal); } else { - return (0.5 * (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum])); + if (binNum == 0) { + tmpFloat = 0.5 * (xVec->data.F32[binNum] + + xVec->data.F32[binNum + 1]); + } else if (binNum == (xVec->n - 1)) { + // 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]); + } } @@ -1083,6 +1097,9 @@ psFree(yErr); psFree(myPoly); - return (0.0); -} + return(tmpFloat); +} + +#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \ +(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1])) /****************************************************************************** @@ -1111,5 +1128,7 @@ *****************************************************************************/ void p_psVectorRobustStats(const psVector* restrict myVector, - const psVector* restrict maskVector, unsigned int maskVal, psStats* stats) + const psVector* restrict maskVector, + unsigned int maskVal, + psStats* stats) { psHistogram* robustHistogram = NULL; @@ -1118,19 +1137,10 @@ int LQBinNum = -1; // Bin num for lower quartile int UQBinNum = -1; // Bin num for upper quartile + int medianBinNum = -1; int i = 0; // Loop index variable - int maxBinNum = 0; - float maxBinCount = 0.0; + int modeBinNum = 0; + float modeBinCount = 0.0; float dL = 0.0; int numBins = 0; - psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN); - - // psImage* domain; - // psVector* errors; - // psVector* data; - // psVector* initialGuess; - // psVector* theParams; - // float chiSq=0.0; - // float max = -HUGE; - // float max_pos; float myMean = 0.0; float myStdev = 0.0; @@ -1139,14 +1149,11 @@ float sumSquares = 0.0; float sumDiffs = 0.0; - psVector* x = psVectorAlloc(3, PS_TYPE_F64); - psVector* y = psVectorAlloc(3, PS_TYPE_F64); - psVector* yErr = psVectorAlloc(3, PS_TYPE_F64); - psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - psVector* cumulativeRobustSumsFullRange = NULL; - psVector* cumulativeRobustSumsDlRange = NULL; + psVector* cumulativeRobustSums = NULL; float sumRobust = 0.0; float sumN50 = 0.0; float sumNfit = 0.0; - float cumulativeMedian = 0.0; + psScalar tmpScalar; + tmpScalar.type.type = PS_TYPE_F32; + psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN); // Compute the initial bin size of the robust histogram. This is done @@ -1188,4 +1195,5 @@ p_psVectorMax(myVector, maskVector, maskVal, stats); } + // Create the histogram structure. NOTE: we can not specify the bin size // precisely since the argument to psHistogramAlloc() is the number of @@ -1198,7 +1206,6 @@ psVectorHistogram(robustHistogram, myVector, maskVector, maskVal); - // Smooth the histogram. - // XXX: is that the right stdev? - robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram, + // Smooth the histogram, Gaussian-style. + robustHistogramVector = p_psVectorSmoothHistGaussian(robustHistogram, tmpStats->clippedStdev / 4.0f); @@ -1208,87 +1215,62 @@ /************************************************************************** - Determine the lower/upper quartiles. + Determine the median/lower/upper quartile bin numbers. + + We define a vector called "cumulativeRobustSums" where the value at + index position i is equal to the sum of bins 0:i. This will be used in + determining the median and lower/upper quartiles. **************************************************************************/ - // We define a vector called "cumulativeRobustSums..." where the value at - // index position i is equal to the sum of bins 0:i. This will be used - // now and later in determining the lower/upper quartiles. - cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0]; + cumulativeRobustSums = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + cumulativeRobustSums->data.F32[0] = robustHistogramVector->data.F32[0]; for (i = 1; i < robustHistogramVector->n; i++) { - cumulativeRobustSumsFullRange->data.F32[i] = - cumulativeRobustSumsFullRange->data.F32[i - 1] + robustHistogramVector->data.F32[i]; - } - sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1]; - - // Determine the bin number containing the lower quartile point. - LQBinNum = -1; - for (i = 0; i < cumulativeRobustSumsFullRange->n; i++) { - if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust / 4.0)) { - LQBinNum = i; - break; - } - } - - // Determine the bin number containing the upper quartile point. - UQBinNum = -1; - for (i = cumulativeRobustSumsFullRange->n - 1; i >= 0; i--) { - if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0 * sumRobust / 4.0)) { - UQBinNum = i; - break; - } - } - - if ((LQBinNum == -1) || (UQBinNum == -1)) { - psAbort(__func__, "Could not determine the robust lower/upper quartiles."); - } - + cumulativeRobustSums->data.F32[i] = cumulativeRobustSums->data.F32[i - 1] + + robustHistogramVector->data.F32[i]; + } + sumRobust = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; + + tmpScalar.data.F32 = sumRobust / 4.0; + LQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); + tmpScalar.data.F32 = 3.0 * sumRobust / 4.0; + UQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); + tmpScalar.data.F32 = sumRobust / 2.0; + medianBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); + + if ((LQBinNum < 0) || (UQBinNum < 0)) { + psError(__func__, "Could not determine the robust lower/upper quartile bin numbers."); + } + if (medianBinNum < 0) { + psError(__func__, "Could not determine the robust lower/upper quartile bin numbers."); + } /************************************************************************** Determine the mode in the range LQ:UQ. **************************************************************************/ // Determine the bin with the peak value in the range LQ to UQ. - maxBinNum = LQBinNum; - maxBinCount = robustHistogramVector->data.F32[LQBinNum]; - sumN50 = (float)robustHistogram->nums->data.S32[LQBinNum]; + modeBinNum = LQBinNum; + modeBinCount = robustHistogramVector->data.F32[LQBinNum]; + sumN50 = (float)robustHistogram->nums->data.U32[LQBinNum]; for (i = LQBinNum + 1; i <= UQBinNum; i++) { - if (robustHistogramVector->data.F32[i] > maxBinCount) { - maxBinNum = i; - maxBinCount = robustHistogramVector->data.F32[i]; - } - sumN50 += (float)robustHistogram->nums->data.S32[i]; - } - - // XXX: is dL defined as the value at the LQ/UQ, or the bin number? + if (robustHistogramVector->data.F32[i] > modeBinCount) { + modeBinNum = i; + modeBinCount = robustHistogramVector->data.F32[i]; + } + sumN50 += (float)robustHistogram->nums->data.U32[i]; + } + dL = (UQBinNum - LQBinNum) / 4; - - printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum); - /************************************************************************** Determine the mean/stdev for the bins in the range mode-dL to mode+dL **************************************************************************/ - cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - for (i = 0; i < robustHistogramVector->n; i++) { - cumulativeRobustSumsDlRange->data.F32[i] = 0.0; - } - sumNfit = 0.0; - cumulativeMedian = 0.0; - for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) { - if ((0 <= i) && (i < robustHistogramVector->n)) { - cumulativeRobustSumsDlRange->data.F32[i] = - cumulativeRobustSumsDlRange->data.F32[i - 1] + robustHistogramVector->data.F32[i]; - cumulativeMedian += robustHistogramVector->data.F32[i]; - sumNfit += (float)robustHistogram->nums->data.S32[i]; - } - } - // Calculate the mean of the smoothed robust histogram in the range // mode-dL to mode+dL. We use the midpoint of each bin as the mean for // that bin (this is a non-exact approximation). + sumNfit = 0.0; myMean = 0.0; - for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) { + for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { if ((0 <= i) && (i < robustHistogramVector->n)) { - myMean += (robustHistogramVector->data.F32[i]) * 0.5 * - (robustHistogram->bounds->data.F32[i + 1] + robustHistogram->bounds->data.F32[i]); + myMean += (robustHistogramVector->data.F32[i]) * PS_BIN_MIDPOINT(robustHistogram, i); countFloat += robustHistogramVector->data.F32[i]; } + + sumNfit += (float)robustHistogram->nums->data.U32[i]; } myMean /= countFloat; @@ -1297,8 +1279,7 @@ // mode-dL to mode+dL. We use the midpoint of each bin as the mean for // that bin. - for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) { + for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { if ((0 <= i) && (i < robustHistogramVector->n)) { - diff = (0.5 * (robustHistogram->bounds->data.F32[i + 1] + - robustHistogram->bounds->data.F32[i])) - myMean; + diff = PS_BIN_MIDPOINT(robustHistogram, i) - myMean; sumSquares += diff * diff * robustHistogramVector->data.F32[i]; sumDiffs += diff * robustHistogramVector->data.F32[i]; @@ -1315,6 +1296,5 @@ if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = 0.5 * - (robustHistogram->bounds->data.F32[maxBinNum] + robustHistogram->bounds->data.F32[maxBinNum + 1]); + stats->robustMode = PS_BIN_MIDPOINT(robustHistogram, modeBinNum); } @@ -1322,140 +1302,42 @@ stats->robustStdev = myStdev; } - // To determine the median (and later, the lower/upper quartile), we fit - // a quadratic to the three bins surrounding the bin containing the median. - // The quadratic y=f(x) with x being the midpoint of each bin, and y being - // the cumulative number of data points in all bins up to, and including, - // this bin. We then solve the quadratic for - - // XXX: Create a "p_psQuadraticFitAndSolveForX() function. - + + // To determine the median, we fit a quadratic y=f(x) to the three bins + // surrounding the bin containing the median (x is the midpoint of each + // bin and y is the value of each bin). Then we figure out what value + // of x corresponds to f(x) being the median (half of all points). if (stats->options & PS_STAT_ROBUST_MEDIAN) { - if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n - 1))) { - x->data.F64[0] = (double)0.5 * - (robustHistogram->bounds->data.F32[maxBinNum - 1] + - robustHistogram->bounds->data.F32[maxBinNum]); - x->data.F64[1] = (double)0.5 * - (robustHistogram->bounds->data.F32[maxBinNum] + - robustHistogram->bounds->data.F32[maxBinNum + 1]); - x->data.F64[2] = (double)0.5 * - (robustHistogram->bounds->data.F32[maxBinNum + 1] + - robustHistogram->bounds->data.F32[maxBinNum + 2]); - - y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum - 1]; - y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum]; - y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum + 1]; - - // Ensure that cumulativeMedian/2 is actually within the range of the bins - // we are using. - cumulativeMedian *= 0.5; - if (!((y->data.F64[0] <= cumulativeMedian) && (cumulativeMedian <= y->data.F64[2]))) { - printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]); - psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n"); - } - // XXX: yErr is not currently used by psVectorFitPolynomial1D(). We - // may have to set this meaningfully later. - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - // Determine the coefficients of the polynomial. - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - // Call p_ps1DPolyMedian(), which does a binary search on the - // polynomial, looking for the value x such that - // f(x) = cumulativeMedian. - stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], cumulativeMedian); - } else { - // If the mode is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum + 1] + - robustHistogram->bounds->data.F32[maxBinNum]); - } - } - // The lower/upper quartile calculations are very similar to the median - // calculations. We fit a quadratic to the array containing the - // cumulative data points in each bin, then search for x such that - // f(x) equals the lower/upper quartile exactly. - // + // Take a psVector. Fit a polynomial y = f(x) to the 3 data elements + // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0. + stats->robustMedian = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, + cumulativeRobustSums, + medianBinNum, + sumRobust/2.0); + } + + // To determine the quartiles, we fit a quadratic y=f(x) to the three bins + // surrounding the bin containing LQ/UQ (x is the midpoint of each + // bin and y is the value of each bin). Then we figure out what value + // of x corresponds to f(x) being the LQ/UQ. if (stats->options & PS_STAT_ROBUST_QUARTILE) { - countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1]; - - if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n - 1))) { - x->data.F64[0] = (double)0.5 * - (robustHistogram->bounds->data.F32[LQBinNum - 1] + - robustHistogram->bounds->data.F32[LQBinNum]); - x->data.F64[1] = (double)0.5 * - (robustHistogram->bounds->data.F32[LQBinNum] + - robustHistogram->bounds->data.F32[LQBinNum + 1]); - x->data.F64[2] = (double)0.5 * - (robustHistogram->bounds->data.F32[LQBinNum + 1] + - robustHistogram->bounds->data.F32[LQBinNum + 2]); - - y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum - 1]; - y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum]; - y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum + 1]; - - if (!((y->data.F64[0] <= (countFloat / 4.0)) && ((countFloat / 4.0) <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - stats->robustLQ = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], countFloat / 4.0); - - } else { - // If the LQ is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum + 1] + - robustHistogram->bounds->data.F32[LQBinNum]); - } - - if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n - 1))) { - x->data.F64[0] = (double)0.5 * - (robustHistogram->bounds->data.F32[UQBinNum - 1] + - robustHistogram->bounds->data.F32[UQBinNum]); - x->data.F64[1] = (double)0.5 * - (robustHistogram->bounds->data.F32[UQBinNum] + - robustHistogram->bounds->data.F32[UQBinNum + 1]); - x->data.F64[2] = (double)0.5 * - (robustHistogram->bounds->data.F32[UQBinNum + 1] + - robustHistogram->bounds->data.F32[UQBinNum + 2]); - - y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum - 1]; - y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum]; - y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum + 1]; - - if (!((y->data.F64[0] <= (3.0 * countFloat / 4.0)) && - ((3.0 * countFloat / 4.0) <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - stats->robustUQ = p_ps1DPolyMedian(myPoly, - x->data.F64[0], x->data.F64[2], 3.0 * countFloat / 4.0); - } else { - // If the UQ is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum + 1] + - robustHistogram->bounds->data.F32[UQBinNum]); - } - } + countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; + + stats->robustLQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, + cumulativeRobustSums, + LQBinNum, + countFloat/4.0); + stats->robustUQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, + cumulativeRobustSums, + UQBinNum, + 3.0 * countFloat/4.0); + } + // XXX: I think sumNfit == sumN50 here. stats->robustNfit = sumNfit; stats->robustN50 = sumN50; - psFree(x); - psFree(y); - psFree(yErr); psFree(tmpStats); psFree(robustHistogram); - psFree(myPoly); - psFree(cumulativeRobustSumsFullRange); - psFree(cumulativeRobustSumsDlRange); + psFree(robustHistogramVector); + psFree(cumulativeRobustSums); } @@ -1468,6 +1350,4 @@ // I am commenting this code out, and replacing it with code which // calculates the mean directly on the robustHistogram. - // - // XXX: Replace this with the LaGrange interpolation stuff. domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32); @@ -1479,6 +1359,5 @@ max_pos = -1; for (i=0;in;i++) { - domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] + - robustHistogram->bounds->data.F32[i]); + domain->data.F32[i][0] = PS_BIN_MIDPOINT(robustHistogram, i);bounds->data.F32[modeBinNum + 1]); data->data.F32[i] = (float) robustHistogramVector->data.F32[i]; errors->data.F32[i] = 1.0; @@ -1695,4 +1574,6 @@ psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; + psVector* inF32; + int mustFreeTmp = 1; // NOTE: Verify that this is the correct action. @@ -1706,5 +1587,5 @@ if (out->nums->type.type != PS_TYPE_U32) { - psAbort(__func__, "Only data type PS_TYPE_U32 for output.nums member."); + psAbort(__func__, "Only data type PS_TYPE_U32 for out->nums member."); } // NOTE: Verify that this is the correct action. @@ -1713,15 +1594,17 @@ } - if (in->type.type != PS_TYPE_F32) { - psAbort(__func__, "Only data type PS_TYPE_F32 is currently supported (0x%x).", in->type.type); - } - if (mask != NULL) { if (in->n != mask->n) { - psAbort(__func__, "Vector data and vector mask are of different sizes."); + psError(__func__, "Vector data and vector mask are of different sizes."); } if (mask->type.type != PS_TYPE_U8) { - psAbort(__func__, "Vector mask must be type PS_TYPE_U8"); - } + psError(__func__, "Vector mask must be type PS_TYPE_U8"); + } + } + + inF32 = p_psConvertToF32((psVector *) in); + if (inF32 == NULL) { + inF32 = (psVector *) in; + mustFreeTmp = 0; } @@ -1730,14 +1613,14 @@ numBins = out->nums->n; - for (i = 0; i < in->n; i++) { + for (i = 0; i < inF32->n; i++) { // Check if this pixel is masked, and if so, skip it. if ((mask == NULL) || ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) { // Check if this pixel is below the minimum value, and if so // count it, then skip it. - if (in->data.F32[i] < out->bounds->data.F32[0]) { + if (inF32->data.F32[i] < out->bounds->data.F32[0]) { out->minNum++; // Check if this pixel is above the maximum value, and if so // count it, then skip it. - } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { + } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) { out->maxNum++; } else { @@ -1746,5 +1629,5 @@ if (out->uniform == true) { binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; - binNum = (int)((in->data.F32[i] - out->bounds->data.F32[0]) / binSize); + binNum = (int)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); // NOTE: This next if-statement really shouldn't be necessary. @@ -1755,16 +1638,15 @@ } - (out->nums->data.S32[binNum])++; + (out->nums->data.U32[binNum])++; // If this is a non-uniform histogram, determining the correct // bin number requires a bit more work. } else { - tmpScalar.data.F32 = in->data.F32[i]; - tmp = p_psVectorBinDisect(out->bounds, - &tmpScalar); - if (tmp == -1) { - // XXX: Generate warning message + tmpScalar.data.F32 = inF32->data.F32[i]; + tmp = p_psVectorBinDisect(out->bounds, &tmpScalar); + if (tmp < 0) { + // XXX: Generate warning message } else { - (out->nums->data.S32[tmp])++; + (out->nums->data.U32[tmp])++; } } @@ -1772,19 +1654,22 @@ } } + + if (mustFreeTmp == 1) { + psFree(inF32); + } return (out); } /****************************************************************************** -p_psConvertToF32(myVector, maskVector, maskVal, stats): this is the cheap -way to support a variety of vector data types: we simply convert the input -vector to F32 at the beginning, and write all of our functions in F32. If -the vast majority of all vector stat operations are F32 (or any other single -type), then this is probably the best way to go. Otherwise, when the -algorithms stablize, we will then macro everything and put type support in -the various stat functions. +p_psConvertToF32(in): this is the cheap way to support a variety of vector +data types: we simply convert the input vector to F32 at the beginning, and +write all of our functions in F32. If the vast majority of all vector stat +operations are F32 (or any other single type), then this is probably the +best way to go. Otherwise, when the algorithms stablize, we will then macro +everything and put type support in the various stat functions. XXX: Should the default data type be F64? Since we are buying Athlons... *****************************************************************************/ -psVector* p_psConvertToF32(psStats* stats, psVector* in, psVector* mask, unsigned int maskVal) +psVector* p_psConvertToF32(psVector* in) { int i = 0; @@ -1800,4 +1685,9 @@ for (i = 0; i < in->n; i++) { tmp->data.F32[i] = (float)in->data.U16[i]; + } + } else if (in->type.type == PS_TYPE_U8) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (float)in->data.U8[i]; } } else if (in->type.type == PS_TYPE_F64) { @@ -1839,20 +1729,18 @@ int mustFreeTmp = 1; - // NOTE: Verify that this is the correct action. if (in == NULL) { - return (stats); - } - if (stats == NULL) { - return (NULL); - } - - inF32 = p_psConvertToF32(stats, in, mask, maskVal); + psError(__func__, "in is NULL\n"); + return(stats); + } + PS_CHECK_NULL_PTR_RETURN_NULL(stats); + + inF32 = p_psConvertToF32((psVector *) in); if (inF32 == NULL) { inF32 = in; mustFreeTmp = 0; } - // XXX: Should we abort if (stats->min == stats->max)? if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) { - psAbort(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max); + psError(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max); + return(stats); }