Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 2406) +++ trunk/psLib/src/math/psStats.c (revision 2411) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.98 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-23 19:35:30 $ + * @version $Revision: 1.99 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-11-24 00:05:54 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -128,10 +128,9 @@ NOTE: it is assumed that any call to these statistical functions will have been preceded by a call to the psVectorStats() function. Various sanity tests -will only be performed in psVectorStats(). - -XXX: Should we perform the sanity checks in each routine anyway? +will only be performed in psVectorStats(). Should we perform the sanity +checks in each routine anyway? XXX: For many of these private stats routines, what should be done if there -are now acceptable elements in the input vector (if no elements lie within +are no acceptable elements in the input vector (if no elements lie within range, or there are no unmasked elements, or the input vector is NULL)? Currently we set the value to NAN. @@ -931,10 +930,22 @@ /***************************************************************************** -psNormalizeVectorRange##TYPE(myData, low, high): this is a private function which -normalizes the elements of a vector to a range between low and high. - -XXX: Should we make the arguments psScalars? - -XXX: Should we make the arguments PS_TYPE_F64? +These macros and functions define the following functions: + + p_psNormalizeVectorRange(myData, low, high) + +That assumes that the low/high arguments are PS_TYPE_F64; the vector myData +can be of any type. Arguments low/high will be converted to the appropriate +type and one of the type-specific functions below will be called: + + p_psNormalizeVectorRangeU8(myData, low, high) + p_psNormalizeVectorRangeU16(myData, low, high) + p_psNormalizeVectorRangeU32(myData, low, high) + p_psNormalizeVectorRangeU64(myData, low, high) + p_psNormalizeVectorRangeS8(myData, low, high) + p_psNormalizeVectorRangeS16(myData, low, high) + p_psNormalizeVectorRangeS32(myData, low, high) + p_psNormalizeVectorRangeS64(myData, low, high) + p_psNormalizeVectorRangeF32(myData, low, high) + p_psNormalizeVectorRangeF64(myData, low, high) *****************************************************************************/ #define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \ @@ -984,4 +995,49 @@ PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64) +void p_psNormalizeVectorRange(psVector* myData, + psF64 outLow, + psF64 outHigh) +{ + switch (myData->type.type) { + case PS_TYPE_U8: + p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh); + break; + case PS_TYPE_U16: + p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh); + break; + case PS_TYPE_U32: + p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh); + break; + case PS_TYPE_U64: + p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh); + break; + case PS_TYPE_S8: + p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh); + break; + case PS_TYPE_S16: + p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh); + break; + case PS_TYPE_S32: + p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh); + break; + case PS_TYPE_S64: + p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh); + break; + case PS_TYPE_F32: + p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh); + break; + case PS_TYPE_F64: + p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh); + break; + case PS_TYPE_C32: + case PS_TYPE_C64: + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + "Unallowable operation: %s has incorrect type.", + myData); + break; + } +} + /****************************************************************************** p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes @@ -1053,6 +1109,4 @@ XXX: After you fit the polynomial, solve for X analytically. - -XXX: Check for errors in psLib routines that we call. *****************************************************************************/ float fitQuadraticSearchForYThenReturnX(psVector *xVec, @@ -1104,4 +1158,10 @@ // polynomial, looking for the value x such that f(x) = yVal tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal); + if (isnan(tmpFloat)) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN); + return(NAN); + } } else { @@ -1135,10 +1195,6 @@ 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: they 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. +breaking it across several functions for one primary reason: they all require +the same basic initial processing steps (calculate the histogram, etc.). Inputs @@ -1185,5 +1241,11 @@ // by computing the clipped standard deviation of the vector, and dividing // that by 10.0; - p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats); + int rc = p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats); + if (rc != 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_ROBUST_STATS_CLIPPED_STATS); + return(1); + } binSize = tmpStats->clippedStdev / 10.0f; @@ -1213,4 +1275,5 @@ return(0); } + // Determine minimum and maximum values in the data vector. if (isnan(stats->min)) { @@ -1221,8 +1284,10 @@ } } - - if (isnan(stats->max)) { - if (0 != p_psVectorMax(myVector, maskVector, maskVal, stats)) {} + if (0 != p_psVectorMax(myVector, maskVector, maskVal, stats)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: p_psVectorMin(): p_psVectorMax() reported a NAN mean.\n"); + return(1); + } } @@ -1594,5 +1659,5 @@ psS32 i = 0; // Loop index variable - float binSize = 0.0; // Histogram bin size + float binSize = 0.0; // Histogram bin size psS32 binNum = 0; // A temporary bin number psS32 numBins = 0; // The total number of bins @@ -1608,16 +1673,14 @@ } - numBins = out->nums->n; 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 ((mask == NULL) || + ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) { if (inF32->data.F32[i] < out->bounds->data.F32[0]) { + // If this pixel is below minimum value, count it, then skip. out->minNum++; - // Check if this pixel is above the maximum value, and if so - // count it, then skip it. } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) { + // If this pixel is above maximum value, count it, then skip. out->maxNum++; } else { @@ -1629,6 +1692,6 @@ // XXX: This if-statement really shouldn't be necessary. - // However, due to numerical lack of precision, we occasionally - // produce a binNum outside the range of bins. + // However, due to numerical lack of precision, we + // occasionally produce a binNum outside the range. if (binNum >= out->nums->n) { binNum = out->nums->n - 1; @@ -1637,11 +1700,12 @@ (out->nums->data.U32[binNum])++; - // If this is a non-uniform histogram, determining the correct - // bin number requires a bit more work. } else { + // If this is a non-uniform histogram, determining the + // correct bin number requires a bit more work. tmpScalar.data.F32 = inF32->data.F32[i]; tmp = p_psVectorBinDisect(out->bounds, &tmpScalar); if (tmp < 0) { - // XXX: Generate warning message + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: psVectorHistogram(): element outside histogram bounds.\n"); } else { (out->nums->data.U32[tmp])++; @@ -1666,5 +1730,5 @@ everything and put type support in the various stat functions. -XXX: Should the default data type be F64? Since we are buying Athlons... +XXX: Should the default data type be F64? Since we are buying Opterons... *****************************************************************************/ psVector* p_psConvertToF32(psVector* in) @@ -1744,18 +1808,30 @@ // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_MEAN) { - p_psVectorSampleMean(inF32, mask, maskVal, stats); + if (0 != p_psVectorSampleMean(inF32, mask, maskVal, stats)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n"); + } } // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_MEDIAN) { - p_psVectorSampleMedian(inF32, mask, maskVal, stats); + if (false == p_psVectorSampleMedian(inF32, mask, maskVal, stats)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: psVectorStats(): p_psVectorSampleMedian() returned an error.\n"); + } } // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_STDEV) { - p_psVectorSampleMean(inF32, mask, maskVal, stats); + if (0 != p_psVectorSampleMean(inF32, mask, maskVal, stats)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n"); + } p_psVectorSampleStdev(inF32, mask, maskVal, stats); } // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_QUARTILE) { - p_psVectorSampleQuartiles(inF32, mask, maskVal, stats); + if (false == p_psVectorSampleQuartiles(inF32, mask, maskVal, stats)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: psVectorStats(): p_psVectorSampleQuartiles() returned an error.\n"); + } } // Since the various robust stats quantities share much computation, they @@ -1765,5 +1841,6 @@ (stats->options & PS_STAT_ROBUST_MEDIAN) || (stats->options & PS_STAT_ROBUST_MODE) || - (stats->options & PS_STAT_ROBUST_STDEV) || (stats->options & PS_STAT_ROBUST_QUARTILE)) { + (stats->options & PS_STAT_ROBUST_STDEV) || + (stats->options & PS_STAT_ROBUST_QUARTILE)) { if (0 != p_psVectorRobustStats(inF32, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false,