Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 796) +++ trunk/psLib/src/math/psStats.c (revision 797) @@ -15,5 +15,5 @@ #include "float.h" -#define ROBUST_SIZE_THRESHOLD 10000 // Vectors that are large than this +#define ROBUST_SIZE_THRESHOLD 30000 // Vectors that are large than this // will use robust statistical methods. #define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar @@ -27,4 +27,9 @@ #define false 0 #define MYMAXFLOAT 1e99 + +void p_psVectorRobustStats(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats); /****************************************************************************** psStatsAlloc(): This routine must create a new psStats data structure. @@ -40,5 +45,5 @@ newStruct->sampleUQ = NAN; newStruct->sampleLQ = NAN; - newStruct->sampleLimit = NAN; + newStruct->sampleLimit = ROBUST_SIZE_THRESHOLD; newStruct->robustMean = NAN; newStruct->robustMedian = NAN; @@ -51,6 +56,6 @@ newStruct->clippedMean = NAN; newStruct->clippedStdev = NAN; - newStruct->clipSigma = NAN; - newStruct->clipIter = -1; + newStruct->clipSigma = 3.0; + newStruct->clipIter = 3.0; newStruct->min = NAN; newStruct->max = NAN; @@ -187,5 +192,5 @@ psVector *maskVector, unsigned int maskVal, - psStats *newStruct) + psStats *stats) { int i = 0; @@ -207,4 +212,10 @@ ****************************************************************************** *****************************************************************************/ + + +/****************************************************************************** + ASSUMPTION: the mean is always calculated exactly. Robust means are never + calculated in this routine. + *****************************************************************************/ void p_psVectorSampleMean(const psVector *restrict myVector, const psVector *restrict maskVector, @@ -368,22 +379,48 @@ /****************************************************************************** - This routine calculates the number of non-masked pixels in the vector. + This routine calculates the number of non-masked pixels in the vector. + that fall within the min/max range, if given. *****************************************************************************/ int p_psVectorNValues(const psVector *restrict myVector, const psVector *restrict maskVector, unsigned int maskVal, - psStats *newStruct) + psStats *stats) { int i = 0; int numData = 0; - - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->vec.ui8[i])) { - numData++; + float rangeMin = 0.0; + float rangeMax = 0.0; + + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->vec.ui8[i]) && + (rangeMin <= myVector->vec.f[i]) && + (myVector->vec.f[i] <= rangeMax)) { + numData++; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->vec.f[i]) && + (myVector->vec.f[i] <= rangeMax)) { + numData++; + } } } } else { - numData = myVector->n; + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->vec.ui8[i])) { + numData++; + } + } + } else { + numData = myVector->n; + } } return(numData); @@ -395,5 +432,5 @@ const psVector *restrict maskVector, unsigned int maskVal, - psStats *newStruct) + psStats *stats) { psVector *unsortedVector = NULL; @@ -401,8 +438,35 @@ int count = 0; int i = 0; - float median = 0.0; int nValues = 0; - - nValues = p_psVectorNValues(myVector, maskVector, maskVal, newStruct); + float rangeMin = 0.0; + float rangeMax = 0.0; + psStats *stats2 = NULL; + + + // Determine if the number of data points exceed a threshold which will + // cause to generate robust stats, as opposed to exact stats. + + if (myVector->n > stats->sampleLimit) { + // 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. + psStatsFree(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 maxed, IF the maskVector is not NULL> + nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); + + // Allocate temporary vectors for the data. unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, nValues); unsortedVector->n = unsortedVector->nalloc; @@ -410,28 +474,63 @@ sortedVector->n = sortedVector->nalloc; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->vec.ui8[i])) { - unsortedVector->vec.f[count++] = myVector->vec.f[i]; - - } - } - psSort(sortedVector, unsortedVector); + // Determine if we must only use data points within a min/max range. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + + // Store all non-masked data points within the min/max range + // into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->vec.ui8[i]) && + (rangeMin <= myVector->vec.f[i]) && + (myVector->vec.f[i] <= rangeMax)) { + unsortedVector->vec.f[count++] = maskVector->vec.f[i]; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->vec.f[i]) && + (myVector->vec.f[i] <= rangeMax)) { + unsortedVector->vec.f[count++] = maskVector->vec.f[i]; + } + } + } } else { - psSort(sortedVector, myVector); - } - + // Store all non-masked data points into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->vec.ui8[i])) { + unsortedVector->vec.f[count++] = maskVector->vec.f[i]; + } + } + } else { + for (i=0;in;i++) { + unsortedVector->vec.f[i] = maskVector->vec.f[i]; + } + } + } + // Sort the temporary vectors. + psSort(sortedVector, unsortedVector); + + // Calculate the median exactly. if (0 == (nValues % 2)) { - median = 0.5 * (sortedVector->vec.f[(nValues/2)-1] + - sortedVector->vec.f[nValues/2]); + stats->sampleMedian = 0.5 * (sortedVector->vec.f[(nValues/2)-1] + + sortedVector->vec.f[nValues/2]); } else { - median = sortedVector->vec.f[nValues/2]; - } - + stats->sampleMedian = sortedVector->vec.f[nValues/2]; + } + + // Free the temporary data structures. psVectorFree(unsortedVector); psVectorFree(sortedVector); - newStruct->sampleMedian = median; -} - +} + +/****************************************************************************** + This routine smoothes the data in the input robustHistogram with a + Gaussian of width sigma. + *****************************************************************************/ void p_psVectorsmoothHistGaussian(psHistogram *robustHistogram, float sigma) @@ -481,11 +580,11 @@ p_psVectorSampleQuartiles() 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. + data set. + *****************************************************************************/ void p_psVectorSampleQuartiles(const psVector *restrict myVector, const psVector *restrict maskVector, unsigned int maskVal, - psStats *newStruct) + psStats *stats) { psVector *unsortedVector = NULL; @@ -495,33 +594,87 @@ int i = 0; int nValues = 0; - - // return is we have already calculated both quartile points. - if ((!isnan(newStruct->sampleLQ)) && - (!isnan(newStruct->sampleUQ))) { + float rangeMin = 0.0; + float rangeMax = 0.0; + psStats *stats2 = NULL; + + // Determine if the number of data points exceed a threshold which will + // cause to generate robust stats, as opposed to exact stats. + if (myVector->n > stats->sampleLimit) { + // Calculate the robust quartiles. + stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE); + p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); + + // Store the robust quartiles into the sample quartile members. + stats->sampleUQ = stats2->robustUQ; + stats->sampleLQ = stats2->robustLQ; + + // Free temporary data buffers. + psStatsFree(stats2); + + // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. + stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; + return; } - nValues = p_psVectorNValues(myVector, maskVector, maskVal, newStruct); - + // Determine how many data points fit inside this min/max range + // and are not maxed, IF the maskVector is not NULL> + nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); + + // Allocate temporary vectors for the data. unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, nValues); + unsortedVector->n = unsortedVector->nalloc; sortedVector = psVectorAlloc(PS_TYPE_FLOAT, 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); + sortedVector->n = sortedVector->nalloc; + + // Determine if we must only use data points within a min/max range. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + // Store all non-masked data points within the min/max range + // into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->vec.ui8[i]) && + (rangeMin <= myVector->vec.f[i]) && + (myVector->vec.f[i] <= rangeMax)) { + unsortedVector->vec.f[count++] = maskVector->vec.f[i]; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->vec.f[i]) && + (myVector->vec.f[i] <= rangeMax)) { + unsortedVector->vec.f[count++] = maskVector->vec.f[i]; + } + } + } } else { - psSort(sortedVector, myVector); - } - + // Store all non-masked data points into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->vec.ui8[i])) { + unsortedVector->vec.f[count++] = maskVector->vec.f[i]; + } + } + } else { + for (i=0;in;i++) { + unsortedVector->vec.f[i] = maskVector->vec.f[i]; + } + } + } + + // Sort the temporary vectors. + psSort(sortedVector, unsortedVector); + + // Calculate the quartile points exactly. ind = 3 * (nValues / 4); - newStruct->sampleUQ = sortedVector->vec.f[ind]; + stats->sampleUQ = sortedVector->vec.f[ind]; ind = (nValues / 4); - newStruct->sampleLQ = sortedVector->vec.f[ind]; - + stats->sampleLQ = sortedVector->vec.f[ind]; + + // Free the temporary data structures. psVectorFree(unsortedVector); psVectorFree(sortedVector); @@ -530,10 +683,11 @@ /****************************************************************************** - p_psVectorRobustStats(): this procedure calcualtes a variety of robust + p_psVectorRobustStats(): this procedure calculates a variety of robust stat measures: - PS_STAT_ROBUST_MEAN + PS_STAT_ROBUST_MEAN PS_STAT_ROBUST_MEDIAN PS_STAT_ROBUST_MODE PS_STAT_ROBUST_STDEV + PS_STAT_ROBUST_QUARTILE I have included all that computation in a single function, as opposed to breaking it across several functions for one primary reason: the all @@ -547,5 +701,5 @@ const psVector *restrict maskVector, unsigned int maskVal, - psStats *newStruct) + psStats *stats) { psHistogram *robustHistogram = NULL; @@ -565,37 +719,38 @@ // is really required. - if (0.0 == newStruct->sampleUQ) { + if (0.0 == stats->sampleUQ) { + // GUS: fix this - // newStruct->options = newStruct->options | PS_STAT_SAMPLE_UQ; + // stats->options = stats->options | PS_STAT_SAMPLE_UQ; p_psVectorSampleQuartiles(myVector, maskVector, maskVal, - newStruct); - } - - if (isnan(newStruct->sampleLQ)) { + stats); + } + + if (isnan(stats->sampleLQ)) { // GUS: fix this - // newStruct->options = newStruct->options | PS_STAT_SAMPLE_LQ; + // stats->options = stats->options | PS_STAT_SAMPLE_LQ; p_psVectorSampleQuartiles(myVector, maskVector, maskVal, - newStruct); + stats); } // Compute the initial bin size of the robust histogram. - sigmaE = (newStruct->sampleUQ - newStruct->sampleLQ) / 1.34f; + sigmaE = (stats->sampleUQ - stats->sampleLQ) / 1.34f; binSize = sigmaE / 10.0f; // Detemine minimum and maximum values in the data vector. - if (isnan(newStruct->min)) { - p_psVectorMin(myVector, maskVector, maskVal, newStruct); - } - if (isnan(newStruct->max)) { - p_psVectorMax(myVector, maskVector, maskVal, newStruct); + if (isnan(stats->min)) { + p_psVectorMin(myVector, maskVector, maskVal, stats); + } + if (isnan(stats->max)) { + p_psVectorMax(myVector, maskVector, maskVal, stats); } // Create the histogram structure. - robustHistogram = psHistogramAlloc(newStruct->min, - newStruct->max, + robustHistogram = psHistogramAlloc(stats->min, + stats->max, binSize); // Populate the histogram arrat. @@ -610,11 +765,11 @@ UQBinNum = -1; for (i=0;inums->n;i++) { - if ((robustHistogram->nums->vec.i32[i] <= newStruct->sampleLQ) && - (newStruct->sampleLQ <= robustHistogram->nums->vec.i32[i])) { + if ((robustHistogram->nums->vec.i32[i] <= stats->sampleLQ) && + (stats->sampleLQ <= robustHistogram->nums->vec.i32[i])) { LQBinNum = i; } - if ((robustHistogram->nums->vec.i32[i] <= newStruct->sampleUQ) && - (newStruct->sampleUQ <= robustHistogram->nums->vec.i32[i])) { + if ((robustHistogram->nums->vec.i32[i] <= stats->sampleUQ) && + (stats->sampleUQ <= robustHistogram->nums->vec.i32[i])) { UQBinNum = i; } @@ -631,5 +786,5 @@ } - dL = (newStruct->robustUQ - newStruct->robustLQ) / 8.0; + dL = (stats->robustUQ - stats->robustLQ) / 8.0; // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL @@ -639,20 +794,20 @@ // What is the mean_r? - if (newStruct->options & PS_STAT_ROBUST_MEAN) { - newStruct->robustMean = 0.0; - } - - if (newStruct->options & PS_STAT_ROBUST_MEDIAN) { - newStruct->robustMedian = 0.0; - } - if (newStruct->options & PS_STAT_ROBUST_MODE) { - newStruct->robustMode = maxBinNum; - } - if (newStruct->options & PS_STAT_ROBUST_STDEV) { - newStruct->robustStdev = 0.0; - } - if (newStruct->options & PS_STAT_ROBUST_QUARTILE) { - newStruct->robustUQ = 0.0; - newStruct->robustLQ = 0.0; + if (stats->options & PS_STAT_ROBUST_MEAN) { + stats->robustMean = 0.0; + } + + if (stats->options & PS_STAT_ROBUST_MEDIAN) { + stats->robustMedian = 0.0; + } + if (stats->options & PS_STAT_ROBUST_MODE) { + stats->robustMode = maxBinNum; + } + if (stats->options & PS_STAT_ROBUST_STDEV) { + stats->robustStdev = 0.0; + } + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + stats->robustUQ = 0.0; + stats->robustLQ = 0.0; } } @@ -661,4 +816,7 @@ NOTE: This function assumes that p_psVectorMean() has already been called and the correct value is stored in stats->sampleMean. + + NOTE: the mean is always calculated exactly. Robust means are never + calculated in this routine. *****************************************************************************/ void p_psVectorSampleStdev(const psVector *restrict myVector, @@ -742,5 +900,5 @@ const psVector *restrict maskVector, unsigned int maskVal, - psStats *newStruct) + psStats *stats) { int i = 0; @@ -750,14 +908,14 @@ psVector *tmpMask = NULL; - if (!((CLIPPED_NUM_ITER_LB <= newStruct->clipIter ) && - (newStruct->clipIter <= CLIPPED_NUM_ITER_UB))) { + if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) && + (stats->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))) { + stats->clipIter); + } + + if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && + (stats->clipSigma <= CLIPPED_SIGMA_UB))) { psAbort(__func__, "Unallowed value for clipSigma (%f).\n", - newStruct->clipSigma); + stats->clipSigma); } @@ -770,22 +928,22 @@ // 1. Compute the sample median. - p_psVectorSampleMedian(myVector, maskVector, maskVal, newStruct); + p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); // 2. Compute the sample standard deviation. - p_psVectorSampleStdev(myVector, maskVector, maskVal, newStruct); + p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); // 3. Use the sample median as the first estimator of the mean X. - clippedMean = newStruct->sampleMean; + clippedMean = stats->sampleMean; // 4. Use the sample stdev as the first estimator of the mean stdev. - clippedStdev = newStruct->sampleStdev; + clippedStdev = stats->sampleStdev; // 5. Repeat N times: - for (i=0;iclipIter;i++) { + 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)) { + (stats->clipSigma * clippedStdev)) { tmpMask->vec.ui8[i] = 0xff; } @@ -793,23 +951,23 @@ // GUS: I should probably create a new struct here since the // following calls will overwrite any old values in sampleMean. - p_psVectorSampleMedian(myVector, tmpMask, maskVal, newStruct); - p_psVectorSampleStdev(myVector, tmpMask, maskVal, newStruct); + p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); + p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); // c) Use the new mean for x - clippedMean = newStruct->sampleMean; + clippedMean = stats->sampleMean; // d) Use the new stdev for stdev - clippedStdev = newStruct->sampleStdev; + clippedStdev = stats->sampleStdev; } } // 7. The last calcuated value of x is the cliped mean. - if (newStruct->options & PS_STAT_CLIPPED_MEAN) { - newStruct->clippedMean = clippedMean; + if (stats->options & PS_STAT_CLIPPED_MEAN) { + stats->clippedMean = clippedMean; } // 8. The last calcuated value of stdev is the cliped stdev. - if (newStruct->options & PS_STAT_CLIPPED_STDEV) { - newStruct->clippedStdev = clippedStdev; + if (stats->options & PS_STAT_CLIPPED_STDEV) { + stats->clippedStdev = clippedStdev; }