Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 1341) +++ trunk/psLib/src/math/psStats.c (revision 1361) @@ -32,37 +32,37 @@ #define MAX_ITERATIONS 10 -void p_psVectorRobustStats(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats); - +void p_psVectorRobustStats( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ); + /** Preprocessor macro to generate error on an incorrect type */ #define PS_CHECK_VECTOR_TYPE(NAME, TYPE) \ if (NAME->type.type != TYPE) { \ - psError(__func__,"Invalid operation: %s has incorrect type.", #NAME); \ -} - + psError(__func__,"Invalid operation: %s has incorrect type.", #NAME); \ + } + /** Preprocessor macro to generate error on a NULL vector */ #define PS_CHECK_NULL_VECTOR(NAME) \ if (NAME == NULL || NAME->data.V == NULL) { \ - psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \ -} - + psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \ + } + /** Preprocessor macro to generate error for zero length vector */ #define PS_CHECK_EMPTY_VECTOR(NAME) \ if (NAME->n < 1) { \ - psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \ -} - + psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \ + } + /** Preprocessor macro to generate error on differing size vectors */ #define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2) \ if (VEC1->n != VEC2->n) { \ - psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \ -} - + psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \ + } + #define PS_PRINT_VECTOR(NAME) \ for (int my_i=0;my_in;my_i++) { \ - printf("%s->data.F32[%d] is %f\n", #NAME, my_i, NAME->data.F32[my_i]); \ -} \ + printf("%s->data.F32[%d] is %f\n", #NAME, my_i, NAME->data.F32[my_i]); \ + } \ printf("\n"); \ @@ -87,56 +87,56 @@ /*****************************************************************************/ -bool p_psGetStatValue(const psStats* stats, double* value) -{ - - switch (stats->options & - ~ (PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE)) { - case PS_STAT_SAMPLE_MEAN: - *value = stats->sampleMean; - return true; - - case PS_STAT_SAMPLE_MEDIAN: - *value = stats->sampleMedian; - return true; - - case PS_STAT_SAMPLE_STDEV: - *value = stats->sampleStdev; - return true; - - case PS_STAT_ROBUST_MEAN: - *value = stats->robustMean; - return true; - - case PS_STAT_ROBUST_MEDIAN: - *value = stats->robustMedian; - return true; - - case PS_STAT_ROBUST_MODE: - *value = stats->robustMode; - return true; - - case PS_STAT_ROBUST_STDEV: - *value = stats->robustStdev; - return true; - - case PS_STAT_CLIPPED_MEAN: - *value = stats->clippedMean; - return true; - - case PS_STAT_CLIPPED_STDEV: - *value = stats->clippedStdev; - return true; - - case PS_STAT_MAX: - *value = stats->max; - return true; - - case PS_STAT_MIN: - *value = stats->min; - return true; - - default: - return false; - } +bool p_psGetStatValue( const psStats* stats, double* value ) +{ + + switch ( stats->options & + ~ ( PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE ) ) { + case PS_STAT_SAMPLE_MEAN: + *value = stats->sampleMean; + return true; + + case PS_STAT_SAMPLE_MEDIAN: + *value = stats->sampleMedian; + return true; + + case PS_STAT_SAMPLE_STDEV: + *value = stats->sampleStdev; + return true; + + case PS_STAT_ROBUST_MEAN: + *value = stats->robustMean; + return true; + + case PS_STAT_ROBUST_MEDIAN: + *value = stats->robustMedian; + return true; + + case PS_STAT_ROBUST_MODE: + *value = stats->robustMode; + return true; + + case PS_STAT_ROBUST_STDEV: + *value = stats->robustStdev; + return true; + + case PS_STAT_CLIPPED_MEAN: + *value = stats->clippedMean; + return true; + + case PS_STAT_CLIPPED_STDEV: + *value = stats->clippedStdev; + return true; + + case PS_STAT_MAX: + *value = stats->max; + return true; + + case PS_STAT_MIN: + *value = stats->min; + return true; + + default: + return false; + } } @@ -147,17 +147,17 @@ *****************************************************************************/ -void p_psVectorPrint(psVector *myVector, - psVector *maskVector, - unsigned int maskVal, - psStats *stats) +void p_psVectorPrint( psVector *myVector, + psVector *maskVector, + unsigned int maskVal, + psStats *stats ) { int i = 0; // Loop index variable. - - for (i=0;in;i++) { - if (maskVector != NULL) - printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]); - else - printf("Element %d is %f\n", i, myVector->data.F32[i]); - } + + for ( i = 0;i < myVector->n;i++ ) { + if ( maskVector != NULL ) + printf( "Element %d is %f (mask is %d)\n", i, myVector->data.F32[ i ], maskVector->data.U8[ i ] ); + else + printf( "Element %d is %f\n", i, myVector->data.F32[ i ] ); + } } @@ -194,8 +194,8 @@ *****************************************************************************/ -void p_psVectorSampleMean(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) +void p_psVectorSampleMean( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ) { int i = 0; // Loop index variable @@ -204,48 +204,48 @@ 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;in;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)) { - mean+= myVector->data.F32[i]; - count++; - } - } - mean/= (float) count; + 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 ) ) { + mean += myVector->data.F32[ i ]; + count++; + } + } + mean /= ( float ) count; + } else { + for ( i = 0;i < myVector->n;i++ ) { + if ( ( rangeMin <= myVector->data.F32[ i ] ) && + ( myVector->data.F32[ i ] <= rangeMax ) ) { + mean += myVector->data.F32[ i ]; + count++; + } + } + mean /= ( float ) count; + } } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - mean+= myVector->data.F32[i]; - count++; - } - } - mean/= (float) count; - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - mean+= myVector->data.F32[i]; - count++; - } - } - mean/= (float) count; - } else { - for (i=0;in;i++) { - mean+= myVector->data.F32[i]; - } - mean/= (float) myVector->n; - } - } - + if ( maskVector != NULL ) { + for ( i = 0;i < myVector->n;i++ ) { + if ( !( maskVal & maskVector->data.U8[ i ] ) ) { + mean += myVector->data.F32[ i ]; + count++; + } + } + mean /= ( float ) count; + } else { + for ( i = 0;i < myVector->n;i++ ) { + mean += myVector->data.F32[ i ]; + } + mean /= ( float ) myVector->n; + } + } + stats->sampleMean = mean; } @@ -262,8 +262,8 @@ NULL *****************************************************************************/ -void p_psVectorMax(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) +void p_psVectorMax( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ) { int i = 0; // Loop index variable @@ -271,46 +271,46 @@ 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 different loop. - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if ((myVector->data.F32[i] > max) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - max = myVector->data.F32[i]; - } - } - } + 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 ) ) { + max = myVector->data.F32[ i ]; + } + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + if ( ( myVector->data.F32[ i ] > max ) && + ( rangeMin <= myVector->data.F32[ i ] ) && + ( myVector->data.F32[ i ] <= rangeMax ) ) { + max = myVector->data.F32[ i ]; + } + } + } } else { - for (i=0;in;i++) { - if ((myVector->data.F32[i] > max) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - max = myVector->data.F32[i]; - } - } - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if (myVector->data.F32[i] > max) { - max = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if (myVector->data.F32[i] > max) { - max = myVector->data.F32[i]; - } - } - } - } - + if ( maskVector != NULL ) { + for ( i = 0;i < myVector->n;i++ ) { + if ( !( maskVal & maskVector->data.U8[ i ] ) ) { + if ( myVector->data.F32[ i ] > max ) { + max = myVector->data.F32[ i ]; + } + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + if ( myVector->data.F32[ i ] > max ) { + max = myVector->data.F32[ i ]; + } + } + } + } + stats->max = max; } @@ -327,8 +327,8 @@ NULL *****************************************************************************/ -void p_psVectorMin(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) +void p_psVectorMin( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ) { int i = 0; // Loop index variable @@ -336,45 +336,45 @@ 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;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if ((myVector->data.F32[i] < min) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - min = myVector->data.F32[i]; - } - } - } + + 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 ) ) { + min = myVector->data.F32[ i ]; + } + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + if ( ( myVector->data.F32[ i ] < min ) && + ( rangeMin <= myVector->data.F32[ i ] ) && + ( myVector->data.F32[ i ] <= rangeMax ) ) { + min = myVector->data.F32[ i ]; + } + } + } } else { - for (i=0;in;i++) { - if ((myVector->data.F32[i] < min) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - min = myVector->data.F32[i]; - } - } - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if (myVector->data.F32[i] < min) { - min = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if (myVector->data.F32[i] < min) { - min = myVector->data.F32[i]; - } - } - } - } - + if ( maskVector != NULL ) { + for ( i = 0;i < myVector->n;i++ ) { + if ( !( maskVal & maskVector->data.U8[ i ] ) ) { + if ( myVector->data.F32[ i ] < min ) { + min = myVector->data.F32[ i ]; + } + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + if ( myVector->data.F32[ i ] < min ) { + min = myVector->data.F32[ i ]; + } + } + } + } + stats->min = min; } @@ -392,8 +392,8 @@ NULL *****************************************************************************/ -int p_psVectorNValues(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) +int p_psVectorNValues( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ) { int i = 0; // Loop index variable @@ -401,38 +401,38 @@ 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;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - numData++; - } - } + + 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 ) ) { + numData++; + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + if ( ( rangeMin <= myVector->data.F32[ i ] ) && + ( myVector->data.F32[ i ] <= rangeMax ) ) { + numData++; + } + } + } } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - numData++; - } - } - } - } else { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - numData++; - } - } - } else { - numData = myVector->n; - } - } - return(numData); + rangeMin = stats->min; + rangeMax = stats->max; + if ( maskVector != NULL ) { + for ( i = 0;i < myVector->n;i++ ) { + if ( !( maskVal & maskVector->data.U8[ i ] ) ) { + numData++; + } + } + } else { + numData = myVector->n; + } + } + return ( numData ); } @@ -450,10 +450,10 @@ NULL *****************************************************************************/ -void p_psVectorSampleMedian(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - psVector *unsortedVector = NULL; // Temporary vector +void p_psVectorSampleMedian( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ) +{ + psVector * unsortedVector = NULL; // Temporary vector psVector *sortedVector = NULL; // Temporary vector int i = 0; // Loop index variable @@ -462,6 +462,6 @@ float rangeMin = 0.0; // Exclude data below this float rangeMax = 0.0; // Exclude date above this - - + + // Determine if the number of data points exceed a threshold which will // cause to generate robust stats, as opposed to exact stats. @@ -488,70 +488,70 @@ } */ - + // Determine how many data points fit inside this min/max range // and are not masked, IF the maskVector is not NULL> - nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); - + nValues = p_psVectorNValues( myVector, maskVector, maskVal, stats ); + // Allocate temporary vectors for the data. - unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + unsortedVector = psVectorAlloc( nValues, PS_TYPE_F32 ); unsortedVector->n = unsortedVector->nalloc; - - sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + + sortedVector = psVectorAlloc( nValues, PS_TYPE_F32 ); 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->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = maskVector->data.F32[i]; - } - } + 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;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 ]; + } + } + } 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 ]; + } + } + } } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } - } else { - // Store all non-masked data points into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - unsortedVector->data.F32[i] = myVector->data.F32[i]; - } - } - } + // Store all non-masked data points into the temporary vectors. + count = 0; + if ( maskVector != NULL ) { + for ( i = 0;i < myVector->n;i++ ) { + if ( !( maskVal & maskVector->data.U8[ i ] ) ) { + unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ]; + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + unsortedVector->data.F32[ i ] = myVector->data.F32[ i ]; + } + } + } // Sort the temporary vectors. - psVectorSort(sortedVector, unsortedVector); - + psVectorSort( sortedVector, unsortedVector ); + // Calculate the median exactly. // XXX: Is this the correct action? - if (0 == (nValues % 2)) { - stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] + - sortedVector->data.F32[nValues/2]); - } else { - stats->sampleMedian = sortedVector->data.F32[nValues/2]; - } - + if ( 0 == ( nValues % 2 ) ) { + stats->sampleMedian = 0.5 * ( sortedVector->data.F32[ ( nValues / 2 ) - 1 ] + + sortedVector->data.F32[ nValues / 2 ] ); + } else { + stats->sampleMedian = sortedVector->data.F32[ nValues / 2 ]; + } + // Free the temporary data structures. - psFree(unsortedVector); - psFree(sortedVector); + psFree( unsortedVector ); + psFree( sortedVector ); } @@ -568,6 +568,6 @@ XXX: use a static variable for gaussianCoefs[] and compute them once. *****************************************************************************/ -psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram, - float sigma) +psVector *p_psVectorsmoothHistGaussian( psHistogram *robustHistogram, + float sigma ) { int i = 0; // Loop index variable @@ -575,50 +575,50 @@ float denom = 0.0; // Temporary variable float expo = 0.0; // Temporary variable - float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients - psVector *smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32); - - for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) { - if (fabs(sigma) >= FLT_EPSILON) { - // If sigma does not equal zero, then we use Gaussian smoothing. - #ifdef DARWIN - denom = (float) sqrt(2.0 * M_PI * sigma * sigma); - #else - - denom = sqrtf(2.0 * M_PI * sigma * sigma); - #endif - - expo = - (float) ((i-GAUSS_WIDTH) * (i-GAUSS_WIDTH)); - expo/= (2.0 * sigma * sigma); - gaussianCoefs[i] = exp(expo/denom); - - // NOTE: Gaussian smoothing just isn't working with low sigma - // values. The problem is that the Gaussian coefficients are - // all zero, except for the middle coefficient, which is exactly - // one. Therefore, I'm using boxcar smoothing. - gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float) GAUSS_WIDTH)); - // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); - } else { - /* If sigma equals zero (all pixels have the same value) - * the above code will divide by zero. Therefore, we don't need - * to smooth the data. - */ - for(i=0;inums->n;i++) { - smooth->data.F32[i] = (float) robustHistogram->nums->data.S32[i]; - } - return(smooth); - } - } - + float gaussianCoefs[ 1 + ( 2 * GAUSS_WIDTH ) ]; // The Gaussian Coefficients + psVector *smooth = psVectorAlloc( robustHistogram->nums->n, PS_TYPE_F32 ); + + for ( i = 0;i < ( 1 + ( 2 * GAUSS_WIDTH ) );i++ ) { + if ( fabs( sigma ) >= FLT_EPSILON ) { + // If sigma does not equal zero, then we use Gaussian smoothing. + #ifdef DARWIN + denom = ( float ) sqrt( 2.0 * M_PI * sigma * sigma ); + #else + + denom = sqrtf( 2.0 * M_PI * sigma * sigma ); + #endif + + expo = - ( float ) ( ( i - GAUSS_WIDTH ) * ( i - GAUSS_WIDTH ) ); + expo /= ( 2.0 * sigma * sigma ); + gaussianCoefs[ i ] = exp( expo / denom ); + + // NOTE: Gaussian smoothing just isn't working with low sigma + // values. The problem is that the Gaussian coefficients are + // all zero, except for the middle coefficient, which is exactly + // one. Therefore, I'm using boxcar smoothing. + gaussianCoefs[ i ] = 1.0 / ( 1.0 + ( 2.0 * ( float ) GAUSS_WIDTH ) ); + // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); + } else { + /* If sigma equals zero (all pixels have the same value) + * the above code will divide by zero. Therefore, we don't need + * to smooth the data. + */ + for ( i = 0;i < robustHistogram->nums->n;i++ ) { + smooth->data.F32[ i ] = ( float ) robustHistogram->nums->data.S32[ i ]; + } + return ( smooth ); + } + } + // Perform the actual smoothing. - for(i=0;inums->n;i++) { - smooth->data.F32[i] = 0.0; - for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) { - if (((j+i) >= 0) && ((j+i) < smooth->n)) { - smooth->data.F32[i]+= (gaussianCoefs[j+GAUSS_WIDTH] * - (float) robustHistogram->nums->data.S32[j+i]); - } - } - } - return(smooth); + for ( i = 0;i < robustHistogram->nums->n;i++ ) { + smooth->data.F32[ i ] = 0.0; + for ( j = -GAUSS_WIDTH;j <= + GAUSS_WIDTH;j++ ) { + if ( ( ( j + i ) >= 0 ) && ( ( j + i ) < smooth->n ) ) { + smooth->data.F32[ i ] += ( gaussianCoefs[ j + GAUSS_WIDTH ] * + ( float ) robustHistogram->nums->data.S32[ j + i ] ); + } + } + } + return ( smooth ); } @@ -634,10 +634,10 @@ NULL *****************************************************************************/ -void p_psVectorSampleQuartiles(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - psVector *unsortedVector = NULL; // Temporary vector +void p_psVectorSampleQuartiles( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ) +{ + psVector * unsortedVector = NULL; // Temporary vector psVector *sortedVector = NULL; // Temporary vector int i = 0; // Loop index variable @@ -646,64 +646,64 @@ 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> - nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); - + nValues = p_psVectorNValues( myVector, maskVector, maskVal, stats ); + // Allocate temporary vectors for the data. - unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + unsortedVector = psVectorAlloc( nValues, PS_TYPE_F32 ); unsortedVector->n = unsortedVector->nalloc; - sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + sortedVector = psVectorAlloc( nValues, PS_TYPE_F32 ); 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->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } + 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;i < myVector->n;i++ ) { + if ( !( maskVal & maskVector->data.U8[ i ] ) && + ( rangeMin <= myVector->data.F32[ i ] ) && + ( myVector->data.F32[ i ] <= rangeMax ) ) { + unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ]; + } + } + } 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 ]; + } + } + } } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } - } else { - // Store all non-masked data points into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - unsortedVector->data.F32[i] = myVector->data.F32[i]; - } - } - } - + // Store all non-masked data points into the temporary vectors. + count = 0; + if ( maskVector != NULL ) { + for ( i = 0;i < myVector->n;i++ ) { + if ( !( maskVal & maskVector->data.U8[ i ] ) ) { + unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ]; + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + unsortedVector->data.F32[ i ] = myVector->data.F32[ i ]; + } + } + } + // Sort the temporary vectors. - psVectorSort(sortedVector, unsortedVector); - + psVectorSort( sortedVector, unsortedVector ); + // Calculate the quartile points exactly. - stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)]; - stats->sampleLQ = sortedVector->data.F32[nValues / 4]; - + stats->sampleUQ = sortedVector->data.F32[ 3 * ( nValues / 4 ) ]; + stats->sampleLQ = sortedVector->data.F32[ nValues / 4 ]; + // Free the temporary data structures. - psFree(unsortedVector); - psFree(sortedVector); + psFree( unsortedVector ); + psFree( sortedVector ); } @@ -721,8 +721,8 @@ *****************************************************************************/ -void p_psVectorSampleStdev(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) +void p_psVectorSampleStdev( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ) { int i = 0; // Loop index variable @@ -735,66 +735,66 @@ 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 // calculated, then call p_psVectorSampleMean() - if (0 != isnan(stats->sampleMean)) { - p_psVectorSampleMean(myVector, maskVector, maskVal, stats); - } + if ( 0 != isnan( stats->sampleMean ) ) { + p_psVectorSampleMean( myVector, maskVector, maskVal, stats ); + } mean = stats->sampleMean; - - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } + + 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 ) ) { + diff = myVector->data.F32[ i ] - mean; + sumSquares += ( diff * diff ); + sumDiffs += diff; + countInt++; + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + if ( ( rangeMin <= myVector->data.F32[ i ] ) && + ( myVector->data.F32[ i ] <= rangeMax ) ) { + diff = myVector->data.F32[ i ] - mean; + sumSquares += ( diff * diff ); + sumDiffs += diff; + countInt++; + } + } + countInt = myVector->n; + } } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - countInt = myVector->n; - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - } else { - for (i=0;in;i++) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - countInt = myVector->n; - } - } - countFloat = (float) countInt; - + if ( maskVector != NULL ) { + for ( i = 0;i < myVector->n;i++ ) { + if ( !( maskVal & maskVector->data.U8[ i ] ) ) { + diff = myVector->data.F32[ i ] - mean; + sumSquares += ( diff * diff ); + sumDiffs += diff; + countInt++; + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + diff = myVector->data.F32[ i ] - mean; + sumSquares += ( diff * diff ); + sumDiffs += diff; + countInt++; + } + countInt = myVector->n; + } + } + 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 } @@ -812,8 +812,8 @@ NULL *****************************************************************************/ -void p_psVectorClippedStats(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) +void p_psVectorClippedStats( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ) { int i = 0; // Loop index variable @@ -824,86 +824,86 @@ float oldStanStdev = 0.0; // Temporary variable psVector *tmpMask = NULL; // Temporary vector - + // Endure that stats->clipIter is within the proper range. - if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) && - (stats->clipIter <= CLIPPED_NUM_ITER_UB))) { - psAbort(__func__, "Unallowed value for clipIter (%d).\n", - stats->clipIter); - } - + if ( !( ( CLIPPED_NUM_ITER_LB <= stats->clipIter ) && + ( stats->clipIter <= CLIPPED_NUM_ITER_UB ) ) ) { + psAbort( __func__, "Unallowed value for clipIter (%d).\n", + stats->clipIter ); + } + // Endure that stats->clipSigma is within the proper range. - if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && - (stats->clipSigma <= CLIPPED_SIGMA_UB))) { - psAbort(__func__, "Unallowed value for clipSigma (%f).\n", - stats->clipSigma); - } - + if ( !( ( CLIPPED_SIGMA_LB <= stats->clipSigma ) && + ( stats->clipSigma <= CLIPPED_SIGMA_UB ) ) ) { + psAbort( __func__, "Unallowed value for clipSigma (%f).\n", + stats->clipSigma ); + } + // We allocate a temporary mask vector since during the iterative // steps that follow, we will be masking off additional data points. // However, we do no want to modify the original mask vector. - tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8); + tmpMask = psVectorAlloc( myVector->n, PS_TYPE_U8 ); tmpMask->n = myVector->n; - + // If we were called with a mask vector, then initialize the temporary // mask vector with those values. - if (maskVector != NULL) { - for (i=0;in;i++) { - tmpMask->data.U8[i] = maskVector->data.U8[i]; - } - } - + if ( maskVector != NULL ) { + for ( i = 0;i < tmpMask->n;i++ ) { + tmpMask->data.U8[ i ] = maskVector->data.U8[ i ]; + } + } + // 1. Compute the sample median. // NOTE: This seems odd. Verify with IfA that we want to calculate the // median here, not the mean. - p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); - + p_psVectorSampleMedian( myVector, maskVector, maskVal, stats ); + // 2. Compute the sample standard deviation. - p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); - + p_psVectorSampleStdev( myVector, maskVector, maskVal, stats ); + // 3. Use the sample median as the first estimator of the mean X. clippedMean = stats->sampleMean; - + // 4. Use the sample stdev as the first estimator of the mean stdev. clippedStdev = stats->sampleStdev; - + // Must save the old sampleMean and sampleStdev since the following code // block overwrites them. oldStanMean = stats->sampleMean; oldStanStdev = stats->sampleStdev; - + // 5. Repeat N times: - 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->data.F32[j] - clippedMean) > - (stats->clipSigma * clippedStdev)) { - 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 - clippedMean = stats->sampleMean; - - // d) Use the new stdev for stdev - clippedStdev = stats->sampleStdev; - } - - } + for ( i = 0;i < stats->clipIter;i++ ) { + for ( j = 0;j < myVector->n;j++ ) { + // a) Exclude all values x_i for which |x_i - x| > K * stdev + if ( fabs( myVector->data.F32[ j ] - clippedMean ) > + ( stats->clipSigma * clippedStdev ) ) { + 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 + clippedMean = stats->sampleMean; + + // d) Use the new stdev for stdev + clippedStdev = stats->sampleStdev; + } + + } stats->sampleMean = oldStanMean; - stats->sampleStdev= oldStanStdev; - + stats->sampleStdev = oldStanStdev; + // 7. The last calcuated value of x is the cliped mean. - if (stats->options & PS_STAT_CLIPPED_MEAN) { - stats->clippedMean = clippedMean; - } - + if ( stats->options & PS_STAT_CLIPPED_MEAN ) { + stats->clippedMean = clippedMean; + } + // 8. The last calcuated value of stdev is the cliped stdev. - if (stats->options & PS_STAT_CLIPPED_STDEV) { - stats->clippedStdev = clippedStdev; - } - - psFree(tmpMask); + if ( stats->options & PS_STAT_CLIPPED_STDEV ) { + stats->clippedStdev = clippedStdev; + } + + psFree( tmpMask ); } @@ -912,24 +912,24 @@ elements of a vector to a range between 0.0 and 1.0. *****************************************************************************/ -void p_psNormalizeVector(psVector *myData) -{ - float min = (float) HUGE; - float max = (float) -HUGE; +void p_psNormalizeVector( psVector *myData ) +{ + float min = ( float ) HUGE; + float max = ( float ) - HUGE; float range = 0.0; int i = 0; - - for (i=0;in;i++) { - if (myData->data.F32[i] < min) { - min = myData->data.F32[i]; - } - if (myData->data.F32[i] > max) { - max = myData->data.F32[i]; - } - } - + + 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 ]; + } + } + range = max - min; - for (i=0;in;i++) { - myData->data.F32[i] = (myData->data.F32[i] - min) / range; - } + for ( i = 0;i < myData->n;i++ ) { + myData->data.F32[ i ] = ( myData->data.F32[ i ] - min ) / range; + } } @@ -940,15 +940,15 @@ specified data point. *****************************************************************************/ -float p_psGaussian(const psVector *restrict myData, - const psVector *restrict myParams) -{ - float x = myData->data.F32[0]; - float mean = myParams->data.F32[0]; - float stdev = myParams->data.F32[1]; - float tmp = exp(-((x-mean) * (x-mean)) / (2.0 * stdev * stdev)); - tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev))); - +float p_psGaussian( const psVector *restrict myData, + const psVector *restrict myParams ) +{ + float x = myData->data.F32[ 0 ]; + float mean = myParams->data.F32[ 0 ]; + float stdev = myParams->data.F32[ 1 ]; + float tmp = exp( -( ( x - mean ) * ( x - mean ) ) / ( 2.0 * stdev * stdev ) ); + tmp /= ( ( float ) sqrt( 2.0 * M_PI * ( stdev * stdev ) ) ); + // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); - return(tmp); + return ( tmp ); } @@ -957,25 +957,25 @@ calculates the specified partial derivative of the above Gaussian function. *****************************************************************************/ -float p_psGaussianDeriv(const psVector *restrict myData, - const psVector *restrict myParams, - int whichParam) -{ - float x = myData->data.F32[0]; - float mean = myParams->data.F32[0]; - float stdev = myParams->data.F32[1]; +float p_psGaussianDeriv( const psVector *restrict myData, + const psVector *restrict myParams, + int whichParam ) +{ + float x = myData->data.F32[ 0 ]; + float mean = myParams->data.F32[ 0 ]; + float stdev = myParams->data.F32[ 1 ]; float tmp = 0.0; - - if (whichParam == 0) { - // Return the derivative w.r.t. the mean. - tmp = (x - mean) * p_psGaussian(myData, myParams); - tmp/= (stdev * stdev); - } else if (whichParam == 1) { - // Return the derivative w.r.t. the stdev. - tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams); - tmp/= (stdev * stdev * stdev); - } - printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp); - - return(tmp); + + if ( whichParam == 0 ) { + // Return the derivative w.r.t. the mean. + tmp = ( x - mean ) * p_psGaussian( myData, myParams ); + tmp /= ( stdev * stdev ); + } else if ( whichParam == 1 ) { + // Return the derivative w.r.t. the stdev. + tmp = ( x - mean ) * ( x - mean ) * p_psGaussian( myData, myParams ); + tmp /= ( stdev * stdev * stdev ); + } + printf( "p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp ); + + return ( tmp ); } @@ -986,15 +986,15 @@ specified data point. *****************************************************************************/ -float p_psQuadratic(const psVector *restrict myParams, - const psVector *restrict myCoords) -{ - float x = myCoords->data.F32[0]; - float A = myParams->data.F32[0]; - float B = myParams->data.F32[1]; - float C = myParams->data.F32[2]; +float p_psQuadratic( const psVector *restrict myParams, + const psVector *restrict myCoords ) +{ + float x = myCoords->data.F32[ 0 ]; + float A = myParams->data.F32[ 0 ]; + float B = myParams->data.F32[ 1 ]; + float C = myParams->data.F32[ 2 ]; float tmp = 0.0; - - tmp = (A * x * x) + (B * x) + C; - return(tmp); + + tmp = ( A * x * x ) + ( B * x ) + C; + return ( tmp ); } @@ -1003,20 +1003,20 @@ calculates the specified partial derivative of the above quadratic function. *****************************************************************************/ -float p_psQuadraticDeriv(const psVector *restrict myParams, - const psVector *restrict myCoords, - int whichParamDeriv) -{ - float x = myCoords->data.F32[0]; +float p_psQuadraticDeriv( const psVector *restrict myParams, + const psVector *restrict myCoords, + int whichParamDeriv ) +{ + float x = myCoords->data.F32[ 0 ]; float tmp = 0.0; - - if (whichParamDeriv == 0) { - tmp = x * x; - } else if (whichParamDeriv == 1) { - tmp = x; - } else if (whichParamDeriv == 2) { - tmp = 1.0; - } - - return(tmp); + + if ( whichParamDeriv == 0 ) { + tmp = x * x; + } else if ( whichParamDeriv == 1 ) { + tmp = x; + } else if ( whichParamDeriv == 2 ) { + tmp = 1.0; + } + + return ( tmp ); } @@ -1030,37 +1030,37 @@ decreasing within that range. *****************************************************************************/ -float p_ps1DPolyMedian(psPolynomial1D *myPoly, - float rangeLow, - float rangeHigh, - float getThisValue) -{ - int numIterations=0; +float p_ps1DPolyMedian( psPolynomial1D *myPoly, + float rangeLow, + float rangeHigh, + float getThisValue ) +{ + int numIterations = 0; float midpoint = 0.0; float oldMidpoint = 1.0; float f = 0.0; - + // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue); - - while (numIterations < MAX_ITERATIONS) { - midpoint = (rangeHigh + rangeLow) / 2.0; - if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) { - return(midpoint); - } - oldMidpoint = midpoint; - - f = psPolynomial1DEval(midpoint, myPoly); - // printf("p_ps1DPolyMedian() iteration %d. f(%f) is %f\n", numIterations, midpoint, f); - if (fabs(f - getThisValue) <= FLT_EPSILON) { - return(midpoint); - } - - if (f > getThisValue) { - rangeHigh = midpoint; - } else { - rangeLow = midpoint; - } - numIterations++; - } - return(midpoint); + + while ( numIterations < MAX_ITERATIONS ) { + midpoint = ( rangeHigh + rangeLow ) / 2.0; + if ( fabs( midpoint - oldMidpoint ) <= FLT_EPSILON ) { + return ( midpoint ); + } + oldMidpoint = midpoint; + + f = psPolynomial1DEval( midpoint, myPoly ); + // printf("p_ps1DPolyMedian() iteration %d. f(%f) is %f\n", numIterations, midpoint, f); + if ( fabs( f - getThisValue ) <= FLT_EPSILON ) { + return ( midpoint ); + } + + if ( f > getThisValue ) { + rangeHigh = midpoint; + } else { + rangeLow = midpoint; + } + numIterations++; + } + return ( midpoint ); } @@ -1073,52 +1073,52 @@ XXX: This function is currently not being used. *****************************************************************************/ -float p_psFitQuadratic(psHistogram *histogram, - psVector *cumulativeSums, - int binNum, - float fitFloat) -{ - psVector *x = psVectorAlloc(3, PS_TYPE_F64); - psVector *y = psVectorAlloc(3, PS_TYPE_F64); - psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); - psPolynomial1D *myPoly = psPolynomial1DAlloc(2); - - if ((binNum > 0) && - (binNum < (histogram->nums->n+1))) { - x->data.F64[0] = (double) 0.5 * - (histogram->bounds->data.F32[binNum-1] + - histogram->bounds->data.F32[binNum]); - x->data.F64[1] = (double) 0.5 * - (histogram->bounds->data.F32[binNum] + - histogram->bounds->data.F32[binNum+1]); - x->data.F64[2] = (double) 0.5 * - (histogram->bounds->data.F32[binNum+1] + - histogram->bounds->data.F32[binNum+2]); - - y->data.F64[0] = cumulativeSums->data.F32[binNum-1]; - y->data.F64[1] = cumulativeSums->data.F32[binNum]; - y->data.F64[2] = cumulativeSums->data.F32[binNum+1]; - - if (!((y->data.F64[0] <= fitFloat) && - (fitFloat <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); - return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], - fitFloat)); - } else { - return(0.5 * (histogram->bounds->data.F32[binNum+1] + - histogram->bounds->data.F32[binNum])); - } - - psFree(x); - psFree(y); - psFree(yErr); - psFree(myPoly); - return(0.0); +float p_psFitQuadratic( psHistogram *histogram, + psVector *cumulativeSums, + int binNum, + float fitFloat ) +{ + psVector * x = psVectorAlloc( 3, PS_TYPE_F64 ); + psVector *y = psVectorAlloc( 3, PS_TYPE_F64 ); + psVector *yErr = psVectorAlloc( 3, PS_TYPE_F64 ); + psPolynomial1D *myPoly = psPolynomial1DAlloc( 2 ); + + if ( ( binNum > 0 ) && + ( binNum < ( histogram->nums->n + 1 ) ) ) { + x->data.F64[ 0 ] = ( double ) 0.5 * + ( histogram->bounds->data.F32[ binNum - 1 ] + + histogram->bounds->data.F32[ binNum ] ); + x->data.F64[ 1 ] = ( double ) 0.5 * + ( histogram->bounds->data.F32[ binNum ] + + histogram->bounds->data.F32[ binNum + 1 ] ); + x->data.F64[ 2 ] = ( double ) 0.5 * + ( histogram->bounds->data.F32[ binNum + 1 ] + + histogram->bounds->data.F32[ binNum + 2 ] ); + + y->data.F64[ 0 ] = cumulativeSums->data.F32[ binNum - 1 ]; + y->data.F64[ 1 ] = cumulativeSums->data.F32[ binNum ]; + y->data.F64[ 2 ] = cumulativeSums->data.F32[ binNum + 1 ]; + + if ( !( ( y->data.F64[ 0 ] <= fitFloat ) && + ( fitFloat <= y->data.F64[ 2 ] ) ) ) { + psAbort( __func__, "p_psVectorRobustStats(0): midpoint not within y-range\n" ); + } + + yErr->data.F64[ 0 ] = 1.0; + yErr->data.F64[ 1 ] = 1.0; + yErr->data.F64[ 2 ] = 1.0; + + myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr ); + return ( p_ps1DPolyMedian( myPoly, x->data.F64[ 0 ], x->data.F64[ 2 ], + fitFloat ) ); + } else { + return ( 0.5 * ( histogram->bounds->data.F32[ binNum + 1 ] + + histogram->bounds->data.F32[ binNum ] ) ); + } + + psFree( x ); + psFree( y ); + psFree( yErr ); + psFree( myPoly ); + return ( 0.0 ); } @@ -1147,10 +1147,10 @@ NULL *****************************************************************************/ -void p_psVectorRobustStats(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - psHistogram *robustHistogram = NULL; +void p_psVectorRobustStats( const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats ) +{ + psHistogram * robustHistogram = NULL; psVector *robustHistogramVector = NULL; float binSize = 0.0; // Size of the histogram bins @@ -1162,5 +1162,5 @@ float dL = 0.0; int numBins = 0; - psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN); + psStats *tmpStats = psStatsAlloc( PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN ); // psImage *domain; // psVector *errors; @@ -1177,8 +1177,8 @@ float sumSquares = 0.0; float sumDiffs = 0.0; - psVector *x = psVectorAlloc(3, PS_TYPE_F64); - psVector *y = psVectorAlloc(3, PS_TYPE_F64); - psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); - psPolynomial1D *myPoly = psPolynomial1DAlloc(2); + psVector *x = psVectorAlloc( 3, PS_TYPE_F64 ); + psVector *y = psVectorAlloc( 3, PS_TYPE_F64 ); + psVector *yErr = psVectorAlloc( 3, PS_TYPE_F64 ); + psPolynomial1D *myPoly = psPolynomial1DAlloc( 2 ); psVector *cumulativeRobustSumsFullRange = NULL; psVector *cumulativeRobustSumsDlRange = NULL; @@ -1187,186 +1187,186 @@ float sumNfit = 0.0; float cumulativeMedian = 0.0; - + // Compute the initial bin size of the robust histogram. This is done // by computing the clipped standard deviation of the vector, and dividing // that by 10.0; - p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats); + p_psVectorClippedStats( myVector, maskVector, maskVal, tmpStats ); binSize = tmpStats->clippedStdev / 10.0f; - + // If stats->clippedStdev == 0.0, then all data elements have the same // value. Therefore, we can set the appropiate results and return. - if (fabs(binSize) <= FLT_EPSILON) { - if (stats->options & PS_STAT_ROBUST_MEAN) { - stats->robustMean = stats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - stats->robustMedian = stats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = stats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_STDEV) { - stats->robustStdev = 0.0; - } - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - stats->robustUQ = stats->clippedMean; - stats->robustLQ = stats->clippedMean; - } - // XXX: Set these to the number of unmasked data points? - stats->robustNfit = 0.0; - stats->robustN50 = 0.0; - psFree(tmpStats); - return; - } - + if ( fabs( binSize ) <= FLT_EPSILON ) { + if ( stats->options & PS_STAT_ROBUST_MEAN ) { + stats->robustMean = stats->clippedMean; + } + if ( stats->options & PS_STAT_ROBUST_MEDIAN ) { + stats->robustMedian = stats->clippedMean; + } + if ( stats->options & PS_STAT_ROBUST_MODE ) { + stats->robustMode = stats->clippedMean; + } + if ( stats->options & PS_STAT_ROBUST_STDEV ) { + stats->robustStdev = 0.0; + } + if ( stats->options & PS_STAT_ROBUST_QUARTILE ) { + stats->robustUQ = stats->clippedMean; + stats->robustLQ = stats->clippedMean; + } + // XXX: Set these to the number of unmasked data points? + stats->robustNfit = 0.0; + stats->robustN50 = 0.0; + psFree( tmpStats ); + return ; + } + // Determine minimum and maximum values in the data vector. - if (isnan(stats->min)) { - p_psVectorMin(myVector, maskVector, maskVal, stats); - } - if (isnan(stats->max)) { - p_psVectorMax(myVector, maskVector, maskVal, stats); - } - + if ( isnan( stats->min ) ) { + p_psVectorMin( myVector, maskVector, maskVal, stats ); + } + if ( isnan( stats->max ) ) { + p_psVectorMax( myVector, maskVector, maskVal, stats ); + } + // Create the histogram structure. NOTE: we can not specify the bin size // precisely since the argument to psHistogramAlloc() is the number of // bins, not the binSize. Also, if we get here, we know that // binSize != 0.0. - numBins = (int) ((stats->max - stats->min) / binSize); - robustHistogram = psHistogramAlloc(stats->min, - stats->max, - numBins); - + numBins = ( int ) ( ( stats->max - stats->min ) / binSize ); + robustHistogram = psHistogramAlloc( stats->min, + stats->max, + numBins ); + // Populate the histogram array. - psVectorHistogram(robustHistogram, myVector, maskVector, maskVal); - + psVectorHistogram( robustHistogram, myVector, maskVector, maskVal ); + // Smooth the histogram. // XXX: is that the right stdev? - robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram, - tmpStats->clippedStdev/4.0f); - + robustHistogramVector = p_psVectorsmoothHistGaussian( robustHistogram, + tmpStats->clippedStdev / 4.0f ); + // The following was necessary to fit a gaussian to the data, since // gaussian functions produce data between 0.0 and 1.0. // p_psNormalizeVector(robustHistogramVector); - + /************************************************************************** Determine the lower/upper quartiles. - **************************************************************************/ + **************************************************************************/ // We define a vector called "cumulativeRobustSums..." where the value at // index position i is equal to the sum of bins 0:i. This will be used // now and later in determining the lower/upper quartiles. - cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0]; - for (i=1;in;i++) { - cumulativeRobustSumsFullRange->data.F32[i] = - cumulativeRobustSumsFullRange->data.F32[i-1] + - robustHistogramVector->data.F32[i]; - } - sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1]; - + cumulativeRobustSumsFullRange = psVectorAlloc( robustHistogramVector->n, PS_TYPE_F32 ); + cumulativeRobustSumsFullRange->data.F32[ 0 ] = robustHistogramVector->data.F32[ 0 ]; + for ( i = 1;i < robustHistogramVector->n;i++ ) { + cumulativeRobustSumsFullRange->data.F32[ i ] = + cumulativeRobustSumsFullRange->data.F32[ i - 1 ] + + robustHistogramVector->data.F32[ i ]; + } + sumRobust = cumulativeRobustSumsFullRange->data.F32[ robustHistogramVector->n - 1 ]; + // Determine the bin number containing the lower quartile point. LQBinNum = -1; - for (i=0;in;i++) { - if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust/4.0)) { - LQBinNum = i; - break; - } - } - + for ( i = 0;i < cumulativeRobustSumsFullRange->n;i++ ) { + if ( cumulativeRobustSumsFullRange->data.F32[ i ] >= ( sumRobust / 4.0 ) ) { + LQBinNum = i; + break; + } + } + // Determine the bin number containing the upper quartile point. UQBinNum = -1; - for (i=cumulativeRobustSumsFullRange->n-1;i>=0;i--) { - if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0*sumRobust/4.0)) { - UQBinNum = i; - break; - } - } - - if ((LQBinNum == -1) || - (UQBinNum == -1)) { - psAbort(__func__, "Could not determine the robust lower/upper quartiles."); - } + for ( i = cumulativeRobustSumsFullRange->n - 1;i >= 0;i-- ) { + if ( cumulativeRobustSumsFullRange->data.F32[ i ] <= ( 3.0 * sumRobust / 4.0 ) ) { + UQBinNum = i; + break; + } + } + + if ( ( LQBinNum == -1 ) || + ( UQBinNum == -1 ) ) { + psAbort( __func__, "Could not determine the robust lower/upper quartiles." ); + } /************************************************************************** Determine the mode in the range LQ:UQ. - **************************************************************************/ + **************************************************************************/ // Determine the bin with the peak value in the range LQ to UQ. maxBinNum = LQBinNum; - maxBinCount = robustHistogramVector->data.F32[LQBinNum]; - sumN50 = (float) robustHistogram->nums->data.S32[LQBinNum]; - for (i=LQBinNum+1;i<=UQBinNum;i++) { - if (robustHistogramVector->data.F32[i] > maxBinCount) { - maxBinNum = i; - maxBinCount = robustHistogramVector->data.F32[i]; - } - sumN50+= (float) robustHistogram->nums->data.S32[i]; - } - + maxBinCount = robustHistogramVector->data.F32[ LQBinNum ]; + sumN50 = ( float ) robustHistogram->nums->data.S32[ LQBinNum ]; + for ( i = LQBinNum + 1;i <= UQBinNum;i++ ) { + if ( robustHistogramVector->data.F32[ i ] > maxBinCount ) { + maxBinNum = i; + maxBinCount = robustHistogramVector->data.F32[ i ]; + } + sumN50 += ( float ) robustHistogram->nums->data.S32[ i ]; + } + // XXX: is dL defined as the value at the LQ/UQ, or the bin number? - dL = (UQBinNum - LQBinNum) / 4; - - printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum); - + dL = ( UQBinNum - LQBinNum ) / 4; + + printf( "(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum ); + /************************************************************************** Determine the mean/stdev for the bins in the range mode-dL to mode+dL **************************************************************************/ - cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - for (i=0;in;i++) { - cumulativeRobustSumsDlRange->data.F32[i] = 0.0; - } + cumulativeRobustSumsDlRange = psVectorAlloc( robustHistogramVector->n, PS_TYPE_F32 ); + for ( i = 0;i < robustHistogramVector->n;i++ ) { + cumulativeRobustSumsDlRange->data.F32[ i ] = 0.0; + } sumNfit = 0.0; cumulativeMedian = 0.0; - for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) { - if ((0 <= i) && (i < robustHistogramVector->n)) { - cumulativeRobustSumsDlRange->data.F32[i] = - cumulativeRobustSumsDlRange->data.F32[i-1] + - robustHistogramVector->data.F32[i]; - cumulativeMedian+= robustHistogramVector->data.F32[i]; - sumNfit+= (float) robustHistogram->nums->data.S32[i]; - } - } - + for ( i = maxBinNum - dL;i <= maxBinNum + dL;i++ ) { + if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) { + cumulativeRobustSumsDlRange->data.F32[ i ] = + cumulativeRobustSumsDlRange->data.F32[ i - 1 ] + + robustHistogramVector->data.F32[ i ]; + cumulativeMedian += robustHistogramVector->data.F32[ i ]; + sumNfit += ( float ) robustHistogram->nums->data.S32[ i ]; + } + } + // Calculate the mean of the smoothed robust histogram in the range // mode-dL to mode+dL. We use the midpoint of each bin as the mean for // that bin (this is a non-exact approximation). myMean = 0.0; - for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) { - if ((0 <= i) && (i < robustHistogramVector->n)) { - myMean+= (robustHistogramVector->data.F32[i]) * 0.5 * - (robustHistogram->bounds->data.F32[i+1] + - robustHistogram->bounds->data.F32[i]); - countFloat+= robustHistogramVector->data.F32[i]; - } - } - myMean/= countFloat; - + for ( i = maxBinNum - dL;i <= maxBinNum + dL;i++ ) { + if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) { + myMean += ( robustHistogramVector->data.F32[ i ] ) * 0.5 * + ( robustHistogram->bounds->data.F32[ i + 1 ] + + robustHistogram->bounds->data.F32[ i ] ); + countFloat += robustHistogramVector->data.F32[ i ]; + } + } + myMean /= countFloat; + // Calculate the stdev of the smoothed robust histogram in the range // mode-dL to mode+dL. We use the midpoint of each bin as the mean for // that bin. - for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) { - if ((0 <= i) && (i < robustHistogramVector->n)) { - diff = (0.5 * (robustHistogram->bounds->data.F32[i+1] + - robustHistogram->bounds->data.F32[i])) - myMean; - sumSquares+= diff * diff * robustHistogramVector->data.F32[i]; - sumDiffs+= diff * robustHistogramVector->data.F32[i]; - } - } - myStdev = sqrt((sumSquares-(sumDiffs * sumDiffs/countFloat))/ (countFloat-1)); - + for ( i = maxBinNum - dL;i <= maxBinNum + dL;i++ ) { + if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) { + diff = ( 0.5 * ( robustHistogram->bounds->data.F32[ i + 1 ] + + robustHistogram->bounds->data.F32[ i ] ) ) - myMean; + sumSquares += diff * diff * robustHistogramVector->data.F32[ i ]; + sumDiffs += diff * robustHistogramVector->data.F32[ i ]; + } + } + myStdev = sqrt( ( sumSquares - ( sumDiffs * sumDiffs / countFloat ) ) / ( countFloat - 1 ) ); + /************************************************************************** Set the appropriate members in the output stats struct. **************************************************************************/ - if (stats->options & PS_STAT_ROBUST_MEAN) { - stats->robustMean = myMean; - } - - if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = 0.5 * - (robustHistogram->bounds->data.F32[maxBinNum] + - robustHistogram->bounds->data.F32[maxBinNum+1]); - } - - if (stats->options & PS_STAT_ROBUST_STDEV) { - stats->robustStdev = myStdev; - } - + if ( stats->options & PS_STAT_ROBUST_MEAN ) { + stats->robustMean = myMean; + } + + if ( stats->options & PS_STAT_ROBUST_MODE ) { + stats->robustMode = 0.5 * + ( robustHistogram->bounds->data.F32[ maxBinNum ] + + robustHistogram->bounds->data.F32[ maxBinNum + 1 ] ); + } + + if ( stats->options & PS_STAT_ROBUST_STDEV ) { + stats->robustStdev = myStdev; + } + // To determine the median (and later, the lower/upper quartile), we fit // a quadratic to the three bins surrounding the bin containing the median. @@ -1374,50 +1374,50 @@ // the cumulative number of data points in all bins up to, and including, // this bin. We then solve the quadratic for - - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n-1))) { - x->data.F64[0] = (double) 0.5 * - (robustHistogram->bounds->data.F32[maxBinNum-1] + - robustHistogram->bounds->data.F32[maxBinNum]); - x->data.F64[1] = (double) 0.5 * - (robustHistogram->bounds->data.F32[maxBinNum] + - robustHistogram->bounds->data.F32[maxBinNum+1]); - x->data.F64[2] = (double) 0.5 * - (robustHistogram->bounds->data.F32[maxBinNum+1] + - robustHistogram->bounds->data.F32[maxBinNum+2]); - - y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum-1]; - y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum]; - y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum+1]; - - // Ensure that cumulativeMedian/2 is actually within the range of the bins - // we are using. - cumulativeMedian*= 0.5; - if (!((y->data.F64[0] <= cumulativeMedian) && - (cumulativeMedian <= y->data.F64[2]))) { - printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]); - psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n"); - } - // XXX: yErr is not currently used by psGetArrayPolynomial(). We - // may have to set this meaningfully later. - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - // Determine the coefficients of the polynomial. - myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); - // Call p_ps1DPolyMedian(), which does a binary search on the - // polynomial, looking for the value x such that - // f(x) = cumulativeMedian. - stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], - x->data.F64[2], cumulativeMedian); - } else { - // If the mode is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum+1] + - robustHistogram->bounds->data.F32[maxBinNum]); - } - } - + + if ( stats->options & PS_STAT_ROBUST_MEDIAN ) { + if ( ( maxBinNum > 0 ) && ( maxBinNum < ( robustHistogram->nums->n - 1 ) ) ) { + x->data.F64[ 0 ] = ( double ) 0.5 * + ( robustHistogram->bounds->data.F32[ maxBinNum - 1 ] + + robustHistogram->bounds->data.F32[ maxBinNum ] ); + x->data.F64[ 1 ] = ( double ) 0.5 * + ( robustHistogram->bounds->data.F32[ maxBinNum ] + + robustHistogram->bounds->data.F32[ maxBinNum + 1 ] ); + x->data.F64[ 2 ] = ( double ) 0.5 * + ( robustHistogram->bounds->data.F32[ maxBinNum + 1 ] + + robustHistogram->bounds->data.F32[ maxBinNum + 2 ] ); + + y->data.F64[ 0 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum - 1 ]; + y->data.F64[ 1 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum ]; + y->data.F64[ 2 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum + 1 ]; + + // Ensure that cumulativeMedian/2 is actually within the range of the bins + // we are using. + cumulativeMedian *= 0.5; + if ( !( ( y->data.F64[ 0 ] <= cumulativeMedian ) && + ( cumulativeMedian <= y->data.F64[ 2 ] ) ) ) { + printf( "((%f), %f, %f)\n", cumulativeMedian, y->data.F64[ 0 ], y->data.F64[ 2 ] ); + psAbort( __func__, "p_psVectorRobustStats(1): midpoint not within y-range\n" ); + } + // XXX: yErr is not currently used by psVectorFitPolynomial1D(). We + // may have to set this meaningfully later. + yErr->data.F64[ 0 ] = 1.0; + yErr->data.F64[ 1 ] = 1.0; + yErr->data.F64[ 2 ] = 1.0; + + // Determine the coefficients of the polynomial. + myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr ); + // Call p_ps1DPolyMedian(), which does a binary search on the + // polynomial, looking for the value x such that + // f(x) = cumulativeMedian. + stats->robustMedian = p_ps1DPolyMedian( myPoly, x->data.F64[ 0 ], + x->data.F64[ 2 ], cumulativeMedian ); + } else { + // If the mode is the first/last histogram bin, then simply use + // the midpoint of that bin. + stats->robustMedian = 0.5 * ( robustHistogram->bounds->data.F32[ maxBinNum + 1 ] + + robustHistogram->bounds->data.F32[ maxBinNum ] ); + } + } + // The lower/upper quartile calculations are very similar to the median // calculations. We fit a quadratic to the array containing the @@ -1425,91 +1425,91 @@ // f(x) equals the lower/upper quartile exactly. // - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1]; - - if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n-1))) { - x->data.F64[0] = (double) 0.5 * - (robustHistogram->bounds->data.F32[LQBinNum-1] + - robustHistogram->bounds->data.F32[LQBinNum]); - x->data.F64[1] = (double) 0.5 * - (robustHistogram->bounds->data.F32[LQBinNum] + - robustHistogram->bounds->data.F32[LQBinNum+1]); - x->data.F64[2] = (double) 0.5 * - (robustHistogram->bounds->data.F32[LQBinNum+1] + - robustHistogram->bounds->data.F32[LQBinNum+2]); - - y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum-1]; - y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum]; - y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum+1]; - - if (!((y->data.F64[0] <= (countFloat/4.0)) && - ((countFloat/4.0) <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); - stats->robustLQ = p_ps1DPolyMedian(myPoly, - x->data.F64[0], - x->data.F64[2], - countFloat/4.0); - - } else { - // If the LQ is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum+1] + - robustHistogram->bounds->data.F32[LQBinNum]); - } - - if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n-1))) { - x->data.F64[0] = (double) 0.5 * - (robustHistogram->bounds->data.F32[UQBinNum-1] + - robustHistogram->bounds->data.F32[UQBinNum]); - x->data.F64[1] = (double) 0.5 * - (robustHistogram->bounds->data.F32[UQBinNum] + - robustHistogram->bounds->data.F32[UQBinNum+1]); - x->data.F64[2] = (double) 0.5 * - (robustHistogram->bounds->data.F32[UQBinNum+1] + - robustHistogram->bounds->data.F32[UQBinNum+2]); - - y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum-1]; - y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum]; - y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum+1]; - - if (!((y->data.F64[0] <= (3.0 * countFloat/4.0)) && - ((3.0 * countFloat/4.0) <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); - stats->robustUQ = p_ps1DPolyMedian(myPoly, - x->data.F64[0], - x->data.F64[2], - 3.0*countFloat/4.0); - } else { - // If the UQ is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum+1] + - robustHistogram->bounds->data.F32[UQBinNum]); - } - } + if ( stats->options & PS_STAT_ROBUST_QUARTILE ) { + countFloat = cumulativeRobustSumsFullRange->data.F32[ robustHistogramVector->n - 1 ]; + + if ( ( LQBinNum > 0 ) && ( LQBinNum < ( robustHistogram->nums->n - 1 ) ) ) { + x->data.F64[ 0 ] = ( double ) 0.5 * + ( robustHistogram->bounds->data.F32[ LQBinNum - 1 ] + + robustHistogram->bounds->data.F32[ LQBinNum ] ); + x->data.F64[ 1 ] = ( double ) 0.5 * + ( robustHistogram->bounds->data.F32[ LQBinNum ] + + robustHistogram->bounds->data.F32[ LQBinNum + 1 ] ); + x->data.F64[ 2 ] = ( double ) 0.5 * + ( robustHistogram->bounds->data.F32[ LQBinNum + 1 ] + + robustHistogram->bounds->data.F32[ LQBinNum + 2 ] ); + + y->data.F64[ 0 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum - 1 ]; + y->data.F64[ 1 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum ]; + y->data.F64[ 2 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum + 1 ]; + + if ( !( ( y->data.F64[ 0 ] <= ( countFloat / 4.0 ) ) && + ( ( countFloat / 4.0 ) <= y->data.F64[ 2 ] ) ) ) { + psAbort( __func__, "p_psVectorRobustStats(2): midpoint not within y-range\n" ); + } + + yErr->data.F64[ 0 ] = 1.0; + yErr->data.F64[ 1 ] = 1.0; + yErr->data.F64[ 2 ] = 1.0; + + myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr ); + stats->robustLQ = p_ps1DPolyMedian( myPoly, + x->data.F64[ 0 ], + x->data.F64[ 2 ], + countFloat / 4.0 ); + + } else { + // If the LQ is the first/last histogram bin, then simply use + // the midpoint of that bin. + stats->robustLQ = 0.5 * ( robustHistogram->bounds->data.F32[ LQBinNum + 1 ] + + robustHistogram->bounds->data.F32[ LQBinNum ] ); + } + + if ( ( UQBinNum > 0 ) && ( UQBinNum < ( robustHistogram->nums->n - 1 ) ) ) { + x->data.F64[ 0 ] = ( double ) 0.5 * + ( robustHistogram->bounds->data.F32[ UQBinNum - 1 ] + + robustHistogram->bounds->data.F32[ UQBinNum ] ); + x->data.F64[ 1 ] = ( double ) 0.5 * + ( robustHistogram->bounds->data.F32[ UQBinNum ] + + robustHistogram->bounds->data.F32[ UQBinNum + 1 ] ); + x->data.F64[ 2 ] = ( double ) 0.5 * + ( robustHistogram->bounds->data.F32[ UQBinNum + 1 ] + + robustHistogram->bounds->data.F32[ UQBinNum + 2 ] ); + + y->data.F64[ 0 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum - 1 ]; + y->data.F64[ 1 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum ]; + y->data.F64[ 2 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum + 1 ]; + + if ( !( ( y->data.F64[ 0 ] <= ( 3.0 * countFloat / 4.0 ) ) && + ( ( 3.0 * countFloat / 4.0 ) <= y->data.F64[ 2 ] ) ) ) { + psAbort( __func__, "p_psVectorRobustStats(3): midpoint not within y-range\n" ); + } + + yErr->data.F64[ 0 ] = 1.0; + yErr->data.F64[ 1 ] = 1.0; + yErr->data.F64[ 2 ] = 1.0; + + myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr ); + stats->robustUQ = p_ps1DPolyMedian( myPoly, + x->data.F64[ 0 ], + x->data.F64[ 2 ], + 3.0 * countFloat / 4.0 ); + } else { + // If the UQ is the first/last histogram bin, then simply use + // the midpoint of that bin. + stats->robustUQ = 0.5 * ( robustHistogram->bounds->data.F32[ UQBinNum + 1 ] + + robustHistogram->bounds->data.F32[ UQBinNum ] ); + } + } stats->robustNfit = sumNfit; stats->robustN50 = sumN50; - - psFree(x); - psFree(y); - psFree(yErr); - psFree(tmpStats); - psFree(robustHistogram); - psFree(myPoly); - psFree(cumulativeRobustSumsFullRange); - psFree(cumulativeRobustSumsDlRange); + + psFree( x ); + psFree( y ); + psFree( yErr ); + psFree( tmpStats ); + psFree( robustHistogram ); + psFree( myPoly ); + psFree( cumulativeRobustSumsFullRange ); + psFree( cumulativeRobustSumsDlRange ); } @@ -1568,14 +1568,14 @@ /*****************************************************************************/ -static void histogramFree(psHistogram *myHist); +static void histogramFree( psHistogram *myHist ); /****************************************************************************** psStatsAlloc(): This routine must create a new psStats data structure. *****************************************************************************/ -psStats *psStatsAlloc(psStatsOptions options) -{ - psStats *newStruct = NULL; - - newStruct = (psStats *) psAlloc(sizeof(psStats)); +psStats *psStatsAlloc( psStatsOptions options ) +{ + psStats * newStruct = NULL; + + newStruct = ( psStats * ) psAlloc( sizeof( psStats ) ); newStruct->sampleMean = NAN; newStruct->sampleMedian = NAN; @@ -1600,6 +1600,6 @@ newStruct->binsize = NAN; newStruct->options = options; - - return(newStruct); + + return ( newStruct ); } @@ -1616,54 +1616,54 @@ The histogram structure *****************************************************************************/ -psHistogram *psHistogramAlloc(float lower, - float upper, - int n) +psHistogram *psHistogramAlloc( float lower, + float upper, + int n ) { int i = 0; // Loop index variable psHistogram *newHist = NULL; // The new histogram structure float binSize = 0.0; // The histogram bin size - + // NOTE: Verify that this is the correct action. - if (n == 0) { - return(NULL); - } - if (n < 0) { - psAbort(__func__, "psHistogramAlloc() called with bin size %d.\n", n); - } - + if ( n == 0 ) { + return ( NULL ); + } + if ( n < 0 ) { + psAbort( __func__, "psHistogramAlloc() called with bin size %d.\n", n ); + } + // NOTE: Verify that this is the correct action. - if (lower > upper) { - return(NULL); - } - + if ( lower > upper ) { + return ( NULL ); + } + // Allocate memory for the new histogram structure. If there are N // bins, then there are N+1 bounds to those bins. - newHist = (psHistogram *) psAlloc(sizeof(psHistogram)); - p_psMemSetDeallocator(newHist,(psFreeFcn)histogramFree); - newHist->bounds = psVectorAlloc(n+1, PS_TYPE_F32); + newHist = ( psHistogram * ) psAlloc( sizeof( psHistogram ) ); + p_psMemSetDeallocator( newHist, ( psFreeFcn ) histogramFree ); + newHist->bounds = psVectorAlloc( n + 1, PS_TYPE_F32 ); newHist->bounds->n = newHist->bounds->nalloc; - + // Calculate the bounds for each bin. - binSize = (upper - lower) / (float) n; + binSize = ( upper - lower ) / ( float ) n; // NOTE: Is the following necessary? It prevents the max data point // from being in a non-existant bin. - binSize+= FLT_EPSILON; - for (i=0;ibounds->data.F32[i] = lower + (binSize * (float) i); - } - + binSize += FLT_EPSILON; + for ( i = 0;i < n + 1;i++ ) { + newHist->bounds->data.F32[ i ] = lower + ( binSize * ( float ) i ); + } + // Allocate the bins, and initialize them to zero. - newHist->nums = psVectorAlloc(n, PS_TYPE_U32); + newHist->nums = psVectorAlloc( n, PS_TYPE_U32 ); newHist->nums->n = newHist->nums->nalloc; - for (i=0;inums->n;i++) { - newHist->nums->data.U32[i] = 0; - } - + for ( i = 0;i < newHist->nums->n;i++ ) { + newHist->nums->data.U32[ i ] = 0; + } + // Initialize the other members. newHist->minNum = 0; newHist->maxNum = 0; newHist->uniform = true; - - return(newHist); + + return ( newHist ); } @@ -1677,55 +1677,55 @@ The histogram structure *****************************************************************************/ -psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds) -{ - psHistogram *newHist = NULL; // The new histogram structure +psHistogram *psHistogramAllocGeneric( const psVector *restrict bounds ) +{ + psHistogram * newHist = NULL; // The new histogram structure int i; // Loop index variable - + // NOTE: Verify that this is the correct action. - if (bounds == NULL) { - // psAbort(__func__, "psHistogram requested with NULL bounds"); - return(NULL); - } - + if ( bounds == NULL ) { + // psAbort(__func__, "psHistogram requested with NULL bounds"); + return ( NULL ); + } + // NOTE: Verify that this is the correct action. - if (bounds->n <= 1) { - // psAbort(__func__, "psHistogram requested with NULL bounds"); - return(NULL); - } - - if (bounds->type.type != PS_TYPE_F32) { - // psAbort(__func__, "psHistogram request a bound which is not type F32"); - return(NULL); - } - + if ( bounds->n <= 1 ) { + // psAbort(__func__, "psHistogram requested with NULL bounds"); + return ( NULL ); + } + + if ( bounds->type.type != PS_TYPE_F32 ) { + // psAbort(__func__, "psHistogram request a bound which is not type F32"); + return ( NULL ); + } + // Allocate memory for the new histogram structure. - newHist = (psHistogram *) psAlloc(sizeof(psHistogram)); - p_psMemSetDeallocator(newHist,(psFreeFcn)histogramFree); - newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32); + newHist = ( psHistogram * ) psAlloc( sizeof( psHistogram ) ); + p_psMemSetDeallocator( newHist, ( psFreeFcn ) histogramFree ); + 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]; - } - + for ( i = 0;i < bounds->n;i++ ) { + newHist->bounds->data.F32[ i ] = bounds->data.F32[ i ]; + } + // Allocate the bins, and initialize them to zero. If there are N bounds, // then there are N-1 bins. - newHist->nums = psVectorAlloc((bounds->n)-1, PS_TYPE_U32); + newHist->nums = psVectorAlloc( ( bounds->n ) - 1, PS_TYPE_U32 ); newHist->nums->n = newHist->nums->nalloc; - for (i=0;inums->n;i++) { - newHist->nums->data.U32[i] = 0; - } - + for ( i = 0;i < newHist->nums->n;i++ ) { + newHist->nums->data.U32[ i ] = 0; + } + // Initialize the other members. newHist->minNum = 0; newHist->maxNum = 0; newHist->uniform = false; - - return(newHist); -} - -static void histogramFree(psHistogram *myHist) -{ - psFree(myHist->bounds); - psFree(myHist->nums); + + return ( newHist ); +} + +static void histogramFree( psHistogram *myHist ) +{ + psFree( myHist->bounds ); + psFree( myHist->nums ); } @@ -1745,8 +1745,8 @@ The histogram structure "out". *****************************************************************************/ -psHistogram *psVectorHistogram(psHistogram *out, - const psVector *restrict in, - const psVector *restrict mask, - unsigned int maskVal) +psHistogram *psVectorHistogram( psHistogram *out, + const psVector *restrict in, + const psVector *restrict mask, + unsigned int maskVal ) { int i = 0; // Loop index variable @@ -1755,90 +1755,90 @@ int binNum = 0; // A temporary bin number int numBins = 0; // The total number of bins - + // NOTE: Verify that this is the correct action. - if (out == NULL) { - return(NULL); - } - + if ( out == NULL ) { + return ( NULL ); + } + // Check the specified output histogram for type psF32 - if (out->bounds->type.type != PS_TYPE_F32) { - psAbort(__func__, - "Only data type PS_TYPE_F32 for the output.bounds member."); - } - - if (out->nums->type.type != PS_TYPE_U32) { - psAbort(__func__, - "Only data type PS_TYPE_U32 for output.nums member."); - } - + if ( out->bounds->type.type != PS_TYPE_F32 ) { + psAbort( __func__, + "Only data type PS_TYPE_F32 for the output.bounds member." ); + } + + if ( out->nums->type.type != PS_TYPE_U32 ) { + psAbort( __func__, + "Only data type PS_TYPE_U32 for output.nums member." ); + } + // NOTE: Verify that this is the correct action. - if (in == NULL) { - return(out); - } - - if (in->type.type != PS_TYPE_F32) { - psAbort(__func__, - "Only data type PS_TYPE_F32 is currently supported (0x%x).", - in->type.type); - } - - if (mask != NULL) { - if (in->n != mask->n) { - psAbort(__func__, - "Vector data and vector mask are of different sizes."); - } - if (mask->type.type != PS_TYPE_U8) { - psAbort(__func__, "Vector mask must be type PS_TYPE_U8"); - } - } + if ( in == NULL ) { + return ( out ); + } + + if ( in->type.type != PS_TYPE_F32 ) { + psAbort( __func__, + "Only data type PS_TYPE_F32 is currently supported (0x%x).", + in->type.type ); + } + + if ( mask != NULL ) { + if ( in->n != mask->n ) { + psAbort( __func__, + "Vector data and vector mask are of different sizes." ); + } + if ( mask->type.type != PS_TYPE_U8 ) { + psAbort( __func__, "Vector mask must be type PS_TYPE_U8" ); + } + } // NOTE: determine the correct action for a variety of other cases: // in vector has 0 elements, and histogram structure has zero bins. - + numBins = out->nums->n; - for (i=0;in;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 (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); - - // NOTE: This next if-statement really shouldn't be necessary. - // However, do to numerical lack of precision, we occasionally - // produce a binNum outside the range of bins. - if (binNum >= out->nums->n) { - binNum = out->nums->n-1; - } - - (out->nums->data.S32[binNum])++; - - // If this is a non-uniform histogram, determining the correct - // bin number requires a bit more work. - } else { - // NOTE: 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])++; - } - } - } - } - } - } - return(out); + for ( i = 0;i < in->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 ( 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 ); + + // NOTE: This next if-statement really shouldn't be necessary. + // However, do to numerical lack of precision, we occasionally + // produce a binNum outside the range of bins. + if ( binNum >= out->nums->n ) { + binNum = out->nums->n - 1; + } + + ( out->nums->data.S32[ binNum ] ) ++; + + // If this is a non-uniform histogram, determining the correct + // bin number requires a bit more work. + } else { + // NOTE: 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 ] ) ++; + } + } + } + } + } + } + return ( out ); } @@ -1852,33 +1852,33 @@ the various stat functions. *****************************************************************************/ -psVector *p_psConvertToF32(psStats *stats, - psVector *in, - psVector *mask, - unsigned int maskVal) +psVector *p_psConvertToF32( psStats *stats, + psVector *in, + psVector *mask, + unsigned int maskVal ) { int i = 0; psVector *tmp = NULL; - - if (in->type.type == PS_TYPE_S32) { - tmp = psVectorAlloc(in->n, PS_TYPE_F32); - for (i=0;in;i++) { - tmp->data.F32[i] = (float) in->data.S32[i]; - } - } else if (in->type.type == PS_TYPE_U32) { - tmp = psVectorAlloc(in->n, PS_TYPE_F32); - for (i=0;in;i++) { - tmp->data.F32[i] = (float) in->data.U32[i]; - } - } else if (in->type.type == PS_TYPE_F64) { - tmp = psVectorAlloc(in->n, PS_TYPE_F32); - for (i=0;in;i++) { - tmp->data.F32[i] = (float) in->data.F64[i]; - } - } else if (in->type.type == PS_TYPE_F32) { - // do nothing - } else { - psAbort(__func__, "unsupported vector type 0x%x\n", in->type.type); - } - return(tmp); + + if ( in->type.type == PS_TYPE_S32 ) { + tmp = psVectorAlloc( in->n, PS_TYPE_F32 ); + for ( i = 0;i < in->n;i++ ) { + tmp->data.F32[ i ] = ( float ) in->data.S32[ i ]; + } + } else if ( in->type.type == PS_TYPE_U32 ) { + tmp = psVectorAlloc( in->n, PS_TYPE_F32 ); + for ( i = 0;i < in->n;i++ ) { + tmp->data.F32[ i ] = ( float ) in->data.U32[ i ]; + } + } else if ( in->type.type == PS_TYPE_F64 ) { + tmp = psVectorAlloc( in->n, PS_TYPE_F32 ); + for ( i = 0;i < in->n;i++ ) { + tmp->data.F32[ i ] = ( float ) in->data.F64[ i ]; + } + } else if ( in->type.type == PS_TYPE_F32 ) { + // do nothing + } else { + psAbort( __func__, "unsupported vector type 0x%x\n", in->type.type ); + } + return ( tmp ); } @@ -1901,93 +1901,93 @@ macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64. *****************************************************************************/ -psStats *psVectorStats(psStats *stats, - psVector *in, - psVector *mask, - unsigned int maskVal) -{ - psVector *inF32; +psStats *psVectorStats( psStats *stats, + psVector *in, + psVector *mask, + unsigned int maskVal ) +{ + psVector * inF32; int mustFreeTmp = 1; - + // NOTE: Verify that this is the correct action. - if (in == NULL) { - return(stats); - } - if (stats == NULL) { - return(NULL); - } - - inF32 = p_psConvertToF32(stats, in, mask, maskVal); - if (inF32 == NULL) { - inF32 = in; - mustFreeTmp = 0; - } - + if ( in == NULL ) { + return ( stats ); + } + if ( stats == NULL ) { + return ( NULL ); + } + + inF32 = p_psConvertToF32( stats, in, mask, maskVal ); + if ( inF32 == NULL ) { + inF32 = in; + mustFreeTmp = 0; + } + // XXX: Should we abort if (stats->min == stats->max)? - if ((stats->options & PS_STAT_USE_RANGE) && - (stats->min >= stats->max)) { - psAbort(__func__, "psVectorStats() called with range: %f to %f\n", - stats->min, stats->max); - } - + if ( ( stats->options & PS_STAT_USE_RANGE ) && + ( stats->min >= stats->max ) ) { + psAbort( __func__, "psVectorStats() called with range: %f to %f\n", + stats->min, stats->max ); + } + // PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32); - if (mask != NULL) { - PS_CHECK_NULL_VECTOR(mask); - PS_CHECK_EMPTY_VECTOR(mask); - PS_CHECK_VECTOR_SIZE_EQUAL(mask, in); - PS_CHECK_VECTOR_TYPE(mask, PS_TYPE_U8); - } - + if ( mask != NULL ) { + PS_CHECK_NULL_VECTOR( mask ); + PS_CHECK_EMPTY_VECTOR( mask ); + PS_CHECK_VECTOR_SIZE_EQUAL( mask, in ); + PS_CHECK_VECTOR_TYPE( mask, PS_TYPE_U8 ); + } + // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_MEAN) { - p_psVectorSampleMean(in, mask, maskVal, stats); - } - + if ( stats->options & PS_STAT_SAMPLE_MEAN ) { + p_psVectorSampleMean( in, mask, maskVal, stats ); + } + // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_MEDIAN) { - p_psVectorSampleMedian(in, mask, maskVal, stats); - } - + if ( stats->options & PS_STAT_SAMPLE_MEDIAN ) { + p_psVectorSampleMedian( in, mask, maskVal, stats ); + } + // ************************************************************************ // NOTE: 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); - } - + if ( stats->options & PS_STAT_SAMPLE_STDEV ) { + p_psVectorSampleMean( in, mask, maskVal, stats ); + p_psVectorSampleStdev( in, mask, maskVal, stats ); + } + // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_QUARTILE) { - p_psVectorSampleQuartiles(in, mask, maskVal, stats); - } - + if ( stats->options & PS_STAT_SAMPLE_QUARTILE ) { + p_psVectorSampleQuartiles( in, mask, maskVal, stats ); + } + // Since the various robust stats quantities share much computation, they // are grouped together in a single private function: // p_psVectorRobustStats() - if ((stats->options & PS_STAT_ROBUST_MEAN) || - (stats->options & PS_STAT_ROBUST_MEDIAN) || - (stats->options & PS_STAT_ROBUST_MODE) || - (stats->options & PS_STAT_ROBUST_STDEV) || - (stats->options & PS_STAT_ROBUST_QUARTILE)) { - p_psVectorRobustStats(in, mask, maskVal, stats); - } - - if ((stats->options & PS_STAT_CLIPPED_MEAN) || - (stats->options & PS_STAT_CLIPPED_STDEV)) { - p_psVectorClippedStats(in, mask, maskVal, stats); - } - + if ( ( stats->options & PS_STAT_ROBUST_MEAN ) || + ( stats->options & PS_STAT_ROBUST_MEDIAN ) || + ( stats->options & PS_STAT_ROBUST_MODE ) || + ( stats->options & PS_STAT_ROBUST_STDEV ) || + ( stats->options & PS_STAT_ROBUST_QUARTILE ) ) { + p_psVectorRobustStats( in, mask, maskVal, stats ); + } + + if ( ( stats->options & PS_STAT_CLIPPED_MEAN ) || + ( stats->options & PS_STAT_CLIPPED_STDEV ) ) { + p_psVectorClippedStats( in, mask, maskVal, stats ); + } + // ************************************************************************ - if (stats->options & PS_STAT_MAX) { - p_psVectorMax(in, mask, maskVal, stats); - } - + if ( stats->options & PS_STAT_MAX ) { + p_psVectorMax( in, mask, maskVal, stats ); + } + // ************************************************************************ - if (stats->options & PS_STAT_MIN) { - p_psVectorMin(in, mask, maskVal, stats); - } - - if (mustFreeTmp == 1) { - psFree(inF32); - } - return(stats); -} + if ( stats->options & PS_STAT_MIN ) { + p_psVectorMin( in, mask, maskVal, stats ); + } + + if ( mustFreeTmp == 1 ) { + psFree( inF32 ); + } + return ( stats ); +}