Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 724) +++ trunk/psLib/src/math/psStats.c (revision 729) @@ -14,4 +14,15 @@ #include "psSort.h" +#include "float.h" +#include +#define ROBUST_SIZE_THRESHOLD 10000 // Vectors that are large than this +// will use robust statistical methods. +#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar +// smoothing. +#define PI 3.141592653 +#define CLIPPED_NUM_ITER_LB 1 +#define CLIPPED_NUM_ITER_UB 10 +#define CLIPPED_SIGMA_LB 1.0 +#define CLIPPED_SIGMA_UB 10.0 /****************************************************************************** @@ -78,5 +89,5 @@ newHist->minNum = 0; newHist->maxNum = 0; - + newHist->numBins = numBins; return(newHist); } @@ -156,4 +167,20 @@ } +void p_printVector(psVector *myVector, + psVector *maskVector, + unsigned int maskVal, + psStats *newStruct) +{ + int i = 0; + + for (i=0;in;i++) { + if (maskVector != NULL) + printf("Element %d is %f (mask is %d)\n", i, myVector->vec.f[i], maskVector->vec.ui8[i]); + else + printf("Element %d is %f\n", i, myVector->vec.f[i]); + } +} + + /****************************************************************************** MISC STATISTICAL FUNCTIONS @@ -262,5 +289,4 @@ -#define MEDIAN_SIZE_THRESHOLD 10000 void p_psArraySampleMedian(const psVector *restrict myVector, const psVector *restrict maskVector, @@ -278,10 +304,13 @@ } unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues); + unsortedVector->n = unsortedVector->nalloc; sortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues); + sortedVector->n = sortedVector->nalloc; if (maskVector != NULL) { for (i=0;in;i++) { if (!(maskVal & maskVector->vec.ui8[i])) { - unsortedVector->vec.f[count++] = maskVector->vec.f[i]; + unsortedVector->vec.f[count++] = myVector->vec.f[i]; + } } @@ -303,46 +332,242 @@ } -void p_psArrayRobustMedian(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *newStruct) +void p_psArraysmoothHistGaussian(psHistogram *robustHistogram, + float sigma) +{ + int i = 0; + int j = 0; + float tmpf = 0.0; + float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; + // The coefficients used in the histogram + // smoothing calculation. + + 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 + tmpf = (float) sqrt(2.0f * PI * sigma * sigma); + #else + + tmpf = sqrtf(2.0f * PI * sigma * sigma); + #endif + + gaussianCoefs[i] = (float) exp( (-((float) (i-GAUSS_WIDTH)) * + ((float) (i-GAUSS_WIDTH))) / + (2.0f * sigma * sigma)) / + tmpf; + } else { + /* If sigma equals zero (all pixels have the same value) + * the above code will divide by zero. Therefore, we instead + * use boxcar smoothing. + */ + gaussianCoefs[i] = 1.0f / (1.0f + (2.0f * (float) GAUSS_WIDTH)); + } + } + + for(i=0;inumBins;i++) { + for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) { + if (((j+i) >= 0) && ((j+i) < robustHistogram->numBins)) { + robustHistogram->nums->vec.i32[j+i]+= + (gaussianCoefs[j+GAUSS_WIDTH] * + (float) robustHistogram->nums->vec.i32[j+i]); + } + } + } +} + +/****************************************************************************** + p_psArraySampleQuartiles() + This procedure calculates the upper and/or lower quartiles of the + data set. It is assumed that the data set is small enough that we + can sort all the data points and calculate the quartiles exactly. + *****************************************************************************/ +void p_psArraySampleQuartiles(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *newStruct) +{ + psVector *unsortedVector = NULL; + psVector *sortedVector = NULL; + int count = 0; + int ind = 0; + int i = 0; + + // return is we have already calculated both quartile points. + if ((!isnan(newStruct->sampleLQ)) && + (!isnan(newStruct->sampleUQ))) { + return; + } + + if (-1 == newStruct->nValues) { + p_psArrayNValues(myVector, maskVector, maskVal, newStruct); + } + + + unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues); + sortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues); + + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->vec.ui8[i])) { + unsortedVector->vec.f[count++] = maskVector->vec.f[i]; + } + } + psSort(sortedVector, unsortedVector); + } else { + psSort(sortedVector, myVector); + } + + if (newStruct->options & PS_STAT_SAMPLE_LQ) { + ind = 3 * (newStruct->nValues / 4); + newStruct->sampleUQ = sortedVector->vec.f[ind]; + } + + if (newStruct->options & PS_STAT_SAMPLE_UQ) { + ind = (newStruct->nValues / 4); + newStruct->sampleLQ = sortedVector->vec.f[ind]; + } + + psVectorFree(unsortedVector); + psVectorFree(sortedVector); +} + + +/****************************************************************************** + p_psArrayRobustStats(): this procedure calcualtes a variety of robust + stat measures: + PS_STAT_ROBUST_MEAN + PS_STAT_ROBUST_MEAN_NVALUES + PS_STAT_ROBUST_MEDIAN + PS_STAT_ROBUST_MEDIAN_NVALUES + PS_STAT_ROBUST_MODE + PS_STAT_ROBUST_MODE_NVALUES + PS_STAT_ROBUST_STDEV + I have included all that computation in a single function, as opposed to + breaking it across several functions for one primary reason: the all + require the same basic initial processing steps (calculate the histogram, + etc.) however there is no currently defined means, in the SDRS, to + specify this computation as a separate step. If the above robust stat + measures were calcualted in separate functiosn, then much of the initial + processing would be duplicated. + *****************************************************************************/ +void p_psArrayRobustStats(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *newStruct) { psHistogram *robustHistogram = NULL; float binSize = 0.0; - - // if (isnan(myVector->robustLQ) || - // isnan(myVector->robustUQ)) { - // p_psArrayRobustQuartiles(myVector, maskVector, maskVal, newStruct); - // } - // binSize = ((myVector->robustUQ - myVector->robustLQ) / 1.34) / 10.0; - + float sigmaE = 0.0; + int LQBinNum = -1; + int UQBinNum = -1; + int i = 0; // Loop index variable. + int maxBinNum = 0; + int maxBinCount = 0; + float dL = 0.0; + + // NOTE: The SDRS states that the sample quartiles must be used to + // determine the initial bin sizes. However, the sample quartiles are + // calculated based on a full sort of the data set, regardless of the + // size of the data set. We should consult with IfA to ensure that this + // is really required. + + if (0.0 == newStruct->sampleUQ) { + newStruct->options = newStruct->options | PS_STAT_SAMPLE_UQ; + p_psArraySampleQuartiles(myVector, + maskVector, + maskVal, + newStruct); + } + + if (isnan(newStruct->sampleLQ)) { + newStruct->options = newStruct->options | PS_STAT_SAMPLE_LQ; + p_psArraySampleQuartiles(myVector, + maskVector, + maskVal, + newStruct); + } + + // Compute the initial bin size of the robust histogram. + sigmaE = (newStruct->sampleUQ - newStruct->sampleLQ) / 1.34f; + binSize = sigmaE / 10.0f; + + // Detemine minimum and maximum values in the data vector. + if (isnan(newStruct->min)) { + p_psArrayMin(myVector, maskVector, maskVal, newStruct); + } + if (isnan(newStruct->max)) { + p_psArrayMax(myVector, maskVector, maskVal, newStruct); + } + + // Create the histogram structure. robustHistogram = psHistogramAlloc(newStruct->min, newStruct->max, binSize); - // p_psArrayNValues(myVector, maskVector, maskVal, newStruct); - // robustHistogram = psGetArrayHistogram(robustHistogram, myVector); - // p_psArraySmooth(robustHistogram, (binSize / 4.0)); - // dL = (myVector->robustUQ - myVector->robustLQ) / 8.0; - - - // BROAD: Calculate the Robust Median - // Determine the LQ of the distribution. - // Determine the UQ of the distribution. - // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0. - // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4 - - // Find the bin with the peak value between LQ and UQ (the MODE) - // dL = (UQ - LQ) / 8 + // Populate the histogram arrat. + robustHistogram = psGetArrayHistogram(robustHistogram, myVector); + + // Smooth the histogram. + p_psArraysmoothHistGaussian(robustHistogram, sigmaE/4.0f); + + LQBinNum = -1; + UQBinNum = -1; + for (i=0;inumBins;i++) { + if ((robustHistogram->nums->vec.i32[i] <= newStruct->sampleLQ) && + (newStruct->sampleLQ <= robustHistogram->nums->vec.i32[i])) { + LQBinNum = i; + } + + if ((robustHistogram->nums->vec.i32[i] <= newStruct->sampleUQ) && + (newStruct->sampleUQ <= robustHistogram->nums->vec.i32[i])) { + UQBinNum = i; + } + } + + // Determine the bin with the peak value in the range LQ to UQ. + maxBinNum = LQBinNum; + maxBinCount = robustHistogram->nums->vec.i32[maxBinNum]; + for (i=LQBinNum;i<=UQBinNum;i++) { + if (robustHistogram->nums->vec.i32[i] > maxBinCount) { + maxBinNum = i; + maxBinCount = robustHistogram->nums->vec.i32[i]; + } + } + + dL = (newStruct->robustUQ - newStruct->robustLQ) / 8.0; + // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL + // What algorithm should I use to do this? + // The resulting fit parameters are the robust mean, mean_r, and sigma -} - - -float p_psArrayXXX(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal) -{ - printf("ERROR: Don't call me: p_psArrayXXX()\n"); - exit(1); + // What is the mean_r? + + if (newStruct->options & PS_STAT_ROBUST_MEAN) { + newStruct->robustMean = 0.0; + } + if (newStruct->options & PS_STAT_ROBUST_MEAN_NVALUES) { + newStruct->robustMeanNvalues = 0.0; + } + if (newStruct->options & PS_STAT_ROBUST_MEDIAN) { + newStruct->robustMedian = 0.0; + } + if (newStruct->options & PS_STAT_ROBUST_MEDIAN_NVALUES) { + newStruct->robustMedianNvalues = 0.0; + } + if (newStruct->options & PS_STAT_ROBUST_MODE) { + newStruct->robustMode = maxBinNum; + } + if (newStruct->options & PS_STAT_ROBUST_MODE_NVALUES) { + newStruct->robustModeNvalues = 0.0; + } + if (newStruct->options & PS_STAT_ROBUST_STDEV) { + newStruct->robustStdev = 0.0; + } + if (newStruct->options & PS_STAT_ROBUST_UQ) { + newStruct->robustUQ = 0.0; + } + if (newStruct->options & PS_STAT_ROBUST_LQ) { + newStruct->robustLQ = 0.0; + } } @@ -385,190 +610,197 @@ } countFloat = (float) countInt; + #ifdef DARWIN - newStruct->sampleStdev = sqrt( - #else - newStruct->sampleStdev = sqrtf( - #endif - (sumSquares-(sumDiffs * sumDiffs/countFloat))/ (countFloat-1)); - } - - /****************************************************************************** - *****************************************************************************/ - void p_psArraySampleQuartiles(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *newStruct) - { - psVector *unsortedVector = NULL; - psVector *sortedVector = NULL; - int count = 0; - int ind = 0; - int i = 0; - - // return is we have already calculated both quartile points. - if ((!isnan(newStruct->sampleLQ)) && - (!isnan(newStruct->sampleUQ))) { - return; - } - - if (-1 == newStruct->nValues) { - p_psArrayNValues(myVector, maskVector, maskVal, newStruct); - } - - - unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues); - sortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues); - - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->vec.ui8[i])) { - unsortedVector->vec.f[count++] = maskVector->vec.f[i]; - } - } - psSort(sortedVector, unsortedVector); - } else { - psSort(sortedVector, myVector); - } - - ind = 3 * (newStruct->nValues / 4); - newStruct->sampleUQ = sortedVector->vec.f[ind]; - ind = (newStruct->nValues / 4); - newStruct->sampleLQ = sortedVector->vec.f[ind]; - - psVectorFree(unsortedVector); - psVectorFree(sortedVector); - } - - /****************************************************************************** - - NOTE: The current strategy is to implement everything assuming that all - input data is of type PS_TYPE_FLOAT. Once the basic code is in place, - we will macro-ize everything and add PS_TYPE_UINT16 and PS_TYPE_DOUBLE. - - *****************************************************************************/ - psStats *psArrayStats(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) - { - psStats *newStruct = NULL; - - if (myVector == NULL) { - psAbort(__func__, - "Input data array (myVector) was NULL."); - } - - if (myVector->type.type != PS_TYPE_FLOAT) { - psAbort(__func__, - "Only data type PS_TYPE_FLOAT is currently supported."); - } - - if (maskVector != NULL) { - if (myVector->n != maskVector->n) { - psAbort(__func__, - "Vector data and vector mask are of different sizes."); - } - if (maskVector->type.type != PS_TYPE_UINT8) { - psAbort(__func__, - "Vector mask must be type PS_TYPE_UINT8"); - } - } - newStruct = psStatsAlloc(stats->options); - - // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_MEAN) { - p_psArraySampleMean(myVector, maskVector, maskVal, newStruct); - } - - // ************************************************************************ - if (stats->options & PS_STAT_MAX) { - p_psArrayMax(myVector, maskVector, maskVal, newStruct); - } - - // ************************************************************************ - if (stats->options & PS_STAT_MIN) { - p_psArrayMin(myVector, maskVector, maskVal, newStruct); - } - - // ************************************************************************ - if (stats->options & PS_STAT_NVALUES) { - p_psArrayNValues(myVector, maskVector, maskVal, newStruct); - } - - // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_MEDIAN) { - p_psArraySampleMedian(myVector, maskVector, maskVal, newStruct); - } - - // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_STDEV) { - p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct); - } - - // ************************************************************************ - if ((stats->options & PS_STAT_SAMPLE_UQ) || - (stats->options & PS_STAT_SAMPLE_LQ)) { - p_psArraySampleQuartiles(myVector, maskVector, maskVal, newStruct); - } - - - - - if (stats->options & PS_STAT_ROBUST_MEAN) { - newStruct->robustMean = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_ROBUST_MEAN_NVALUES) { - newStruct->robustMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - newStruct->robustMedian = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_ROBUST_MEDIAN_NVALUES) { - newStruct->robustMedianNvalues = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_ROBUST_MODE) { - newStruct->robustMode = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_ROBUST_MODE_NVALUES) { - newStruct->robustModeNvalues = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_ROBUST_STDEV) { - newStruct->robustStdev = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if ((stats->options & PS_STAT_ROBUST_UQ) || - (stats->options & PS_STAT_ROBUST_LQ)) { - newStruct->robustLQ = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_CLIPPED_MEAN) { - newStruct->clippedMean = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_CLIPPED_MEAN_NVALUES) { - newStruct->clippedMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_CLIPPED_MEAN_NSIGMA) { - newStruct->clipSigma = p_psArrayXXX(myVector, maskVector, maskVal); - } - - if (stats->options & PS_STAT_CLIPPED_STDEV) { - newStruct->clippedStdev = p_psArrayXXX(myVector, maskVector, maskVal); - } - - - // OLD CODE: Should we check for an unknown option? - // default: - // psAbort(__func__, "Unknown options 0x%x.\n", stats->options); - - return(newStruct); - } + newStruct->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs * + sumDiffs/countFloat))/ (countFloat-1)); + #else + + newStruct->sampleStdev = sqrtf( (sumSquares-(sumDiffs * + sumDiffs/countFloat))/ (countFloat-1)); + #endif +} + +void p_psArrayClippedStats(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *newStruct) +{ + int i = 0; + int j = 0; + float clippedMean = 0.0; + float clippedStdev = 0.0; + psVector *tmpMask = NULL; + + if (!((CLIPPED_NUM_ITER_LB <= newStruct->clipIter ) && + (newStruct->clipIter <= CLIPPED_NUM_ITER_UB))) { + psAbort(__func__, "Unallowed value for clipIter (%d).\n", + newStruct->clipIter); + } + + if (!((CLIPPED_SIGMA_LB <= newStruct->clipSigma ) && + (newStruct->clipSigma <= CLIPPED_SIGMA_UB))) { + psAbort(__func__, "Unallowed value for clipSigma (%f).\n", + newStruct->clipSigma); + } + + tmpMask = psVectorAlloc(maskVector->type.type, myVector->nalloc); + + tmpMask->n = maskVector->n; + for (i=0;in;i++) { + tmpMask->vec.ui8[i] = maskVector->vec.ui8[i]; + } + + // 1. Compute the sample median. + p_psArraySampleMedian(myVector, maskVector, maskVal, newStruct); + + // 2. Compute the sample standard deviation. + p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct); + + // 3. Use the sample median as the first estimator of the mean X. + clippedMean = newStruct->sampleMean; + + // 4. Use the sample stdev as the first estimator of the mean stdev. + clippedStdev = newStruct->sampleStdev; + + // 5. Repeat N times: + + for (i=0;iclipIter;i++) { + for (j=0;jn;j++) { + // a) Exclude all values x_i for which |x_i - x| > K * stdev + if ( fabs(myVector->vec.f[j] - clippedMean) > + (newStruct->clipSigma * clippedStdev)) { + tmpMask->vec.ui8[i] = 0xff; + } + // b) compute new mean and stdev + // GUS: I should probably create a new struct here since the + // following calls will overwrite any old values in sampleMean. + p_psArraySampleMedian(myVector, tmpMask, maskVal, newStruct); + p_psArraySampleStdev(myVector, tmpMask, maskVal, newStruct); + + // c) Use the new mean for x + clippedMean = newStruct->sampleMean; + + // d) Use the new stdev for stdev + clippedStdev = newStruct->sampleStdev; + } + } + + // 7. The last calcuated value of x is the cliped mean. + if (newStruct->options & PS_STAT_CLIPPED_MEAN) { + newStruct->clippedMean = clippedMean; + } + + // 8. The last calcuated value of stdev is the cliped stdev. + if (newStruct->options & PS_STAT_CLIPPED_STDEV) { + newStruct->clippedStdev = clippedStdev; + } + + if (newStruct->options & PS_STAT_CLIPPED_MEAN_NVALUES) { + p_psArrayNValues(myVector, tmpMask, maskVal, newStruct); + } + + if (newStruct->options & PS_STAT_CLIPPED_MEAN_NSIGMA) { + // GUS: What to do here? + } + + psVectorFree(tmpMask); +} + + +/****************************************************************************** + + NOTE: The current strategy is to implement everything assuming that all + input data is of type PS_TYPE_FLOAT. Once the basic code is in place, + we will macro-ize everything and add PS_TYPE_UINT16 and PS_TYPE_DOUBLE. + + *****************************************************************************/ +psStats *psArrayStats(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + psStats *newStruct = NULL; + + if (myVector == NULL) { + psAbort(__func__, + "Input data array (myVector) was NULL."); + } + + if (myVector->type.type != PS_TYPE_FLOAT) { + psAbort(__func__, + "Only data type PS_TYPE_FLOAT is currently supported."); + } + if (maskVector != NULL) { + if (myVector->n != maskVector->n) { + psAbort(__func__, + "Vector data and vector mask are of different sizes."); + } + if (maskVector->type.type != PS_TYPE_UINT8) { + psAbort(__func__, + "Vector mask must be type PS_TYPE_UINT8"); + } + } + newStruct = psStatsAlloc(stats->options); + + // ************************************************************************ + if (stats->options & PS_STAT_SAMPLE_MEAN) { + p_psArraySampleMean(myVector, maskVector, maskVal, newStruct); + } + + // ************************************************************************ + if (stats->options & PS_STAT_MAX) { + p_psArrayMax(myVector, maskVector, maskVal, newStruct); + } + + // ************************************************************************ + if (stats->options & PS_STAT_MIN) { + p_psArrayMin(myVector, maskVector, maskVal, newStruct); + } + + // ************************************************************************ + if (stats->options & PS_STAT_NVALUES) { + p_psArrayNValues(myVector, maskVector, maskVal, newStruct); + } + + // ************************************************************************ + if (stats->options & PS_STAT_SAMPLE_MEDIAN) { + p_psArraySampleMedian(myVector, maskVector, maskVal, newStruct); + } + + // ************************************************************************ + if (stats->options & PS_STAT_SAMPLE_STDEV) { + p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct); + } + + // ************************************************************************ + if ((stats->options & PS_STAT_SAMPLE_UQ) || + (stats->options & PS_STAT_SAMPLE_LQ)) { + p_psArraySampleQuartiles(myVector, maskVector, maskVal, newStruct); + } + + if ((stats->options & PS_STAT_ROBUST_MEAN) || + (stats->options & PS_STAT_ROBUST_MEAN_NVALUES) || + (stats->options & PS_STAT_ROBUST_MEDIAN) || + (stats->options & PS_STAT_ROBUST_MEDIAN_NVALUES) || + (stats->options & PS_STAT_ROBUST_MODE) || + (stats->options & PS_STAT_ROBUST_MODE_NVALUES) || + (stats->options & PS_STAT_ROBUST_STDEV) || + (stats->options & PS_STAT_ROBUST_UQ) || + (stats->options & PS_STAT_ROBUST_LQ)) { + p_psArrayClippedStats(myVector, maskVector, maskVal, newStruct); + } + + if ((stats->options & PS_STAT_CLIPPED_MEAN) || + (stats->options & PS_STAT_CLIPPED_MEAN_NVALUES) || + (stats->options & PS_STAT_CLIPPED_MEAN_NSIGMA) || + (stats->options & PS_STAT_CLIPPED_STDEV)) { + p_psArrayClippedStats(myVector, maskVector, maskVal, newStruct); + } + + // OLD CODE: Should we check for an unknown option? + // default: + // psAbort(__func__, "Unknown options 0x%x.\n", stats->options); + + return(newStruct); +}