Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 1406) +++ trunk/psLib/src/math/psStats.c (revision 1407) @@ -1,2 +1,3 @@ + /** @file psStats.c * \brief basic statistical operations @@ -9,6 +10,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.51 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-06 22:34:05 $ + * @version $Revision: 1.52 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-07 00:06:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -21,6 +22,9 @@ #include #include + /*****************************************************************************/ + /* INCLUDE FILES */ + /*****************************************************************************/ #include "psMemory.h" @@ -35,8 +39,10 @@ /*****************************************************************************/ + /* DEFINE STATEMENTS */ + /*****************************************************************************/ // will use robust statistical methods. -#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. +#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. #define CLIPPED_NUM_ITER_LB 1 #define CLIPPED_NUM_ITER_UB 10 @@ -48,8 +54,6 @@ #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 */ @@ -83,10 +87,15 @@ printf("\n"); \ + /*****************************************************************************/ + /* TYPE DEFINITIONS */ + /*****************************************************************************/ /*****************************************************************************/ + /* GLOBAL VARIABLES */ + /*****************************************************************************/ @@ -94,5 +103,7 @@ /*****************************************************************************/ + /* FILE STATIC VARIABLES */ + /*****************************************************************************/ @@ -100,12 +111,13 @@ /*****************************************************************************/ + /* FUNCTION IMPLEMENTATION - LOCAL */ + /*****************************************************************************/ -bool p_psGetStatValue( const psStats* stats, double* value ) -{ - - switch ( stats->options & - ~ ( PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE ) ) { +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; @@ -163,16 +175,13 @@ *****************************************************************************/ -void p_psVectorPrint( psVector *myVector, - psVector *maskVector, - unsigned int maskVal, - psStats *stats ) -{ - int i = 0; // Loop index variable. - - 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 ] ); +void p_psVectorPrint(psVector * myVector, psVector * maskVector, unsigned int maskVal, psStats * stats) +{ + int i = 0; // Loop index variable. + + 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 ] ); + printf("Element %d is %f\n", i, myVector->data.F32[i]); } } @@ -195,5 +204,4 @@ *****************************************************************************/ - /****************************************************************************** p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the @@ -210,55 +218,51 @@ *****************************************************************************/ -void p_psVectorSampleMean( const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats ) -{ - int i = 0; // Loop index variable - float mean = 0.0; // The mean - int count = 0; // # of points in this mean? - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this +void p_psVectorSampleMean(const psVector * restrict myVector, + const psVector * restrict maskVector, unsigned int maskVal, psStats * stats) +{ + int i = 0; // Loop index variable + float mean = 0.0; // The mean + int 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 ) { + if (stats->options & PS_STAT_USE_RANGE) { rangeMin = stats->min; rangeMax = stats->max; - if ( maskVector != NULL ) { - for ( i = 0;i < myVector->n;i++ ) { + 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 ]; + 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; + 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 ]; + 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; + mean /= (float)count; } } else { - if ( maskVector != NULL ) { - for ( i = 0;i < myVector->n;i++ ) { - if ( !( maskVal & maskVector->data.U8[ i ] ) ) { - mean += myVector->data.F32[ i ]; + 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; + mean /= (float)count; } else { - for ( i = 0;i < myVector->n;i++ ) { - mean += myVector->data.F32[ i ]; - } - mean /= ( float ) myVector->n; + for (i = 0; i < myVector->n; i++) { + mean += myVector->data.F32[i]; + } + mean /= (float)myVector->n; } } @@ -278,50 +282,46 @@ NULL *****************************************************************************/ -void p_psVectorMax( const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats ) -{ - int i = 0; // Loop index variable - float max = -MY_MAX_FLOAT; // The calculated maximum - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this +void p_psVectorMax(const psVector * restrict myVector, + const psVector * restrict maskVector, unsigned int maskVal, psStats * stats) +{ + int i = 0; // Loop index variable + float max = -MY_MAX_FLOAT; // The calculated maximum + 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 ) { + 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 ]; + 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 ]; + 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 { - 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 ]; + 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 ]; + for (i = 0; i < myVector->n; i++) { + if (myVector->data.F32[i] > max) { + max = myVector->data.F32[i]; } } @@ -343,49 +343,45 @@ NULL *****************************************************************************/ -void p_psVectorMin( const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats ) -{ - int i = 0; // Loop index variable - float min = MY_MAX_FLOAT; // The calculated maximum - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - if ( stats->options & PS_STAT_USE_RANGE ) { +void p_psVectorMin(const psVector * restrict myVector, + const psVector * restrict maskVector, unsigned int maskVal, psStats * stats) +{ + int i = 0; // Loop index variable + float min = MY_MAX_FLOAT; // The calculated maximum + 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 ] ) ) { - if ( ( myVector->data.F32[ i ] < min ) && - ( rangeMin <= myVector->data.F32[ i ] ) && - ( myVector->data.F32[ i ] <= rangeMax ) ) { - 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) && + (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 ]; + 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 { - 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 ]; + 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 ]; + for (i = 0; i < myVector->n; i++) { + if (myVector->data.F32[i] < min) { + min = myVector->data.F32[i]; } } @@ -408,29 +404,25 @@ NULL *****************************************************************************/ -int p_psVectorNValues( const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats ) -{ - int i = 0; // Loop index variable - int 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 ) { +int p_psVectorNValues(const psVector * restrict myVector, + const psVector * restrict maskVector, unsigned int maskVal, psStats * stats) +{ + int i = 0; // Loop index variable + int 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 ) ) { + 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 ) ) { + for (i = 0; i < myVector->n; i++) { + if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { numData++; } @@ -440,7 +432,7 @@ rangeMin = stats->min; rangeMax = stats->max; - if ( maskVector != NULL ) { - for ( i = 0;i < myVector->n;i++ ) { - if ( !( maskVal & maskVector->data.U8[ i ] ) ) { + if (maskVector != NULL) { + for (i = 0; i < myVector->n; i++) { + if (!(maskVal & maskVector->data.U8[i])) { numData++; } @@ -450,8 +442,6 @@ } } - return ( numData ); -} - - + return (numData); +} /****************************************************************************** @@ -466,17 +456,14 @@ NULL *****************************************************************************/ -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 - int count = 0; // # of points in this mean? - int nValues = 0; // # of points in vector - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - +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 + int count = 0; // # of points in this mean? + int nValues = 0; // # of points in vector + 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 @@ -485,37 +472,33 @@ // regardless of the vector size. /* - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleMedian = stats2->robustMedian; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } - */ + * if (myVector->n > stats->sampleLimit) { psAbort(__func__, "Robust Statistic Algorithms have not yet + * been defined or implemented."); + * + * // Calculate the robust quartiles. stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + * p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); + * + * // Store the robust quartiles into the sample quartile members. stats->sampleMedian = + * stats2->robustMedian; + * + * // Free temporary data buffers. psFree(stats2); + * + * // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. stats->options = stats->options | + * PS_STAT_ROBUST_FOR_SAMPLE; + * + * return; } */ // Determine how many data points fit inside this min/max range // 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 ) { + if (stats->options & PS_STAT_USE_RANGE) { rangeMin = stats->min; rangeMax = stats->max; @@ -524,17 +507,15 @@ // 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 ]; + 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 ]; + 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]; } } @@ -543,31 +524,31 @@ // 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 ]; + 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 ]; + 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 ] ); + 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 ]; + stats->sampleMedian = sortedVector->data.F32[nValues / 2]; } // Free the temporary data structures. - psFree( unsortedVector ); - psFree( sortedVector ); + psFree(unsortedVector); + psFree(sortedVector); } @@ -584,27 +565,26 @@ XXX: use a static variable for gaussianCoefs[] and compute them once. *****************************************************************************/ -psVector *p_psVectorsmoothHistGaussian( psHistogram *robustHistogram, - float sigma ) -{ - int i = 0; // Loop index variable - int j = 0; // Loop index variable - 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 ) { +psVector *p_psVectorsmoothHistGaussian(psHistogram * robustHistogram, float sigma) +{ + int i = 0; // Loop index variable + int j = 0; // Loop index variable + 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 ); + denom = (float)sqrt(2.0 * M_PI * sigma * sigma); #else - denom = sqrtf( 2.0 * M_PI * sigma * sigma ); + 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 ); + 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 @@ -612,29 +592,27 @@ // 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]); + 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 ); + /* 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;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 ); + 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); } @@ -650,29 +628,27 @@ NULL *****************************************************************************/ -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 - int count = 0; // # of points in this mean? - int nValues = 0; // # data points - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this +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 + int count = 0; // # of points in this mean? + int nValues = 0; // # data points + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this // Determine how many data points fit inside this min/max range // and are not maxed, IF the maskVector is not NULL> - 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 ) { + if (stats->options & PS_STAT_USE_RANGE) { rangeMin = stats->min; rangeMax = stats->max; @@ -680,17 +656,15 @@ // 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 ]; + 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 ]; + 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]; } } @@ -699,13 +673,13 @@ // 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 ]; + 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 ]; + for (i = 0; i < myVector->n; i++) { + unsortedVector->data.F32[i] = myVector->data.F32[i]; } } @@ -713,15 +687,14 @@ // 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); +} /****************************************************************************** @@ -737,34 +710,31 @@ *****************************************************************************/ -void p_psVectorSampleStdev( const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats ) -{ - int i = 0; // Loop index variable - int countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this +void p_psVectorSampleStdev(const psVector * restrict myVector, + const psVector * restrict maskVector, unsigned int maskVal, psStats * stats) +{ + int i = 0; // Loop index variable + int countInt = 0; // # of data points being used + float countFloat = 0.0; // # of data points being used + float mean = 0.0; // The mean + float diff = 0.0; // Used in calculating stdev + float sumSquares = 0.0; // temporary variable + float sumDiffs = 0.0; // temporary variable + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this // This procedure requires the mean. If it has not been already // 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;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 ); + 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++; @@ -772,9 +742,8 @@ } } 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 ); + 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++; @@ -784,9 +753,9 @@ } } else { - 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 ); + 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++; @@ -794,7 +763,7 @@ } } else { - for ( i = 0;i < myVector->n;i++ ) { - diff = myVector->data.F32[ i ] - mean; - sumSquares += ( diff * diff ); + for (i = 0; i < myVector->n; i++) { + diff = myVector->data.F32[i] - mean; + sumSquares += (diff * diff); sumDiffs += diff; countInt++; @@ -803,14 +772,12 @@ } } - countFloat = ( float ) countInt; + 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 } @@ -828,52 +795,43 @@ NULL *****************************************************************************/ -void p_psVectorClippedStats( const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats ) -{ - int i = 0; // Loop index variable - int j = 0; // Loop index variable - float clippedMean = 0.0; // self-explanatory - float clippedStdev = 0.0; // self-explanatory - float oldStanMean = 0.0; // Temporary variable - float oldStanStdev = 0.0; // Temporary variable - psVector *tmpMask = NULL; // Temporary vector +void p_psVectorClippedStats(const psVector * restrict myVector, + const psVector * restrict maskVector, unsigned int maskVal, psStats * stats) +{ + int i = 0; // Loop index variable + int j = 0; // Loop index variable + float clippedMean = 0.0; // self-explanatory + float clippedStdev = 0.0; // self-explanatory + float oldStanMean = 0.0; // Temporary variable + 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;i < tmpMask->n;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. @@ -889,14 +847,13 @@ // 5. Repeat N times: - for ( i = 0;i < stats->clipIter;i++ ) { - for ( j = 0;j < myVector->n;j++ ) { + 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; + 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 ); + p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); + p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); // c) Use the new mean for x @@ -912,14 +869,13 @@ // 7. The last calcuated value of x is the cliped mean. - if ( stats->options & PS_STAT_CLIPPED_MEAN ) { + 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 ) { + if (stats->options & PS_STAT_CLIPPED_STDEV) { stats->clippedStdev = clippedStdev; } - psFree( tmpMask ); + psFree(tmpMask); } @@ -928,26 +884,25 @@ 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;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 ]; + 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;i < myData->n;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; + } +} /***************************************************************************** @@ -956,15 +911,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 ) ) ); - - // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); - return ( tmp ); +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); } @@ -973,28 +928,24 @@ 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 ) { + 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 ); -} - + 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); +} /***************************************************************************** @@ -1003,15 +954,14 @@ 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); } @@ -1020,22 +970,19 @@ 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 ) { + if (whichParamDeriv == 0) { tmp = x * x; - } else - if ( whichParamDeriv == 1 ) { - tmp = x; - } else - if ( whichParamDeriv == 2 ) { - tmp = 1.0; - } - - return ( tmp ); + } else if (whichParamDeriv == 1) { + tmp = x; + } else if (whichParamDeriv == 2) { + tmp = 1.0; + } + + return (tmp); } @@ -1049,8 +996,5 @@ decreasing within that range. *****************************************************************************/ -float p_ps1DPolyMedian( psPolynomial1D *myPoly, - float rangeLow, - float rangeHigh, - float getThisValue ) +float p_ps1DPolyMedian(psPolynomial1D * myPoly, float rangeLow, float rangeHigh, float getThisValue) { int numIterations = 0; @@ -1059,20 +1003,20 @@ 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 ); + // 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 ) { + 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 { @@ -1081,5 +1025,5 @@ numIterations++; } - return ( midpoint ); + return (midpoint); } @@ -1092,52 +1036,42 @@ 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 = psVectorFitPolynomial1D( myPoly, x, y, yErr ); - return ( p_ps1DPolyMedian( myPoly, x->data.F64[ 0 ], x->data.F64[ 2 ], - fitFloat ) ); +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 ); + 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); } @@ -1166,28 +1100,27 @@ 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 - int LQBinNum = -1; // Bin num for lower quartile - int UQBinNum = -1; // Bin num for upper quartile - int i = 0; // Loop index variable + float binSize = 0.0; // Size of the histogram bins + int LQBinNum = -1; // Bin num for lower quartile + int UQBinNum = -1; // Bin num for upper quartile + int i = 0; // Loop index variable int maxBinNum = 0; float maxBinCount = 0.0; float dL = 0.0; int numBins = 0; - psStats *tmpStats = psStatsAlloc( PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN ); - // psImage *domain; - // psVector *errors; - // psVector *data; - // psVector *initialGuess; - // psVector *theParams; - // float chiSq=0.0; - // float max = -HUGE; - // float max_pos; + psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN); + + // psImage *domain; + // psVector *errors; + // psVector *data; + // psVector *initialGuess; + // psVector *theParams; + // float chiSq=0.0; + // float max = -HUGE; + // float max_pos; float myMean = 0.0; float myStdev = 0.0; @@ -1196,8 +1129,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; @@ -1210,23 +1143,23 @@ // 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 ) { + if (fabs(binSize) <= FLT_EPSILON) { + if (stats->options & PS_STAT_ROBUST_MEAN) { stats->robustMean = stats->clippedMean; } - if ( stats->options & PS_STAT_ROBUST_MEDIAN ) { + if (stats->options & PS_STAT_ROBUST_MEDIAN) { stats->robustMedian = stats->clippedMean; } - if ( stats->options & PS_STAT_ROBUST_MODE ) { + if (stats->options & PS_STAT_ROBUST_MODE) { stats->robustMode = stats->clippedMean; } - if ( stats->options & PS_STAT_ROBUST_STDEV ) { + if (stats->options & PS_STAT_ROBUST_STDEV) { stats->robustStdev = 0.0; } - if ( stats->options & PS_STAT_ROBUST_QUARTILE ) { + if (stats->options & PS_STAT_ROBUST_QUARTILE) { stats->robustUQ = stats->clippedMean; stats->robustLQ = stats->clippedMean; @@ -1235,32 +1168,27 @@ stats->robustNfit = 0.0; stats->robustN50 = 0.0; - psFree( tmpStats ); - return ; - } - + 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 @@ -1274,17 +1202,16 @@ // 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;i < robustHistogramVector->n;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;i < cumulativeRobustSumsFullRange->n;i++ ) { - if ( cumulativeRobustSumsFullRange->data.F32[ i ] >= ( sumRobust / 4.0 ) ) { + for (i = 0; i < cumulativeRobustSumsFullRange->n; i++) { + if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust / 4.0)) { LQBinNum = i; break; @@ -1294,6 +1221,6 @@ // 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 ) ) { + for (i = cumulativeRobustSumsFullRange->n - 1; i >= 0; i--) { + if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0 * sumRobust / 4.0)) { UQBinNum = i; break; @@ -1301,8 +1228,8 @@ } - if ( ( LQBinNum == -1 ) || - ( UQBinNum == -1 ) ) { - psAbort( __func__, "Could not determine the robust lower/upper quartiles." ); - } + if ((LQBinNum == -1) || (UQBinNum == -1)) { + psAbort(__func__, "Could not determine the robust lower/upper quartiles."); + } + /************************************************************************** Determine the mode in the range LQ:UQ. @@ -1310,35 +1237,34 @@ // 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 ) { + 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[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;i < robustHistogramVector->n;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]; } } @@ -1348,10 +1274,9 @@ // 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 ]; + 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]; } } @@ -1361,31 +1286,29 @@ // 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 ) { + if (stats->options & PS_STAT_ROBUST_MEAN) { stats->robustMean = myMean; } - if ( stats->options & PS_STAT_ROBUST_MODE ) { + 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 ) { + (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. @@ -1394,49 +1317,46 @@ // 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 ]; + 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" ); + 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; + 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 ); + 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 ); + 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 ] ); - } - } - + 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 @@ -1444,78 +1364,72 @@ // 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 = psVectorFitPolynomial1D( myPoly, x, y, yErr ); - stats->robustLQ = p_ps1DPolyMedian( myPoly, - x->data.F64[ 0 ], - x->data.F64[ 2 ], - countFloat / 4.0 ); + 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 ); + 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->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum + 1] + + robustHistogram->bounds->data.F32[UQBinNum]); } } @@ -1523,15 +1437,13 @@ 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); +} /* @@ -1582,19 +1494,20 @@ */ - /*****************************************************************************/ + /* FUNCTION IMPLEMENTATION - PUBLIC */ + /*****************************************************************************/ -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; @@ -1620,5 +1533,5 @@ newStruct->options = options; - return ( newStruct ); + return (newStruct); } @@ -1635,46 +1548,42 @@ The histogram structure *****************************************************************************/ -psHistogram *psHistogramAlloc( float lower, - float upper, - int n ) -{ - int i = 0; // Loop index variable +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 + 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; - // NOTE: Is the following necessary? It prevents the max data point + 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;i < n + 1;i++ ) { - newHist->bounds->data.F32[ i ] = lower + ( binSize * ( float ) i ); + 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;i < newHist->nums->n;i++ ) { - newHist->nums->data.U32[ i ] = 0; + for (i = 0; i < newHist->nums->n; i++) { + newHist->nums->data.U32[i] = 0; } @@ -1684,5 +1593,5 @@ newHist->uniform = true; - return ( newHist ); + return (newHist); } @@ -1696,41 +1605,39 @@ The histogram structure *****************************************************************************/ -psHistogram *psHistogramAllocGeneric( const psVector *restrict bounds ) -{ - psHistogram * newHist = NULL; // The new histogram structure - int i; // Loop index variable +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 ) { + if (bounds == NULL) { // psAbort(__func__, "psHistogram requested with NULL bounds"); - return ( NULL ); - } - + return (NULL); + } // NOTE: Verify that this is the correct action. - if ( bounds->n <= 1 ) { + if (bounds->n <= 1) { // psAbort(__func__, "psHistogram requested with NULL bounds"); - return ( NULL ); - } - - if ( bounds->type.type != PS_TYPE_F32 ) { + return (NULL); + } + + if (bounds->type.type != PS_TYPE_F32) { // psAbort(__func__, "psHistogram request a bound which is not type F32"); - return ( NULL ); - } - + 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;i < bounds->n;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;i < newHist->nums->n;i++ ) { - newHist->nums->data.U32[ i ] = 0; + for (i = 0; i < newHist->nums->n; i++) { + newHist->nums->data.U32[i] = 0; } @@ -1740,13 +1647,12 @@ 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); +} /***************************************************************************** @@ -1764,49 +1670,41 @@ The histogram structure "out". *****************************************************************************/ -psHistogram *psVectorHistogram( psHistogram *out, - const psVector *restrict in, - const psVector *restrict mask, - unsigned int maskVal ) -{ - int i = 0; // Loop index variable - int j = 0; // Loop index variable - float binSize = 0.0; // Histogram bin size - int binNum = 0; // A temporary bin number - int numBins = 0; // The total number of bins +psHistogram *psVectorHistogram(psHistogram * out, + const psVector * restrict in, + const psVector * restrict mask, unsigned int maskVal) +{ + int i = 0; // Loop index variable + int j = 0; // Loop index variable + float binSize = 0.0; // Histogram bin size + 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"); } } @@ -1815,50 +1713,47 @@ numBins = out->nums->n; - for ( i = 0;i < in->n;i++ ) { + 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 ) ) ) ) { + 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 ] ) { + 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 (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 { - // 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 ] ) ++; - } + // 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 ); + } + } + } + return (out); } @@ -1872,38 +1767,31 @@ 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;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 ); -} - + 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); +} /****************************************************************************** @@ -1924,93 +1812,77 @@ 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 ) { + 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 ); - } - - // 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 ((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 ( 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 ); - } - + // 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_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); +}