Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 2329) +++ trunk/psLib/src/math/psStats.c (revision 2339) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.89 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-10 23:22:32 $ + * @version $Revision: 1.90 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-11-11 20:13:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -53,5 +53,4 @@ /*****************************************************************************/ psVector* p_psConvertToF32(psVector* in); - /*****************************************************************************/ /* GLOBAL VARIABLES */ @@ -124,15 +123,20 @@ /****************************************************************************** - MISC PRIVATE STATISTICAL FUNCTIONS +MISC PRIVATE STATISTICAL FUNCTIONS - 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(). - Is the mask vector the same length as the data vector? - Is the mask vector of type PS_TYPE_U8? - Is the stats data structure NULL? - Is the in data structure NULL? - Is the in data structure of type PS_TYPE_F32? +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? + +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 +range, or there are no unmasked elements, or the input vector is NULL)? +Currently we set the value to NAN. + +XXX: Optimization: many routines have an "empty" boolean variable which keeps +track of whether or not the vector has any valid elements. This code can +possibly be optimized away. *****************************************************************************/ /****************************************************************************** @@ -146,6 +150,6 @@ Returns NULL - *****************************************************************************/ - + + *****************************************************************************/ void p_psVectorSampleMean(const psVector* restrict myVector, const psVector* restrict maskVector, @@ -156,18 +160,14 @@ float mean = 0.0; // The mean psS32 count = 0; // # of points in this mean - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this // If PS_STAT_USE_RANGE is requested, then we enter a slightly different // loop. if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; if (maskVector != NULL) { for (i = 0; i < myVector->n; i++) { // Check if the data is with the specified range 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)) { mean += myVector->data.F32[i]; count++; @@ -182,5 +182,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)) { mean += myVector->data.F32[i]; count++; @@ -234,18 +235,17 @@ { psS32 i = 0; // Loop index variable - float max = -PS_MAX_F32; // The calculated maximum - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this + float max = -PS_MAX_F32; // The calculated maximum + psS32 empty = true; // Does this vector have valid elements? // If PS_STAT_USE_RANGE is requested, then we enter a different loop. if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; if (maskVector != NULL) { for (i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i])) { if ((myVector->data.F32[i] > max) && - (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { + (stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { max = myVector->data.F32[i]; + empty = false; } } @@ -254,6 +254,8 @@ for (i = 0; i < myVector->n; i++) { if ((myVector->data.F32[i] > max) && - (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { + (stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { max = myVector->data.F32[i]; + empty = false; } } @@ -265,4 +267,5 @@ if (myVector->data.F32[i] > max) { max = myVector->data.F32[i]; + empty = false; } } @@ -272,10 +275,14 @@ if (myVector->data.F32[i] > max) { max = myVector->data.F32[i]; - } - } - } - } - - stats->max = max; + empty = false; + } + } + } + } + if (empty == false) { + stats->max = max; + } else { + stats->max = NAN; + } } @@ -298,16 +305,15 @@ psS32 i = 0; // Loop index variable float min = PS_MAX_F32; // The calculated maximum - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this + psS32 empty = true; // Does this vector have valid elements? if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; if (maskVector != NULL) { for (i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i])) { if ((myVector->data.F32[i] < min) && - (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { + (stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { min = myVector->data.F32[i]; + empty = false; } } @@ -316,6 +322,8 @@ for (i = 0; i < myVector->n; i++) { if ((myVector->data.F32[i] < min) && - (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { + (stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { min = myVector->data.F32[i]; + empty = false; } } @@ -327,4 +335,5 @@ if (myVector->data.F32[i] < min) { min = myVector->data.F32[i]; + empty = false; } } @@ -334,10 +343,15 @@ if (myVector->data.F32[i] < min) { min = myVector->data.F32[i]; - } - } - } - } - - stats->min = min; + empty = false; + } + } + } + } + + if (empty == false) { + stats->min = min; + } else { + stats->min = NAN; + } } @@ -345,5 +359,5 @@ p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the number of non-masked pixels in the vector that fall within the min/max -range, if given. +range, if specified. Inputs myVector @@ -363,14 +377,11 @@ psS32 i = 0; // Loop index variable psS32 numData = 0; // The number of data points - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; 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)) { numData++; } @@ -378,5 +389,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)) { numData++; } @@ -394,4 +406,5 @@ } } + return (numData); } @@ -408,5 +421,5 @@ NULL *****************************************************************************/ -void p_psVectorSampleMedian(const psVector* restrict myVector, +bool p_psVectorSampleMedian(const psVector* restrict myVector, const psVector* restrict maskVector, psU32 maskVal, @@ -418,6 +431,4 @@ psS32 count = 0; // # of points in this mean? psS32 nValues = 0; // # of points in vector - 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 @@ -425,13 +436,17 @@ nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); + // XXX: Is a warning message appropriate? What value should we set the + // sampleMedian to? Should we generate an error? + if (nValues <= 0) { + printf("WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n"); + return(true); + stats->sampleMedian = 0; + } + // 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. @@ -440,12 +455,16 @@ for (i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = maskVector->data.F32[i]; + (stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { + unsortedVector->data.F32[count] = myVector->data.F32[i]; + count++; } } } else { for (i = 0; i < myVector->n; i++) { - if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; + if ((stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { + unsortedVector->data.F32[count] = myVector->data.F32[i]; + count++; } } @@ -457,5 +476,6 @@ for (i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; + unsortedVector->data.F32[count] = myVector->data.F32[i]; + count++; } } @@ -467,5 +487,11 @@ } // 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 median exactly. @@ -480,4 +506,7 @@ psFree(unsortedVector); psFree(sortedVector); + + // Return "true" on success. + return(true); } @@ -485,8 +514,16 @@ p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input robustHistogram with a Gaussian of width sigma. + +XXX: Only PS_TYPE_F32 is supported. + +XXX: Write a general routine which smoothes a psVector. This routine should +call that. Is that possible? *****************************************************************************/ psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, float sigma) { + PS_PTR_CHECK_NULL(robustHistogram, NULL); + PS_PTR_CHECK_NULL(robustHistogram->bounds, NULL); + psS32 i = 0; // Loop index variable psS32 j = 0; // Loop index variable @@ -504,6 +541,12 @@ x.type.type = PS_TYPE_F32; for (i = 0; i < numBins; i++) { + // Determine the midpoint of bin i. iMid = (robustHistogram->bounds->data.F32[i] + robustHistogram->bounds->data.F32[i+1]) / 2.0; + + + // We determine the bin numbers corresponding to a range of data + // values surrounding iMid. The ranges is of size + // s*PS_GAUSS_WIDTH*sigma // YYY: The p_psVectorBinDisect() routine does much of the work of @@ -516,4 +559,10 @@ if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { jMin = p_psVectorBinDisect(robustHistogram->bounds, &x); + if (jMin < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } } else if (x.data.F32 <= firstBound) { jMin = 0;