Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 831) +++ trunk/psLib/src/math/psStats.c (revision 889) @@ -85,6 +85,12 @@ float binSize = 0.0; + if (n == 0) { + psAbort(__func__, "psHistogram requested with 0 bins"); + } + newHist = (psHistogram *) psAlloc(sizeof(psHistogram)); newHist->bounds = psVectorAlloc(n+1, PS_TYPE_F32); + newHist->bounds->n = newHist->bounds->nalloc; + binSize = (upper - lower) / (float) n; for (i=0;inums = psVectorAlloc(n, PS_TYPE_S32); + newHist->nums->n = newHist->nums->nalloc; + for (i=0;inums->n;i++) { + newHist->nums->data.S32[i] = 0; + } newHist->minNum = 0; newHist->maxNum = 0; @@ -99,5 +109,6 @@ } - +// When this is called, the number of data elements in bounds +// is n. Therefore, the number of bins is n-1. psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds) { @@ -107,8 +118,13 @@ newHist = (psHistogram *) psAlloc(sizeof(psHistogram)); newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32); + newHist->bounds->n = newHist->bounds->nalloc; for (i=0;in;i++) { newHist->bounds->data.F32[i] = bounds->data.F32[i]; } newHist->nums = psVectorAlloc((bounds->n)-1, PS_TYPE_S32); + newHist->nums->n = newHist->nums->nalloc; + for (i=0;inums->n;i++) { + newHist->nums->data.S32[i] = 0; + } newHist->minNum = 0; @@ -119,6 +135,8 @@ } -void psHistogramFree(psHistogram *restrict myHist) -{ +void psHistogramFree(psHistogram *myHist) +{ + psVectorFree(myHist->bounds); + psVectorFree(myHist->nums); psFree(myHist); } @@ -150,15 +168,51 @@ numBins = out->nums->n; - for (i=0;in;i++) { - // Check if this pixel is masked, and if so, skip it. - if (!(mask->data.S32[i] & maskVal)) { - // Check if this pixel is below the minimum value, and if so - // count it, then skip it. + + if (mask != NULL) { + for (i=0;in;i++) { + // Check if this pixel is masked, and if so, skip it. + if (!(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]) { + 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]) { + out->maxNum++; + } else { + // If this is a uniform histogram, determining the correct + // number is trivial. + 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); + (out->nums->data.S32[binNum])++; + // If this is a non-uniform histogram, determining the correct + // bin number requires a bit more work. + } else { + // GUS: This is slow. Put a smarter algorithm here to + // find the correct bin number (bin search, probably) + for (j=0;j<(out->bounds->n)-1;j++) { + if ((out->bounds->data.S32[j] <= in->data.F32[i]) && + (in->data.F32[i] <= out->bounds->data.S32[j+1])) { + (out->nums->data.S32[j])++; + } + } + } + } + } + } + } else { + for (i=0;in;i++) { if (in->data.F32[i] < out->bounds->data.F32[0]) { + // Check if this pixel is below the minimum value, and if so + // count it, then skip it. out->minNum++; - + } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { // 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]) { out->maxNum++; } else { @@ -171,4 +225,5 @@ binSize); (out->nums->data.S32[binNum])++; + // If this is a non-uniform histogram, determining the correct // bin number requires a bit more work. @@ -185,4 +240,5 @@ } } + } return(out); @@ -471,4 +527,5 @@ unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); unsortedVector->n = unsortedVector->nalloc; + sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); sortedVector->n = sortedVector->nalloc; @@ -494,5 +551,5 @@ if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = maskVector->data.F32[i]; + unsortedVector->data.F32[count++] = myVector->data.F32[i]; } } @@ -504,10 +561,10 @@ for (i=0;in;i++) { if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count++] = maskVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - unsortedVector->data.F32[i] = maskVector->data.F32[i]; + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + unsortedVector->data.F32[i] = myVector->data.F32[i]; } } @@ -596,25 +653,27 @@ 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; - } + /* GUS: When IfA provides the algorithm, insert is here. + if (myVector->n > stats->sampleLimit) { + psStats *stats2 = NULL; + // 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; + } + */ // Determine how many data points fit inside this min/max range @@ -640,5 +699,5 @@ (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = maskVector->data.F32[i]; + unsortedVector->data.F32[count++] = myVector->data.F32[i]; } } @@ -647,5 +706,5 @@ if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = maskVector->data.F32[i]; + unsortedVector->data.F32[count++] = myVector->data.F32[i]; } } @@ -657,10 +716,10 @@ for (i=0;in;i++) { if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count++] = maskVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - unsortedVector->data.F32[i] = maskVector->data.F32[i]; + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + unsortedVector->data.F32[i] = myVector->data.F32[i]; } } @@ -679,4 +738,5 @@ psVectorFree(unsortedVector); psVectorFree(sortedVector); + // GUS: This is the } @@ -691,5 +751,5 @@ 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 + breaking it across several functions for one primary reason: they all require the same basic initial processing steps (calculate the histogram, etc.) however there is no currently defined means, in the SDRS, to @@ -719,8 +779,7 @@ // is really required. - if (0.0 == stats->sampleUQ) { - - // GUS: fix this - // stats->options = stats->options | PS_STAT_SAMPLE_UQ; + if (isnan(stats->sampleUQ) || + isnan(stats->sampleLQ)) { + stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE; p_psVectorSampleQuartiles(myVector, maskVector, @@ -729,13 +788,4 @@ } - if (isnan(stats->sampleLQ)) { - // GUS: fix this - // stats->options = stats->options | PS_STAT_SAMPLE_LQ; - p_psVectorSampleQuartiles(myVector, - maskVector, - maskVal, - stats); - } - // Compute the initial bin size of the robust histogram. sigmaE = (stats->sampleUQ - stats->sampleLQ) / 1.34f; @@ -753,6 +803,8 @@ robustHistogram = psHistogramAlloc(stats->min, stats->max, - binSize); - // Populate the histogram arrat. + 1 + (int) ((stats->max - stats->min) / binSize)); + + + // Populate the histogram array. // GUS: fix this // robustHistogram = psHistogramVector(robustHistogram, myVector); @@ -811,4 +863,7 @@ stats->robustLQ = 0.0; } + stats->robustNfit = 0.0; + stats->robustN50 = 0.0; + psHistogramFree(robustHistogram); } @@ -836,5 +891,5 @@ if (0 != isnan(stats->sampleMean)) { - psAbort(__func__, "stats->sampleMean == NAN"); + p_psVectorSampleMean(myVector, maskVector, maskVal, stats); } mean = stats->sampleMean; @@ -920,12 +975,17 @@ } - tmpMask = psVectorAlloc(myVector->nalloc, maskVector->type.type); - - tmpMask->n = maskVector->n; - for (i=0;in;i++) { - tmpMask->data.U8[i] = maskVector->data.U8[i]; + + tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8); + tmpMask->n = myVector->n; + + if (maskVector != NULL) { + for (i=0;in;i++) { + tmpMask->data.U8[i] = maskVector->data.U8[i]; + } } // 1. Compute the sample median. + // GUS: This seems odd. Verify with IfA that we want to calculate the + // median here, not the mean. p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); @@ -942,7 +1002,9 @@ for (i=0;iclipIter;i++) { + // printf("p_psVectorClippedStats(): Iteration %d (%f %f)\n", i, + // clippedMean, clippedStdev); for (j=0;jn;j++) { // a) Exclude all values x_i for which |x_i - x| > K * stdev - if ( fabs(myVector->data.F32[j] - clippedMean) > + if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { tmpMask->data.U8[i] = 0xff; @@ -994,6 +1056,8 @@ if (in->type.type != PS_TYPE_F32) { psAbort(__func__, - "Only data type PS_TYPE_F32 is currently supported."); - } + "Only data type PS_TYPE_F32 is currently supported (0x%x).", + in->type.type); + } + if (mask != NULL) { if (in->n != mask->n) { @@ -1018,5 +1082,8 @@ // ************************************************************************ + // GUS: The Stdev calculation requires the mean. Should we assume the + // mean has already been calculated? Or should we always calculate it? if (stats->options & PS_STAT_SAMPLE_STDEV) { + p_psVectorSampleMean(in, mask, maskVal, stats); p_psVectorSampleStdev(in, mask, maskVal, stats); }