Index: trunk/psLib/src/dataManip/psStats.c =================================================================== --- trunk/psLib/src/dataManip/psStats.c (revision 1361) +++ trunk/psLib/src/dataManip/psStats.c (revision 1385) @@ -36,33 +36,33 @@ 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"); \ @@ -92,51 +92,51 @@ 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; - } + 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; + } } @@ -153,11 +153,11 @@ { 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 ] ); - } + 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 ] ); + } } @@ -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;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; - } + 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 { - 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; - } - } - + 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 { + 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; } @@ -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;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 ]; - } - } - } + 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 { - 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 ) && + ( 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 ]; + } + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + if ( myVector->data.F32[ i ] > max ) { + max = myVector->data.F32[ i ]; + } + } + } + } + stats->max = max; } @@ -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;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 ]; - } - } - } + 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 { - 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 ) && + ( 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 ]; + } + } + } + } else { + for ( i = 0;i < myVector->n;i++ ) { + if ( myVector->data.F32[ i ] < min ) { + min = myVector->data.F32[ i ]; + } + } + } + } + stats->min = min; } @@ -401,37 +401,37 @@ 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 ) ) { - numData++; - } - } - } else { - for ( i = 0;i < myVector->n;i++ ) { - if ( ( rangeMin <= myVector->data.F32[ i ] ) && - ( myVector->data.F32[ i ] <= rangeMax ) ) { - numData++; - } - } - } + 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 { - 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; - } - } + for ( i = 0;i < myVector->n;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;i < myVector->n;i++ ) { + if ( !( maskVal & maskVector->data.U8[ i ] ) ) { + numData++; + } + } + } else { + numData = myVector->n; + } + } return ( numData ); } @@ -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,67 +488,67 @@ } */ - + // 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 ); - + // Allocate temporary vectors for the data. unsortedVector = psVectorAlloc( nValues, PS_TYPE_F32 ); unsortedVector->n = unsortedVector->nalloc; - + 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;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 ]; - } - } - } + 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 { - // 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 ]; - } - } - } + 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 { + // 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 ); - + // 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 ]; - } - + 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 ); @@ -577,47 +577,47 @@ 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 ); - } - } - + 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;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 ] ); - } - } - } + 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 ); } @@ -646,9 +646,9 @@ 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 ); - + // Allocate temporary vectors for the data. unsortedVector = psVectorAlloc( nValues, PS_TYPE_F32 ); @@ -656,51 +656,51 @@ 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;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 ]; - } - } - } + 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 { - // 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 ]; - } - } - } - + 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 { + // 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 ); - + // Calculate the quartile points exactly. stats->sampleUQ = sortedVector->data.F32[ 3 * ( nValues / 4 ) ]; stats->sampleLQ = sortedVector->data.F32[ nValues / 4 ]; - + // Free the temporary data structures. psFree( unsortedVector ); @@ -735,64 +735,64 @@ 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 ); - } + 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 ); - 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; - } + 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 { - 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; - } - } + 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 { + 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 ) ); #else - + stats->sampleStdev = sqrtf( ( sumSquares - ( sumDiffs * sumDiffs / countFloat ) ) / ( countFloat - 1 ) ); @@ -824,19 +824,19 @@ 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 ); - } - + 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 ); - } - + 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. @@ -844,65 +844,65 @@ 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 ]; - } - } - + 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 ); - + // 2. Compute the sample standard deviation. 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;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; - } - - } + 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; - + // 7. The last calcuated value of x is the cliped mean. if ( stats->options & PS_STAT_CLIPPED_MEAN ) { - stats->clippedMean = clippedMean; - } - + stats->clippedMean = clippedMean; + } + // 8. The last calcuated value of stdev is the cliped stdev. if ( stats->options & PS_STAT_CLIPPED_STDEV ) { - stats->clippedStdev = clippedStdev; - } - + stats->clippedStdev = clippedStdev; + } + psFree( tmpMask ); } @@ -918,18 +918,18 @@ 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 ]; - } - } - + 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; - } + myData->data.F32[ i ] = ( myData->data.F32[ i ] - min ) / range; + } } @@ -948,5 +948,5 @@ 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 ); @@ -965,10 +965,11 @@ 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 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 ); @@ -976,5 +977,5 @@ } printf( "p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp ); - + return ( tmp ); } @@ -994,5 +995,5 @@ float C = myParams->data.F32[ 2 ]; float tmp = 0.0; - + tmp = ( A * x * x ) + ( B * x ) + C; return ( tmp ); @@ -1009,13 +1010,15 @@ float x = myCoords->data.F32[ 0 ]; float tmp = 0.0; - + if ( whichParamDeriv == 0 ) { - tmp = x * x; - } else if ( whichParamDeriv == 1 ) { + tmp = x * x; + } else + if ( whichParamDeriv == 1 ) { tmp = x; - } else if ( whichParamDeriv == 2 ) { - tmp = 1.0; - } - + } else + if ( whichParamDeriv == 2 ) { + tmp = 1.0; + } + return ( tmp ); } @@ -1039,27 +1042,27 @@ 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++; - } + 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 ); } @@ -1082,38 +1085,38 @@ 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 ] ) ); - } - + 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 ); @@ -1187,5 +1190,5 @@ 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 @@ -1193,39 +1196,39 @@ 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 ( 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 ); - } + p_psVectorMin( myVector, maskVector, maskVal, stats ); + } if ( isnan( stats->max ) ) { - p_psVectorMax( myVector, maskVector, maskVal, stats ); - } - + 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 @@ -1236,20 +1239,20 @@ stats->max, numBins ); - + // Populate the histogram array. psVectorHistogram( robustHistogram, myVector, maskVector, maskVal ); - + // Smooth the histogram. // XXX: is that the right stdev? 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 @@ -1258,35 +1261,35 @@ 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 ]; - } + 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 ) ) { - LQBinNum = i; - break; - } - } - + 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 ( 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." ); - } + 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; @@ -1294,16 +1297,16 @@ 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 ]; - } - + 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 ); - + /************************************************************************** Determine the mean/stdev for the bins in the range mode-dL to mode+dL @@ -1311,18 +1314,18 @@ cumulativeRobustSumsDlRange = psVectorAlloc( robustHistogramVector->n, PS_TYPE_F32 ); for ( i = 0;i < robustHistogramVector->n;i++ ) { - cumulativeRobustSumsDlRange->data.F32[ i ] = 0.0; - } + 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 ]; - } - } - + 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 @@ -1330,43 +1333,43 @@ 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 ]; - } - } + 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 ]; - } - } + 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; - } - + 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 ] ); - } - + 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; - } - + 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 +1377,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 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 ] ); - } - } - + 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 @@ -1426,82 +1429,82 @@ // 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 ] ); - } - } + 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 ); @@ -1576,5 +1579,5 @@ { psStats * newStruct = NULL; - + newStruct = ( psStats * ) psAlloc( sizeof( psStats ) ); newStruct->sampleMean = NAN; @@ -1600,5 +1603,5 @@ newStruct->binsize = NAN; newStruct->options = options; - + return ( newStruct ); } @@ -1623,18 +1626,18 @@ 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 ); - } + return ( NULL ); + } if ( n < 0 ) { - psAbort( __func__, "psHistogramAlloc() called with bin size %d.\n", n ); - } - + psAbort( __func__, "psHistogramAlloc() called with bin size %d.\n", n ); + } + // NOTE: Verify that this is the correct action. if ( lower > upper ) { - return ( NULL ); - } - + return ( NULL ); + } + // Allocate memory for the new histogram structure. If there are N // bins, then there are N+1 bounds to those bins. @@ -1643,5 +1646,5 @@ 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; @@ -1650,19 +1653,19 @@ binSize += FLT_EPSILON; for ( i = 0;i < n + 1;i++ ) { - newHist->bounds->data.F32[ i ] = lower + ( binSize * ( float ) 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->n = newHist->nums->nalloc; for ( i = 0;i < newHist->nums->n;i++ ) { - newHist->nums->data.U32[ i ] = 0; - } - + newHist->nums->data.U32[ i ] = 0; + } + // Initialize the other members. newHist->minNum = 0; newHist->maxNum = 0; newHist->uniform = true; - + return ( newHist ); } @@ -1681,22 +1684,22 @@ 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 ); - } - + // 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 ); - } - + // 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 ); - } - + // 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 ) ); @@ -1705,7 +1708,7 @@ newHist->bounds->n = newHist->bounds->nalloc; for ( i = 0;i < bounds->n;i++ ) { - newHist->bounds->data.F32[ i ] = bounds->data.F32[ 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. @@ -1713,12 +1716,12 @@ newHist->nums->n = newHist->nums->nalloc; for ( i = 0;i < newHist->nums->n;i++ ) { - newHist->nums->data.U32[ i ] = 0; - } - + newHist->nums->data.U32[ i ] = 0; + } + // Initialize the other members. newHist->minNum = 0; newHist->maxNum = 0; newHist->uniform = false; - + return ( newHist ); } @@ -1755,89 +1758,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 ); - } - + 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." ); - } - + 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." ); - } - + 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 ); - } - + 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__, - "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" ); - } - } + "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;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 ] ) ++; - } - } - } + // 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 ); } @@ -1859,25 +1863,28 @@ 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.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.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 ); - } + } else + if ( in->type.type == PS_TYPE_F32 ) { + // do nothing + } else { + psAbort( __func__, "unsupported vector type 0x%x\n", in->type.type ); + } return ( tmp ); } @@ -1908,57 +1915,57 @@ psVector * inF32; int mustFreeTmp = 1; - + // NOTE: Verify that this is the correct action. if ( in == NULL ) { - return ( stats ); - } + return ( stats ); + } if ( stats == NULL ) { - return ( NULL ); - } - + return ( NULL ); + } + inF32 = p_psConvertToF32( stats, in, mask, maskVal ); if ( inF32 == NULL ) { - inF32 = in; - mustFreeTmp = 0; - } - + 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 ); - } - + 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 ); - } - + 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 ); - } - + p_psVectorSampleMean( in, mask, maskVal, stats ); + } + // ************************************************************************ if ( stats->options & PS_STAT_SAMPLE_MEDIAN ) { - p_psVectorSampleMedian( in, mask, maskVal, stats ); - } - + 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 ); - } - + 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 ); - } - + p_psVectorSampleQuartiles( in, mask, maskVal, stats ); + } + // Since the various robust stats quantities share much computation, they // are grouped together in a single private function: @@ -1969,25 +1976,25 @@ ( stats->options & PS_STAT_ROBUST_STDEV ) || ( stats->options & PS_STAT_ROBUST_QUARTILE ) ) { - p_psVectorRobustStats( in, mask, maskVal, stats ); - } - + 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 ); - } - + p_psVectorClippedStats( in, mask, maskVal, stats ); + } + // ************************************************************************ if ( stats->options & PS_STAT_MAX ) { - p_psVectorMax( in, mask, maskVal, stats ); - } - + p_psVectorMax( in, mask, maskVal, stats ); + } + // ************************************************************************ if ( stats->options & PS_STAT_MIN ) { - p_psVectorMin( in, mask, maskVal, stats ); - } - + p_psVectorMin( in, mask, maskVal, stats ); + } + if ( mustFreeTmp == 1 ) { - psFree( inF32 ); - } + psFree( inF32 ); + } return ( stats ); }