Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 2339) +++ trunk/psLib/src/math/psStats.c (revision 2340) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.90 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-11 20:13:57 $ + * @version $Revision: 1.91 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-11-12 19:08:17 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -29,4 +29,5 @@ #include "psVector.h" #include "psTrace.h" +#include "psLogMsg.h" #include "psError.h" #include "psStats.h" @@ -142,13 +143,6 @@ /****************************************************************************** p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the -mean of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - +mean of the input vector. If there was a problem with the mean calculation, +this routine sets stats->sampleMean to NAN. *****************************************************************************/ void p_psVectorSampleMean(const psVector* restrict myVector, @@ -220,12 +214,6 @@ /****************************************************************************** p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the -max of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL +max of the input vector. If there was a problem with the max calculation, +this routine sets stats->max to NAN. *****************************************************************************/ void p_psVectorMax(const psVector* restrict myVector, @@ -289,12 +277,6 @@ /****************************************************************************** p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the -minimum of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL +minimum of the input vector. If there was a problem with the min calculation, +this routine sets stats->min to NAN. *****************************************************************************/ void p_psVectorMin(const psVector* restrict myVector, @@ -357,14 +339,56 @@ /****************************************************************************** +p_psVectorNValues(myVector, maskVector, maskVal, stats): This routine returns +"true" if the inputPsVector as 1 or more valid elements (a valid element is an +unmasked element within the specifined min/max range). Otherwise, return +"false". + *****************************************************************************/ +bool p_psVectorCheckNonEmpty(const psVector* restrict myVector, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) +{ + PS_VECTOR_CHECK_NULL(myVector, -1); + PS_PTR_CHECK_NULL(stats, -1); + psS32 i = 0; // Loop index variable + + if (stats->options & PS_STAT_USE_RANGE) { + if (maskVector != NULL) { + for (i = 0; i < myVector->n; i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { + return(true); + } + } + } else { + for (i = 0; i < myVector->n; i++) { + if ((stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { + return(true); + } + } + } + } else { + if (maskVector != NULL) { + for (i = 0; i < myVector->n; i++) { + if (!(maskVal & maskVector->data.U8[i])) { + return(true); + } + } + } else { + if (myVector->n > 0) { + return(true); + } + } + } + return(false); +} + + +/****************************************************************************** p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the number of non-masked pixels in the vector that fall within the min/max range, if specified. -Inputs - myVector - maskVector - maskVal - stats -Returns - int *****************************************************************************/ psS32 p_psVectorNValues(const psVector* restrict myVector, @@ -412,12 +436,6 @@ /****************************************************************************** p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the -median of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL +median of the input vector. Returns true on success (including if there were +no valid input vector elements). *****************************************************************************/ bool p_psVectorSampleMedian(const psVector* restrict myVector, @@ -429,9 +447,9 @@ psVector* sortedVector = NULL; // Temporary vector psS32 i = 0; // Loop index variable - psS32 count = 0; // # of points in this mean? - psS32 nValues = 0; // # of points in vector + psS32 count = 0; + psS32 nValues = 0; // Determine how many data points fit inside this min/max range - // and are not masked, if the maskVector is not NULL> + // and are not masked, if the maskVector is not NULL. nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); @@ -439,7 +457,7 @@ // sampleMedian to? Should we generate an error? if (nValues <= 0) { - printf("WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n"); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n"); return(true); - stats->sampleMedian = 0; + stats->sampleMedian = NAN; } @@ -569,5 +587,4 @@ } else if (x.data.F32 >= lastBound) { jMin = robustHistogram->bounds->n - 1; - } @@ -575,9 +592,14 @@ if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { jMax = p_psVectorBinDisect(robustHistogram->bounds, &x); + if (jMax < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } } else if (x.data.F32 <= firstBound) { jMax = 0; } else if (x.data.F32 >= lastBound) { jMax = robustHistogram->bounds->n - 1; - } @@ -606,5 +628,5 @@ NULL *****************************************************************************/ -void p_psVectorSampleQuartiles(const psVector* restrict myVector, +bool p_psVectorSampleQuartiles(const psVector* restrict myVector, const psVector* restrict maskVector, psU32 maskVal, @@ -616,19 +638,14 @@ psS32 count = 0; // # of points in this mean. psS32 nValues = 0; // # data points - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this // Determine how many data points fit inside this min/max range - // and are not maxed, IF the maskVector is not NULL> + // 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(nValues, PS_TYPE_F32); - sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); // 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. @@ -637,5 +654,5 @@ for (i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { + (stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { unsortedVector->data.F32[count++] = myVector->data.F32[i]; } @@ -643,5 +660,5 @@ } else { for (i = 0; i < myVector->n; i++) { - if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { + if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { unsortedVector->data.F32[count++] = myVector->data.F32[i]; } @@ -665,7 +682,15 @@ // Sort the temporary vectors. - psVectorSort(sortedVector, unsortedVector); + sortedVector = psVectorSort(sortedVector, unsortedVector); + if (sortedVector == NULL) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_SAMPLE_MEDIAN_SORT_PROBLEM); + return(false); + } // Calculate the quartile points exactly. + // XXX: We should probably do a bin-midpoint if the quartile is not an + // integer. stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)]; stats->sampleLQ = sortedVector->data.F32[nValues / 4]; @@ -674,4 +699,5 @@ psFree(unsortedVector); psFree(sortedVector); + return(true); } @@ -700,6 +726,4 @@ float sumSquares = 0.0; // temporary variable float sumDiffs = 0.0; // temporary variable - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this // This procedure requires the mean. If it has not been already @@ -708,11 +732,18 @@ p_psVectorSampleMean(myVector, maskVector, maskVal, stats); } + // If the mean is NAN, then generate a warning and set the stdev to NAN. + if (0 != isnan(stats->sampleMean)) { + stats->sampleStdev = NAN; + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n"); + return; + } + mean = stats->sampleMean; - if (stats->options & PS_STAT_USE_RANGE) { if (maskVector != NULL) { for (i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { + (stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { diff = myVector->data.F32[i] - mean; sumSquares += (diff * diff); @@ -723,5 +754,6 @@ } else { for (i = 0; i < myVector->n; i++) { - if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { + if ((stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { diff = myVector->data.F32[i] - mean; sumSquares += (diff * diff); @@ -752,15 +784,17 @@ } } - if (countInt < 2) { + if (countInt <= 1) { stats->sampleStdev = NAN; - // XXX PS WARNING + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements. Setting stats->sampleStdev = NAN.\n"); } else { countFloat = (float)countInt; #ifdef DARWIN - stats->sampleStdev = (float)sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); + stats->sampleStdev = (float)sqrt((sumSquares - + (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); #else - stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); + stats->sampleStdev = sqrtf((sumSquares - + (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); #endif @@ -779,4 +813,6 @@ Returns 0 for success. + -1: error + -2: warning *****************************************************************************/ int p_psVectorClippedStats(const psVector* restrict myVector, @@ -794,8 +830,12 @@ // Ensure that stats->clipIter is within the proper range. - PS_INT_CHECK_RANGE(stats->clipIter,PS_CLIPPED_NUM_ITER_LB,PS_CLIPPED_NUM_ITER_UB,-1); + PS_INT_CHECK_RANGE(stats->clipIter, + PS_CLIPPED_NUM_ITER_LB, + PS_CLIPPED_NUM_ITER_UB, -1); // Ensure that stats->clipSigma is within the proper range. - PS_INT_CHECK_RANGE(stats->clipSigma,PS_CLIPPED_SIGMA_LB,PS_CLIPPED_SIGMA_UB,-1); + PS_INT_CHECK_RANGE(stats->clipSigma, + PS_CLIPPED_SIGMA_LB, + PS_CLIPPED_SIGMA_UB, -1); // We allocate a temporary mask vector since during the iterative @@ -803,5 +843,4 @@ // However, we do no want to modify the original mask vector. tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8); - tmpMask->n = myVector->n; // If we were called with a mask vector, then initialize the temporary @@ -814,7 +853,19 @@ // 1. Compute the sample median. p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); + if (isnan(stats->sampleMean)) { + psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleMean returned NAN\n"); + stats->clippedMean = NAN; + stats->clippedStdev = NAN; + return(-2); + } // 2. Compute the sample standard deviation. p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); + if (isnan(stats->sampleStdev)) { + psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleStdev returned NAN\n"); + stats->clippedMean = NAN; + stats->clippedStdev = NAN; + return(-2); + } // 3. Use the sample median as the first estimator of the mean X. @@ -836,16 +887,21 @@ tmpMask->data.U8[i] = 0xff; } - // b) compute new mean and stdev - p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); - p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); - - // c) Use the new mean for x + } + + // b) compute new mean and stdev + p_psVectorSampleMean(myVector, tmpMask, maskVal, stats); + p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); + + // If the new mean and stdev are not NAN, then we use them. + // Otherwise, keep the old ones, and exit the loop. + if (! (isnan(stats->sampleMean) || isnan(stats->sampleStdev))) { clippedMean = stats->sampleMean; - - // d) Use the new stdev for stdev clippedStdev = stats->sampleStdev; - } - - } + + // Exit loop + i = stats->clipIter; + } + } + stats->sampleMean = oldStanMean; stats->sampleStdev = oldStanStdev; @@ -865,89 +921,54 @@ /***************************************************************************** - *****************************************************************************/ -void p_psNormalizeVectorRangeF32(psVector* myData, - float outLow, - float outHigh) -{ - float min = PS_MAX_F32; - float max = -PS_MAX_F32; - psS32 i = 0; - - for (i = 0; i < myData->n; i++) { - if (myData->data.F32[i] < min) { - min = myData->data.F32[i]; - } - if (myData->data.F32[i] > max) { - max = myData->data.F32[i]; - } - } - - // Ensure that max!=min before we divide by (max-min) - if (FLT_EPSILON < fabs(max - min)) { - for (i = 0; i < myData->n; i++) { - myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * - (outHigh - outLow) / (max - min); - } - } else { - // XXX: PS_WARNING - for (i = 0; i < myData->n; i++) { - myData->data.F32[i] = outLow; - } - } - - for (i = 0; i < myData->n; i++) { - myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * (outHigh - outLow) / (max - min); - } -} -void p_psNormalizeVectorRangeF64(psVector* myData, - float outLow, - float outHigh) -{ - float min = PS_MAX_F32; - float max = -PS_MAX_F32; - psS32 i = 0; - - for (i = 0; i < myData->n; i++) { - if (myData->data.F32[i] < min) { - min = myData->data.F64[i]; - } - if (myData->data.F32[i] > max) { - max = myData->data.F64[i]; - } - } - - // Ensure that max!=min before we divide by (max-min) - if (FLT_EPSILON < fabs(max - min)) { - for (i = 0; i < myData->n; i++) { - myData->data.F64[i] = outLow + (myData->data.F64[i] - min) * (outHigh - outLow) / (max - min); - } - } else { - // XXX: PS_WARNING - for (i = 0; i < myData->n; i++) { - myData->data.F64[i] = outLow; - } - } -} - -/***************************************************************************** -psNormalizeVectorRange(myData, low, high): this is a private function which +psNormalizeVectorRange##TYPE(myData, low, high): this is a private function which normalizes the elements of a vector to a range between low and high. - *****************************************************************************/ -void psNormalizeVectorRange(psVector* myData, - float low, - float high) -{ - if (myData->type.type == PS_TYPE_F32) { - p_psNormalizeVectorRangeF32(myData, low, high); - } else if (myData->type.type == PS_TYPE_F64) { - p_psNormalizeVectorRangeF64(myData, low, high); - } else { - char* strType; - PS_TYPE_NAME(strType,myData->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psStats_NOT_F32_F64, - strType); - } -} + +XXX: Should we make the arguments psScalars? + *****************************************************************************/ +#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \ +void p_psNormalizeVectorRange##TYPE(psVector* myData, \ + ps##TYPE outLow, \ + ps##TYPE outHigh) \ +{ \ + ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \ + ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \ + psS32 i = 0; \ + \ + for (i = 0; i < myData->n; i++) { \ + if (myData->data.TYPE[i] < min) { \ + min = myData->data.TYPE[i]; \ + } \ + if (myData->data.TYPE[i] > max) { \ + max = myData->data.TYPE[i]; \ + } \ + } \ + \ + /* Ensure that max!=min before we divide by (max-min) */ \ + if (max != min) { \ + for (i = 0; i < myData->n; i++) { \ + myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \ + (outHigh - outLow) / (max - min)); \ + } \ + } else { \ + psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min). Setting all elements to min.\n"); + for (i = 0; i < myData->n; i++) + { + \ + myData->data.TYPE[i] = outLow; + \ + } \ + } \ +} \ + +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8) +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16) +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32) +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64) +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8) +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16) +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32) +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64) +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32) +PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64) /****************************************************************************** @@ -969,4 +990,5 @@ float getThisValue) { + PS_FLOAT_COMPARE(rangeLow, rangeHigh, NAN); psS32 numIterations = 0; float midpoint = 0.0; @@ -1005,4 +1027,6 @@ XXX: After you fit the polynomial, solve for X analytically. + +XXX: Check for errors in psLib routines that we call. *****************************************************************************/ float fitQuadraticSearchForYThenReturnX(psVector *xVec, @@ -1085,4 +1109,6 @@ Returns 0 on success. + +XXX: Check for errors in psLib routines that we call. *****************************************************************************/ int p_psVectorRobustStats(const psVector* restrict myVector, @@ -1256,5 +1282,5 @@ psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - psNormalizeVectorRange(robustHistogramVector, 0.0, 1.0); + psNormalizeVectorRangeF32(robustHistogramVector, 0.0, 1.0); for (i=0;in;i++) { myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); @@ -1282,6 +1308,9 @@ if (stats->options & PS_STAT_ROBUST_MEAN) { if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) { - printf("WARNING: the fitted Gaussian has more than 10% error for the mean.\n"); - printf("Using the calculated mean (calc, fit) is (%f, %f)\n", myMean, myParams->data.F32[0]); + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: the fitted Gaussian has more than 10% error for the mean.\n"); + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: Using the calculated mean (calc, fit) is (%f, %f)\n", + myMean, myParams->data.F32[0]); stats->robustMean = myMean; } else { @@ -1296,6 +1325,9 @@ if (stats->options & PS_STAT_ROBUST_STDEV) { if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) { - printf("WARNING: the fitted Gaussian has more than 10% error for the stdev.\n"); - printf("Using the calculated stdev (calc, fit) is (%f, %f)\n", myStdev, myParams->data.F32[1]); + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: the fitted Gaussian has more than 10% error for the mean.\n"); + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: Using the calculated mean (calc, fit) is (%f, %f)\n", + myMean, myParams->data.F32[0]); stats->robustStdev = myStdev; } else { @@ -1691,11 +1723,18 @@ psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED); + // XXX: Set to NAN } } if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { - if (0 != p_psVectorClippedStats(inF32, mask, maskVal, stats)) { + if (-1 == p_psVectorClippedStats(inF32, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, - "Failed to calculate statistics for input psVector."); + "Failed to calculate clipped statistics for input psVector."); + stats->clippedMean = NAN; + stats->clippedStdev = NAN; + } else if (-2 == p_psVectorClippedStats(inF32, mask, maskVal, stats)) { + psLogMsg(__func__, PS_LOG_WARN, "Failed to calculate clipped statistics for input psVector."); + stats->clippedMean = NAN; + stats->clippedStdev = NAN; } }