Index: /trunk/psLib/src/math/psMathUtils.c =================================================================== --- /trunk/psLib/src/math/psMathUtils.c (revision 7131) +++ /trunk/psLib/src/math/psMathUtils.c (revision 7132) @@ -1,8 +1,8 @@ /** @file psMathUtils.c * - * This file contains standard math rotines. + * This file contains standard math routines. * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-02 21:09:07 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-18 01:20:51 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -48,31 +48,25 @@ /***************************************************************************** -vectorBinDisectTYPE(): This is a macro for a private function which takes as -input an array of data as well as a single value for that data. The input -vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i). -This routine does a binary disection of the vector and returns "i" such that -(v[i] <= x <= v[i+1). If x lies outside the range of v[], then this routine -prints a warning message and returns (-2 or -1). +This is a macro covering the various types for the below function. *****************************************************************************/ -#define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \ -static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \ - psS32 numBins, \ - ps##TYPE x) \ -{ \ - psS32 min; \ - psS32 max; \ - psS32 mid; \ +#define VECTOR_BIN_DISECT_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + ps##TYPE *bounds = bins->data.TYPE; \ + ps##TYPE value = x->data.TYPE; \ + long min; \ + long max; \ + long mid; \ psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \ /* psTrace(__func__, 6, "Determining the bin for: %f\n", x); */\ - if (x < bins[0]) { \ + if (value < bounds[0]) { \ psLogMsg(__func__, PS_LOG_WARN, \ "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ - #TYPE, x, bins[0], bins[numBins-1]); \ + #TYPE, value, bounds[0], bounds[numBins-1]); \ return(-2); \ } \ - if (x > bins[numBins-1]) { \ + if (value > bounds[numBins-1]) { \ psLogMsg(__func__, PS_LOG_WARN, \ "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ - #TYPE, x, bins[0], bins[numBins-1]); \ + #TYPE, value, bounds[0], bounds[numBins-1]); \ return(-1); \ } \ @@ -82,11 +76,9 @@ mid = ((max+1)-min)/2; \ while (min != max) { \ - /*psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins[mid]) is (%f, %f)\n", min, mid, max, x, bins[mid]); */\ \ - if (x == bins[mid]) { \ - /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[mid], bins[mid+1]); */\ + if (value == bounds[mid]) { \ psTrace(__func__, 4, "---- %s(%d) end (1) ----\n", __func__, mid); \ return(mid); \ - } else if (x < bins[mid]) { \ + } else if (value < bounds[mid]) { \ max = mid-1; \ } else { \ @@ -95,26 +87,17 @@ mid = ((max+1)+min)/2; \ } \ - /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[min], bins[min+1]);*/ \ psTrace(__func__, 4, "---- %s(%d) end (2) ----\n", __func__, min); \ return(min); \ -} \ - -FUNC_MACRO_VECTOR_BIN_DISECT(S8) -FUNC_MACRO_VECTOR_BIN_DISECT(S16) -FUNC_MACRO_VECTOR_BIN_DISECT(S32) -FUNC_MACRO_VECTOR_BIN_DISECT(S64) -FUNC_MACRO_VECTOR_BIN_DISECT(U8) -FUNC_MACRO_VECTOR_BIN_DISECT(U16) -FUNC_MACRO_VECTOR_BIN_DISECT(U32) -FUNC_MACRO_VECTOR_BIN_DISECT(U64) -FUNC_MACRO_VECTOR_BIN_DISECT(F32) -FUNC_MACRO_VECTOR_BIN_DISECT(F64) +} /***************************************************************************** -p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). +p_psVectorBinDisect(): This function takes as input an array of data as well as a single value for that data. +The input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i). This routine does a +binary disection of the vector and returns "i" such that (v[i] <= x <= v[i+1). If x lies outside the range of +v[], then this routine prints a warning message and returns (-2 or -1). *****************************************************************************/ psS32 p_psVectorBinDisect( - psVector *bins, - psScalar *x) + const psVector *bins, + const psScalar *x) { PS_ASSERT_VECTOR_NON_NULL(bins, -4); @@ -122,31 +105,23 @@ PS_ASSERT_PTR_NON_NULL(x, -6); PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3); - char* strType; + long numBins = bins->n; // Number of bins switch (x->type.type) { - case PS_TYPE_U8: - return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8)); - case PS_TYPE_U16: - return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16)); - case PS_TYPE_U32: - return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32)); - case PS_TYPE_U64: - return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64)); - case PS_TYPE_S8: - return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8)); - case PS_TYPE_S16: - return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16)); - case PS_TYPE_S32: - return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); - case PS_TYPE_S64: - return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64)); - case PS_TYPE_F32: - return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); - case PS_TYPE_F64: - return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64)); - default: - PS_TYPE_NAME(strType, x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); - return -1; + VECTOR_BIN_DISECT_CASE(S8); + VECTOR_BIN_DISECT_CASE(S16); + VECTOR_BIN_DISECT_CASE(S32); + VECTOR_BIN_DISECT_CASE(S64); + VECTOR_BIN_DISECT_CASE(U8); + VECTOR_BIN_DISECT_CASE(U16); + VECTOR_BIN_DISECT_CASE(U32); + VECTOR_BIN_DISECT_CASE(U64); + VECTOR_BIN_DISECT_CASE(F32); + VECTOR_BIN_DISECT_CASE(F64); + default: { + char* strType; + PS_TYPE_NAME(strType, x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); + return -1; + } } return(-3); @@ -154,23 +129,9 @@ - -/***************************************************************************** -fullInterpolateTYPE(): This routine will take as input n-element psVectors -domain and range, and the x value, assumed to lie with the domain vector. It -produces as output the (order-1)-order LaGrange interpolated value of x. - -XXX: do we error check for non-distinct domain values? - -This routine will only be called inside this file. So, we do not have to -do PS_ASSERTS. - *****************************************************************************/ -#define FUNC_MACRO_FULL_INTERPOLATE(TYPE) \ -static psScalar *vectorInterpolate##TYPE( \ - psScalar *out, \ - psVector *domain, \ - psVector *range, \ - psS32 order, \ - psScalar *x) \ -{ \ +/************************************************************************************************************* +Helper macro for p_psVectorInterpolate: handles each of the types +*************************************************************************************************************/ +#define VECTOR_INTERPOLATE_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ psTrace(__func__, 4, "---- %s() begin %u-order.) (%d data points) ----\n", __func__, order, order+1); \ if (x->data.TYPE < domain->data.TYPE[0]) { \ @@ -213,16 +174,5 @@ psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \ return(out); \ -} \ - -FUNC_MACRO_FULL_INTERPOLATE(U8) -FUNC_MACRO_FULL_INTERPOLATE(U16) -FUNC_MACRO_FULL_INTERPOLATE(U32) -FUNC_MACRO_FULL_INTERPOLATE(U64) -FUNC_MACRO_FULL_INTERPOLATE(S8) -FUNC_MACRO_FULL_INTERPOLATE(S16) -FUNC_MACRO_FULL_INTERPOLATE(S32) -FUNC_MACRO_FULL_INTERPOLATE(S64) -FUNC_MACRO_FULL_INTERPOLATE(F32) -FUNC_MACRO_FULL_INTERPOLATE(F64) +} /***************************************************************************** @@ -236,8 +186,8 @@ psScalar *p_psVectorInterpolate( psScalar *out, - psVector *domain, - psVector *range, + const psVector *domain, + const psVector *range, psS32 order, - psScalar *x) + const psScalar *x) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -249,36 +199,27 @@ PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL); PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL); - if (out == NULL) { + if (!out) { out = psScalarAlloc(0, x->type.type); } else { PS_ASSERT_PTR_TYPE_EQUAL(domain, out, NULL); } - char* strType; switch (x->type.type) { - case PS_TYPE_U8: - return(vectorInterpolateU8(out, domain, range, order, x)); - case PS_TYPE_U16: - return(vectorInterpolateU16(out, domain, range, order, x)); - case PS_TYPE_U32: - return(vectorInterpolateU32(out, domain, range, order, x)); - case PS_TYPE_U64: - return(vectorInterpolateU64(out, domain, range, order, x)); - case PS_TYPE_S8: - return(vectorInterpolateS8(out, domain, range, order, x)); - case PS_TYPE_S16: - return(vectorInterpolateS16(out, domain, range, order, x)); - case PS_TYPE_S32: - return(vectorInterpolateS32(out, domain, range, order, x)); - case PS_TYPE_S64: - return(vectorInterpolateS64(out, domain, range, order, x)); - case PS_TYPE_F32: - return(vectorInterpolateF32(out, domain, range, order, x)); - case PS_TYPE_F64: - return(vectorInterpolateF64(out, domain, range, order, x)); - default: - PS_TYPE_NAME(strType, x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); - return NULL; + VECTOR_INTERPOLATE_CASE(U8); + VECTOR_INTERPOLATE_CASE(U16); + VECTOR_INTERPOLATE_CASE(U32); + VECTOR_INTERPOLATE_CASE(U64); + VECTOR_INTERPOLATE_CASE(S8); + VECTOR_INTERPOLATE_CASE(S16); + VECTOR_INTERPOLATE_CASE(S32); + VECTOR_INTERPOLATE_CASE(S64); + VECTOR_INTERPOLATE_CASE(F32); + VECTOR_INTERPOLATE_CASE(F64); + default: { + char* strType; + PS_TYPE_NAME(strType, x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); + return NULL; + } } @@ -287,34 +228,14 @@ /***************************************************************************** -These macros and functions define the following functions: - - p_psNormalizeVectorRange(myData, low, high) - -That assumes that the low/high arguments are PS_TYPE_F64; the vector myData -can be of any type. Arguments low/high will be converted to the appropriate -type and one of the type-specific functions below will be called: - - p_psNormalizeVectorRangeU8(myData, low, high) - p_psNormalizeVectorRangeU16(myData, low, high) - p_psNormalizeVectorRangeU32(myData, low, high) - p_psNormalizeVectorRangeU64(myData, low, high) - p_psNormalizeVectorRangeS8(myData, low, high) - p_psNormalizeVectorRangeS16(myData, low, high) - p_psNormalizeVectorRangeS32(myData, low, high) - p_psNormalizeVectorRangeS64(myData, low, high) - p_psNormalizeVectorRangeF32(myData, low, high) - p_psNormalizeVectorRangeF64(myData, low, high) +Helper macro for p_psNormalizeVectorRange: handles each of the types *****************************************************************************/ -#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \ -static void p_psNormalizeVectorRange##TYPE( \ - psVector* myData, \ - ps##TYPE outLow, \ - ps##TYPE outHigh) \ -{ \ - ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \ - ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \ - psS32 i = 0; \ - \ - for (i = 0; i < myData->n; i++) { \ +#define NORMALIZE_VECTOR_RANGE_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + ps##TYPE low = outLow; \ + ps##TYPE high = outHigh; \ + ps##TYPE min = (ps##TYPE)PS_MAX_##TYPE; \ + ps##TYPE max = (ps##TYPE)-PS_MAX_##TYPE; \ + \ + for (long i = 0; i < myData->n; i++) { \ if (myData->data.TYPE[i] < min) { \ min = myData->data.TYPE[i]; \ @@ -327,79 +248,53 @@ /* Ensure that max!=min before we divide by (max-min) */ \ if (max != min) { \ - for (i = 0; i < myData->n; i++) { \ - myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \ - (outHigh - outLow) / (max - min)); \ + for (long i = 0; i < myData->n; i++) { \ + myData->data.TYPE[i] = (low + (myData->data.TYPE[i] - min) * \ + (high - low) / (max - min)); \ } \ } else { \ psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min). Setting all elements to min.\n"); \ - for (i = 0; i < myData->n; i++) \ + for (long i = 0; i < myData->n; i++) \ { \ \ - myData->data.TYPE[i] = outLow; \ + myData->data.TYPE[i] = low; \ \ } \ } \ + break; \ } \ -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64) - -psS32 p_psNormalizeVectorRange(psVector* myData, - psF64 outLow, - psF64 outHigh) +/************************************************************************************************************* +p_psNormalizeVectorRange(myData, low, high): this function normalises the vector (myData) to the range +low: +high. +*************************************************************************************************************/ +bool p_psNormalizeVectorRange(psVector* myData, + psF64 outLow, + psF64 outHigh) { - PS_ASSERT_VECTOR_NON_NULL(myData, -1); - char* strType; - + PS_ASSERT_VECTOR_NON_NULL(myData, false); psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + switch (myData->type.type) { - case PS_TYPE_U8: - p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh); - break; - case PS_TYPE_U16: - p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh); - break; - case PS_TYPE_U32: - p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh); - break; - case PS_TYPE_U64: - p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh); - break; - case PS_TYPE_S8: - p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh); - break; - case PS_TYPE_S16: - p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh); - break; - case PS_TYPE_S32: - p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh); - break; - case PS_TYPE_S64: - p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh); - break; - case PS_TYPE_F32: - p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh); - break; - case PS_TYPE_F64: - p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh); - break; - case PS_TYPE_C32: - case PS_TYPE_C64: - default: - PS_TYPE_NAME(strType, myData->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); - break; + NORMALIZE_VECTOR_RANGE_CASE(U8); + NORMALIZE_VECTOR_RANGE_CASE(U16); + NORMALIZE_VECTOR_RANGE_CASE(U32); + NORMALIZE_VECTOR_RANGE_CASE(U64); + NORMALIZE_VECTOR_RANGE_CASE(S8); + NORMALIZE_VECTOR_RANGE_CASE(S16); + NORMALIZE_VECTOR_RANGE_CASE(S32); + NORMALIZE_VECTOR_RANGE_CASE(S64); + NORMALIZE_VECTOR_RANGE_CASE(F32); + NORMALIZE_VECTOR_RANGE_CASE(F64); + default: { + char* strType; + PS_TYPE_NAME(strType, myData->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); + break; + } } psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(0); -} - - + return true; +} + + Index: /trunk/psLib/src/math/psMathUtils.h =================================================================== --- /trunk/psLib/src/math/psMathUtils.h (revision 7131) +++ /trunk/psLib/src/math/psMathUtils.h (revision 7132) @@ -7,6 +7,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-17 00:56:48 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-18 01:20:51 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -34,6 +34,6 @@ */ psS32 p_psVectorBinDisect( - psVector *bins, ///< Array of non-decreasing values - psScalar *x ///< Target value to find + const psVector *bins, ///< Array of non-decreasing values + const psScalar *x ///< Target value to find ); @@ -44,14 +44,14 @@ */ psScalar *p_psVectorInterpolate( - psScalar *out, - psVector *domain, ///< Domain (x coords) for interpolation - psVector *range, ///< Range (y coords) for interpolation - psS32 order, ///< Order of interpolation function - psScalar *x ///< Location at which to evaluate + psScalar *out, ///< Output scalar, or NULL + const psVector *domain, ///< Domain (x coords) for interpolation + const psVector *range, ///< Range (y coords) for interpolation + psS32 order, ///< Order of interpolation function + const psScalar *x ///< Location at which to evaluate ); -psS32 p_psNormalizeVectorRange(psVector* myData, - psF64 outLow, - psF64 outHigh); +bool p_psNormalizeVectorRange(psVector* myData, + psF64 outLow, + psF64 outHigh); /** \} */ // End of MathGroup Functions Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 7131) +++ /trunk/psLib/src/math/psStats.c (revision 7132) @@ -16,6 +16,6 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.171 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-05-03 23:16:16 $ + * @version $Revision: 1.172 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-18 01:20:50 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -28,4 +28,5 @@ #include #include +#include /*****************************************************************************/ @@ -50,5 +51,5 @@ /* DEFINE STATEMENTS */ /*****************************************************************************/ -#define PS_GAUSS_WIDTH 5 // The width of the Gaussian smoothing. +#define PS_GAUSS_WIDTH 3 // The width of the Gaussian smoothing. // This corresponds to N in the ADD. #define PS_CLIPPED_NUM_ITER_LB 1 @@ -63,5 +64,7 @@ /* TYPE DEFINITIONS */ /*****************************************************************************/ -psVector* p_psConvertToF32(psVector* in); + +// None + /*****************************************************************************/ /* GLOBAL VARIABLES */ @@ -154,26 +157,25 @@ *****************************************************************************/ /****************************************************************************** -p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the +vectorSampleMean(myVector, maskVector, maskVal, stats): calculates the mean of the input vector. If there was a problem with the mean calculation, this routine sets stats->sampleMean to NAN. *****************************************************************************/ -psS32 p_psVectorSampleMean(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) +static bool vectorSampleMean(const psVector* myVector, + const psVector* errors, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psS32 i = 0; // Loop index variable - psF32 mean = 0.0; // The mean - psS32 count = 0; // # of points in this mean - - // If PS_STAT_USE_RANGE is requested, then we enter a slightly different - // loop. - if (errors == NULL) { + psF32 mean = 0.0; // The mean + long count = 0; // Number of points contributing to this mean + + // If PS_STAT_USE_RANGE is requested, then we enter a slightly different loop. + if (!errors) { if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { + if (maskVector) { + // No errors, check range, check mask + for (long i = 0; i < myVector->n; i++) { // Check if the data is with the specified range if (!(maskVal & maskVector->data.U8[i]) && @@ -184,12 +186,8 @@ } } - if (count != 0) { - mean /= (psF32)count; - } else { - mean = NAN; - } } else { - for (i = 0; i < myVector->n; i++) { + // No errors, check range, no mask + for (long i = 0; i < myVector->n; i++) { if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { @@ -198,13 +196,9 @@ } } - if (count != 0) { - mean /= (psF32)count; - } else { - mean = NAN; - } - } - } else { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { + } + } else { + if (maskVector) { + // No errors, no range check, check mask + for (long i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i])) { mean += myVector->data.F32[i]; @@ -212,77 +206,70 @@ } } - if (count != 0) { - mean /= (psF32)count; - } else { - mean = NAN; - } + } else { - for (i = 0; i < myVector->n; i++) { + // No errors, no range check, no mask + for (long i = 0; i < myVector->n; i++) { mean += myVector->data.F32[i]; } - mean /= (psF32)myVector->n; - } - } + count = myVector->n; + } + } + + if (count > 0) { + mean /= (psF32)count; + } else { + mean = NAN; + } + } else { - psF32 errorDivisor = 0.0; - psF32 errorSqr = 0.0; + psF32 sumWeights = 0.0; // The sum of the weights if (stats->options & PS_STAT_USE_RANGE) { if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { + for (long i = 0; i < myVector->n; i++) { // Check if the data is with the specified range if (!(maskVal & maskVector->data.U8[i]) && (stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { - errorSqr = errors->data.F32[i]*errors->data.F32[i]; - mean += myVector->data.F32[i]/errorSqr; - errorDivisor+= (1.0 / errorSqr); + float weight = 1.0 / PS_SQR(errors->data.F32[i]); + mean += myVector->data.F32[i] * weight; + sumWeights += weight; count++; } } - if (count != 0) { - mean /= errorDivisor; - } else { - mean = NAN; - } - } else { - for (i = 0; i < myVector->n; i++) { + for (long i = 0; i < myVector->n; i++) { if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { - errorSqr = errors->data.F32[i]*errors->data.F32[i]; - mean += myVector->data.F32[i]/errorSqr; - errorDivisor+= (1.0 / errorSqr); + float weight = 1.0 / PS_SQR(errors->data.F32[i]); + mean += myVector->data.F32[i] * weight; + sumWeights += weight; count++; } } - if (count != 0) { - mean /= errorDivisor; - } else { - mean = NAN; - } } } else { if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { + for (long i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i])) { - errorSqr = errors->data.F32[i]*errors->data.F32[i]; - mean += myVector->data.F32[i]/errorSqr; - errorDivisor+= (1.0 / errorSqr); + float weight = 1.0 / PS_SQR(errors->data.F32[i]); + mean += myVector->data.F32[i] * weight; + sumWeights += weight; count++; } } - if (count != 0) { - mean /= errorDivisor; - } else { - mean = NAN; - } } else { - for (i = 0; i < myVector->n; i++) { - errorSqr = errors->data.F32[i]*errors->data.F32[i]; - mean += myVector->data.F32[i]/errorSqr; - errorDivisor+= (1.0 / errorSqr); - } - mean /= errorDivisor; - } + for (long i = 0; i < myVector->n; i++) { + float weight = 1.0 / PS_SQR(errors->data.F32[i]); + mean += myVector->data.F32[i] * weight; + sumWeights += weight; + } + count = myVector->n; + } + } + + if (count > 0) { + mean /= sumWeights; + } else { + mean = NAN; } } @@ -290,182 +277,135 @@ stats->sampleMean = mean; if (isnan(mean)) { - psTrace(__func__, 4, "---- %s(-1) end ----\n", __func__); - return(-1); - } else { - psTrace(__func__, 4, "---- %s(0) end ----\n", __func__); - return(0); - } - + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; + } + + psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); + return true; } /****************************************************************************** -p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the -max of the input vector. If there was a problem with the max calculation, -this routine sets stats->max to NAN. +vectorMinMax(myVector, maskVector, maskVal, stats): calculates the minimum and maximum of the input vector. +If there was a problem with the calculation, this routine sets stats->max and stats->min to NAN. Return the +number of valid values in the vector (those not masked or outside the specified range. *****************************************************************************/ -psS32 p_psVectorMax(const psVector* myVector, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) +static long vectorMinMax(const psVector* myVector, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats + ) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psS32 i = 0; // Loop index variable - psF32 max = -PS_MAX_F32; // The calculated maximum - psS32 empty = true; // Does this vector have valid elements? + psF32 max = -PS_MAX_F32; // The calculated maximum + psF32 min = PS_MAX_F32; // The calculated minimum + psF32 *vector = myVector->data.F32; // Dereference the vector + psU8 *mask = NULL; // Dereference the mask + if (maskVector) { + mask = maskVector->data.U8; + } + int num = myVector->n; // Number of values + int numValid = 0; // Number of valid values // If PS_STAT_USE_RANGE is requested, then we enter a different loop. if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if ((myVector->data.F32[i] > max) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - max = myVector->data.F32[i]; - empty = false; + if (mask) { + // Need to check range and mask + for (long i = 0; i < num; i++) { + if (!(maskVal & mask[i]) && vector[i] >= stats->min && vector[i] <= stats->max) { + numValid++; + if (vector[i] > max) { + max = vector[i]; } - } - } - } else { - for (i = 0; i < myVector->n; i++) { - if ((myVector->data.F32[i] > max) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - max = myVector->data.F32[i]; - empty = false; + if (vector[i] < min) { + min = vector[i]; + } + } + } + } else { + // Need to check range + for (long i = 0; i < num; i++) { + if (vector[i] >= stats->min && vector[i] <= stats->max) { + numValid++; + if (vector[i] > max) { + max = vector[i]; + } + if (vector[i] < min) { + min = vector[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]; - empty = false; + if (mask) { + // Need to check mask + for (long i = 0; i < num; i++) { + if (!(maskVal & mask[i])) { + numValid++; + if (vector[i] > max) { + max = vector[i]; } - } - } - } else { - for (i = 0; i < myVector->n; i++) { - if (myVector->data.F32[i] > max) { - max = myVector->data.F32[i]; - empty = false; - } - } - } - } - if (empty == false) { + if (vector[i] < min) { + min = vector[i]; + } + } + } + } else { + // No checks + for (long i = 0; i < num; i++) { + if (vector[i] > max) { + max = vector[i]; + } + if (vector[i] < min) { + min = vector[i]; + } + } + numValid = num; + } + } + + + if (numValid == 0) { + stats->max = NAN; + stats->min = NAN; + } else { stats->max = max; - } else { - stats->max = NAN; - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - psTrace(__func__, 4, "---- %s(0) end ----\n", __func__); - return(0); + stats->min = min; + } + psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, numValid); + return numValid; } +#if 0 /****************************************************************************** -p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the -minimum of the input vector. If there was a problem with the min calculation, -this routine sets stats->min to NAN. - *****************************************************************************/ -psS32 p_psVectorMin(const psVector* myVector, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psS32 i = 0; // Loop index variable - psF32 min = PS_MAX_F32; // The calculated maximum - psS32 empty = true; // Does this vector have valid elements? - - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if ((myVector->data.F32[i] < min) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - min = myVector->data.F32[i]; - empty = false; - } - } - } - } else { - for (i = 0; i < myVector->n; i++) { - if ((myVector->data.F32[i] < min) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - min = myVector->data.F32[i]; - empty = false; - } - } - } - } 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]; - empty = false; - } - } - } - } else { - for (i = 0; i < myVector->n; i++) { - if (myVector->data.F32[i] < min) { - min = myVector->data.F32[i]; - empty = false; - } - } - } - } - - if (empty == false) { - stats->min = min; - } else { - stats->min = NAN; - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - psTrace(__func__, 4, "---- %s(0) end ----\n", __func__); - return(0); -} - -/****************************************************************************** -p_psVectorNValues(myVector, maskVector, maskVal, stats): This routine returns +vectorCheckNonEmpty(myVector, maskVector, maskVal, stats): This routine returns "true" if the inputPsVector has 1 or more valid elements (a valid element is an unmasked element within the specified min/max range). Otherwise, return "false". *****************************************************************************/ -bool p_psVectorCheckNonEmpty( - const psVector* myVector, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) +static bool vectorCheckNonEmpty(const psVector* myVector, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_VECTOR_NON_NULL(myVector, -1); - PS_ASSERT_PTR_NON_NULL(stats, -1); - psS32 i = 0; // Loop index variable + PS_ASSERT_VECTOR_NON_NULL(myVector, false); + PS_ASSERT_PTR_NON_NULL(stats, false); if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { + if (maskVector) { + for (long i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i]) && (stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return(true); - } - } - } else { - for (i = 0; i < myVector->n; i++) { + return true; + } + } + } else { + for (long i = 0; i < myVector->n; i++) { if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return(true); + return true; } } @@ -473,8 +413,8 @@ } else { if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { + for (long i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i])) { psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return(true); + return true; } } @@ -482,5 +422,5 @@ if (myVector->n > 0) { psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return(true); + return true; } } @@ -489,85 +429,21 @@ return(false); } - +#endif /****************************************************************************** -p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the -number of non-masked pixels in the vector that fall within the min/max -range, if specified. +vectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the median +and quartiles of the input vector. Returns true on success (including if there +were no valid input vector elements). *****************************************************************************/ -psS32 p_psVectorNValues(const psVector* myVector, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) +static bool vectorSampleMedian(const psVector* myVector, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_VECTOR_NON_NULL(myVector, -1); - PS_ASSERT_PTR_NON_NULL(stats, -1); - psS32 i = 0; // Loop index variable - psS32 numData = 0; // The number of data points - - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - numData++; - } - } - } else { - for (i = 0; i < myVector->n; i++) { - if ((stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - numData++; - } - } - } - } else { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i])) { - numData++; - } - } - } else { - numData = myVector->n; - } - } - - psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, numData); - return (numData); -} - -/****************************************************************************** -p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the -median of the input vector. Returns true on success (including if there were -no valid input vector elements). - *****************************************************************************/ -bool p_psVectorSampleMedian(const psVector* myVector, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psVector* unsortedVector = NULL; // Temporary vector - psVector* sortedVector = NULL; // Temporary vector - psS32 i = 0; // Loop index variable - psS32 count = 0; - psS32 nValues = 0; - - // Determine how many data points fit inside this min/max range - // and are not masked, if the maskVector is not NULL. - nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); - - if (nValues <= 0) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n"); - return(true); - stats->sampleMedian = NAN; - } // Allocate temporary vectors for the data. - unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = nValues; + psVector* vector = psVectorAlloc(myVector->n, PS_TYPE_F32); // Vector containing the unmasked values + long count = 0; // Number of valid entries // Determine if we must only use data points within a min/max range. @@ -575,19 +451,18 @@ // 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++) { + for (long i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i]) && (stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { - unsortedVector->data.F32[count] = myVector->data.F32[i]; + vector->data.F32[count] = myVector->data.F32[i]; count++; } } } else { - for (i = 0; i < myVector->n; i++) { + for (long i = 0; i < myVector->n; i++) { if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { - unsortedVector->data.F32[count] = myVector->data.F32[i]; + vector->data.F32[count] = myVector->data.F32[i]; count++; } @@ -596,55 +471,63 @@ } else { // Store all non-masked data points into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { + if (maskVector) { + for (long i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count] = myVector->data.F32[i]; + vector->data.F32[count] = myVector->data.F32[i]; count++; } } } else { - for (i = 0; i < myVector->n; i++) { - unsortedVector->data.F32[i] = myVector->data.F32[i]; - } - } - } - // Sort the temporary vectors. - sortedVector = psVectorSort(sortedVector, unsortedVector); - if (sortedVector == NULL) { + vector = psVectorCopy(vector, myVector, PS_TYPE_F32); + count = myVector->n; + } + } + vector->n = count; + if (count == 0) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No valid data in input vector.\n"); + psFree(vector); + stats->sampleMedian = NAN; + return false; + } + + // Sort the temporary vector. + vector = psVectorSort(vector, vector); // Sort in-place (since it's a copy, it's OK) + if (!vector) { psError(PS_ERR_UNEXPECTED_NULL, false, PS_ERRORTEXT_psStats_STATS_SAMPLE_MEDIAN_SORT_PROBLEM); psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); - psFree(unsortedVector); - return(false); + psFree(vector); + return false; } // Calculate the median exactly. Use the average if the number of samples // is even. - if (0 == (nValues % 2)) { - stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues / 2) - 1] + - sortedVector->data.F32[nValues / 2]); + if (count % 2 == 0) { + stats->sampleMedian = 0.5 * (vector->data.F32[(count / 2) - 1] + + vector->data.F32[count / 2]); } else { - stats->sampleMedian = sortedVector->data.F32[nValues / 2]; - } + stats->sampleMedian = vector->data.F32[count / 2]; + } + + // Calculate the quartile points exactly. + stats->sampleUQ = vector->data.F32[3 * (count / 4)]; + stats->sampleLQ = vector->data.F32[count / 4]; // Free the temporary data structures. - psFree(unsortedVector); - psFree(sortedVector); + psFree(vector); // Return "true" on success. psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return(true); + return true; } /****************************************************************************** -p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input +vectorSmoothHistGaussian(): This routine smoothes the data in the input robustHistogram with a Gaussian of width sigma. It returns a psVector of the smoothed data. *****************************************************************************/ -psVector *p_psVectorSmoothHistGaussian( - psHistogram *histogram, - psF32 sigma) +psVector *vectorSmoothHistGaussian(psHistogram *histogram, + psF32 sigma) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -657,20 +540,19 @@ } - psS32 numBins = histogram->nums->n; - psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); + long numBins = histogram->nums->n; // Number of histogram bins + psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins smooth->n = smooth->nalloc; - psS32 jMin = 0; - psS32 jMax = 0; - - if (histogram->uniform == false) { + const psVector *bounds = histogram->bounds; // The bounds for the histogram bins + + if (!histogram->uniform) { // // We get here if the histogram is non-uniform. This code is not tested. // However, it is also not used anywhere, yet. // - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n"); - - for (psS32 i = 0; i < numBins; i++) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n"); + + for (long i = 0; i < numBins; i++) { // Determine the midpoint of bin i. - psS32 iMid = PS_BIN_MIDPOINT(histogram, i); + float iMid = PS_BIN_MIDPOINT(histogram, i); // @@ -678,24 +560,11 @@ // range of data values surrounding iMid. The range is of size: // 2*PS_GAUSS_WIDTH*sigma + long jMin, jMax; psF32 x = iMid - (PS_GAUSS_WIDTH * sigma); - jMin = i; - while (jMin > 0) { - if ((histogram->bounds->data.F32[jMin] <= x) && - (histogram->bounds->data.F32[jMin+1] >= x)) { - break; - } - jMin--; - } - + for (jMin = i; jMin > 0 && bounds->data.F32[jMin] > x; jMin--) + ; x = iMid + (PS_GAUSS_WIDTH * sigma); - jMax = i; - while (jMax < (histogram->bounds->n - 1)) { - if ((histogram->bounds->data.F32[jMax] <= x) && - (histogram->bounds->data.F32[jMax+1] >= x)) { - break; - } - jMax++; - } - + for (jMax = i; jMax < bounds->n - 1 && bounds->data.F32[jMax + 1] > x; jMax++) + ; // @@ -703,6 +572,6 @@ // smooth->data.F32[i] = 0.0; - for (psS32 j = jMin ; j <= jMax ; j++) { - psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + for (long j = jMin ; j <= jMax ; j++) { + float jMid = PS_BIN_MIDPOINT(histogram, j); smooth->data.F32[i] += histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); @@ -713,7 +582,7 @@ // We get here if the histogram is uniform. // - for (psS32 i = 0; i < numBins; i++) { - psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0]; - psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize); + for (long i = 0; i < numBins; i++) { + psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0]; + psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize); // @@ -723,5 +592,8 @@ // took way too long. // + #if 0 + gaussWidth = 10; + #endif // // We determine the bin numbers (jMin:jMax) corresponding to a @@ -729,12 +601,6 @@ // 2*PS_GAUSS_WIDTH*sigma // - psS32 jMin = i - gaussWidth; - if (jMin < 0 ) { - jMin = 0; - } - psS32 jMax = i + gaussWidth; - if (jMax > (histogram->bounds->n - 1)) { - jMax = (histogram->bounds->n - 1); - } + psS32 jMin = PS_MAX(i - gaussWidth, 0); + psS32 jMax = PS_MIN(i + gaussWidth, bounds->n - 1); // @@ -742,7 +608,7 @@ // smooth->data.F32[i] = 0.0; - psS32 iMid = PS_BIN_MIDPOINT(histogram, i); - for (psS32 j = jMin ; j <= jMax ; j++) { - psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + float iMid = PS_BIN_MIDPOINT(histogram, i); + for (long j = jMin ; j <= jMax ; j++) { + float jMid = PS_BIN_MIDPOINT(histogram, j); smooth->data.F32[i] += histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); @@ -757,92 +623,7 @@ return(smooth); } + /****************************************************************************** -p_psVectorSampleQuartiles(myVector, maskVector, maskVal, statsalculates -the upper and/or lower quartiles of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -bool p_psVectorSampleQuartiles(const psVector* myVector, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psVector* unsortedVector = NULL; - psVector* sortedVector = NULL; - psS32 i = 0; - psS32 count = 0; - psS32 nValues = 0; - - // Determine how many data points fit inside this min/max range - // and are not masked, if the maskVector is not NULL. - nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); - - // Allocate temporary vectors for the data. - unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = unsortedVector->nalloc; - - if (stats->options & PS_STAT_USE_RANGE) { - // 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]) && - (stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } else { - for (i = 0; i < myVector->n; i++) { - if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { - 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. - sortedVector = psVectorSort(sortedVector, unsortedVector); - if (sortedVector == NULL) { - psError(PS_ERR_UNEXPECTED_NULL, - false, - PS_ERRORTEXT_psStats_STATS_SAMPLE_MEDIAN_SORT_PROBLEM); - psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); - return(false); - } - - // 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); - psFree(sortedVector); - psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return(true); -} - -/****************************************************************************** -p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the +vectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the stdev of the input vector. Inputs @@ -855,119 +636,123 @@ *****************************************************************************/ -void p_psVectorSampleStdev(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) +static bool vectorSampleStdev(const psVector* myVector, + const psVector* errors, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psS32 i = 0; // Loop index variable - psS32 countInt = 0; // # of data points being used - psF32 countFloat = 0.0; // # of data points being used - psF32 mean = 0.0; // The mean - psF32 diff = 0.0; // Used in calculating stdev - psF32 sumSquares = 0.0; // temporary variable - psF32 sumDiffs = 0.0; // temporary variable - psF32 errorDivisor = 0.0f; // This procedure requires the mean. If it has not been already - // calculated, then call p_psVectorSampleMean() - if (0 != isnan(stats->sampleMean)) { - p_psVectorSampleMean(myVector, errors, maskVector, maskVal, stats); + // calculated, then call vectorSampleMean() + if (isnan(stats->sampleMean)) { + vectorSampleMean(myVector, errors, maskVector, maskVal, stats); } // If the mean is NAN, then generate a warning and set the stdev to NAN. - if (0 != isnan(stats->sampleMean)) { + if (isnan(stats->sampleMean)) { stats->sampleStdev = NAN; - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n"); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleStdev(): vectorSampleMean() reported a NAN mean.\n"); psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return; - } - - mean = stats->sampleMean; + return false; + } + + // Accumulate the sums + double mean = stats->sampleMean; // The mean + double sumSquares = 0.0; // Sum of the squares + double sumDiffs = 0.0; // Sum of the differences + double errorDivisor = 0.0; // Division for errors + long count = 0; // Number of data points being used if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { + if (maskVector) { + for (long i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i]) && (stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { - diff = myVector->data.F32[i] - mean; - sumSquares += (diff * diff); + double diff = myVector->data.F32[i] - mean; + sumSquares += PS_SQR(diff); sumDiffs += diff; - countInt++; - if (errors != NULL) { - errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i])); + count++; + if (errors) { + errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]); } } } } else { - for (i = 0; i < myVector->n; i++) { + for (long i = 0; i < myVector->n; i++) { if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { - diff = myVector->data.F32[i] - mean; - sumSquares += (diff * diff); + double diff = myVector->data.F32[i] - mean; + sumSquares += PS_SQR(diff); sumDiffs += diff; - countInt++; - if (errors != NULL) { - errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i])); + count++; + if (errors) { + errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]); } } } - countInt = myVector->n; + count = myVector->n; } } else { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { + if (maskVector) { + for (long i = 0; i < myVector->n; i++) { if (!(maskVal & maskVector->data.U8[i])) { - diff = myVector->data.F32[i] - mean; - sumSquares += (diff * diff); + double diff = myVector->data.F32[i] - mean; + sumSquares += PS_SQR(diff); sumDiffs += diff; - countInt++; - if (errors != NULL) { - errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i])); + count++; + if (errors) { + errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]); } } } } else { - for (i = 0; i < myVector->n; i++) { - diff = myVector->data.F32[i] - mean; - sumSquares += (diff * diff); + for (long i = 0; i < myVector->n; i++) { + double diff = myVector->data.F32[i] - mean; + sumSquares += PS_SQR(diff); sumDiffs += diff; - countInt++; - } - countInt = myVector->n; - if (errors != NULL) { - errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i])); - } - } - } - - if (countInt == 0) { + count++; + if (errors) { + errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]); + } + } + count = myVector->n; + } + } + + if (count == 0) { stats->sampleStdev = NAN; - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements (%d). Setting stats->sampleStdev = NAN.\n", countInt); - } else if (countInt == 1) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleStdev(): no valid psVector elements (%d). Setting stats->sampleStdev = NAN.\n", count); + return false; + } + if (count == 1) { stats->sampleStdev = 0.0; - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d). Setting stats->sampleStdev = 0.0.\n", countInt); - } else { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleStdev(): only one valid psVector elements (%d). Setting stats->sampleStdev = 0.0.\n", count); + return false; + } + + if (errors) { + #if 1 // XXX: The ADD specifies this as the definition of the standard // deviation if the errors are known. Verify this with IfA: none of // the data points in the vector are used. Verify that the masks and // data ranges are used correctly. - if (errors != NULL) { - if (0) { - stats->sampleStdev = (1.0 / sqrtf(errorDivisor)); - } else { - countFloat = (psF32)countInt; - stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); - } - } else { - countFloat = (psF32)countInt; - stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); - - } + stats->sampleStdev = (1.0 / sqrtf(errorDivisor)); + #else + + stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) / + (float)(count - 1)); + #endif + + } else { + stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) / + (float)(count - 1)); + } psTrace(__func__, 4, "---- %s() end ----\n", __func__); + + return true; } /****************************************************************************** -p_psVectorClippedStats(myVector, errors, maskVector, maskVal, stats): calculates the +vectorClippedStats(myVector, errors, maskVector, maskVal, stats): calculates the clipped stats (mean or stdev) of the input vector. @@ -978,21 +763,15 @@ stats Returns - 0 for success. - -1: error - -2: warning + true for success; false otherwise *****************************************************************************/ -psS32 p_psVectorClippedStats(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) +static bool vectorClippedStats(const psVector* myVector, + const psVector* errors, + psVector* maskVector, + psMaskType maskVal, + psStats* stats + ) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); psTrace(__func__, 4, "Trace level is %d\n", psTraceGetLevel(__func__)); - psF32 clippedMean = 0.0; - psF32 clippedStdev = 0.0; - psVector* tmpMask = NULL; - static psStats *statsTmp = NULL; - psS32 rc = 0; // Ensure that stats->clipIter is within the proper range. @@ -1008,94 +787,69 @@ // Allocate a psStats structure for calculating the mean, median, and // stdev. - if (statsTmp == NULL) { - statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - p_psMemSetPersistent(statsTmp, true); + psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + + // We've already copied the mask vector, so it's safe to simply increment the reference counter. + psVector *tmpMask = NULL; + if (maskVector) { + tmpMask = psMemIncrRefCounter(maskVector); } else { - // Initialize structure if already allocated - statsTmp->sampleMean = NAN; - statsTmp->sampleMedian = NAN; - statsTmp->sampleStdev = NAN; - statsTmp->sampleUQ = NAN; - statsTmp->sampleLQ = NAN; - statsTmp->robustMedian = NAN; - statsTmp->robustStdev = NAN; - statsTmp->robustUQ = NAN; - statsTmp->robustLQ = NAN; - statsTmp->robustN50 = -1; - statsTmp->fittedMean = NAN; - statsTmp->fittedStdev = NAN; - statsTmp->fittedNfit = -1; - statsTmp->clippedMean = NAN; - statsTmp->clippedStdev = NAN; - statsTmp->clippedNvalues = -1; - statsTmp->clipSigma = 3.0; - statsTmp->clipIter = 3; - statsTmp->min = NAN; - statsTmp->max = NAN; - statsTmp->binsize = NAN; - statsTmp->nSubsample = 100000; - } - - // 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->n = tmpMask->nalloc; - - // If we were called with a mask vector, then initialize the temporary - // mask vector with those values. Otherwise, initialize to zero. - if (maskVector != NULL) { - for (psS32 i = 0; i < tmpMask->n; i++) { - tmpMask->data.U8[i] = maskVector->data.U8[i]; - } - } else { - for (psS32 i = 0; i < tmpMask->n; i++) { - tmpMask->data.U8[i] = 0; - } + tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8); + tmpMask->n = tmpMask->nalloc; + psVectorInit(tmpMask, 0); } // 1. Compute the sample median. - p_psVectorSampleMedian(myVector, maskVector, maskVal, statsTmp); + vectorSampleMedian(myVector, tmpMask, maskVal, statsTmp); if (isnan(statsTmp->sampleMedian)) { - psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleMedian returned NAN\n"); + psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleMedian returned NAN\n"); stats->clippedMean = NAN; stats->clippedStdev = NAN; - psTrace(__func__, 4, "---- %s(-2) end ----\n", __func__); - return(-2); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + psFree(tmpMask); + psFree(statsTmp); + return false; } psTrace(__func__, 6, "The initial sample median is %f\n", statsTmp->sampleMedian); // 2. Compute the sample standard deviation. - p_psVectorSampleStdev(myVector, errors, maskVector, maskVal, statsTmp); + vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp); if (isnan(statsTmp->sampleStdev)) { - psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleStdev returned NAN\n"); + psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleStdev returned NAN\n"); stats->clippedMean = NAN; stats->clippedStdev = NAN; - psTrace(__func__, 4, "---- %s(-2) end ----\n", __func__); - return(-2); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + psFree(tmpMask); + psFree(statsTmp); + return false; } psTrace(__func__, 6, "The initial sample stdev is %f\n", statsTmp->sampleStdev); // 3. Use the sample median as the first estimator of the mean X. - clippedMean = statsTmp->sampleMedian; + psF32 clippedMean = statsTmp->sampleMedian; // 4. Use the sample stdev as the first estimator of the mean stdev. - clippedStdev = statsTmp->sampleStdev; + psF32 clippedStdev = statsTmp->sampleStdev; // 5. Repeat N (stats->clipIter) times: - for (psS32 iter = 0; iter < stats->clipIter; iter++) { + long numClipped = LONG_MAX; // Number of values clipped in this iteration + for (int iter = 0; iter < stats->clipIter && numClipped > 0; iter++) { + numClipped = 0; psTrace(__func__, 6, "------------ Iteration %d ------------\n", iter); // a) Exclude all values x_i for which |x_i - x| > K * stdev - if (errors != NULL) { - for (psS32 j = 0; j < myVector->n; j++) { - if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * errors->data.F32[j])) { + if (errors) { + for (long j = 0; j < myVector->n; j++) { + if (fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * errors->data.F32[j])) { tmpMask->data.U8[j] = 0xff; - } - } - } else { - for (psS32 j = 0; j < myVector->n; j++) { - if (fabs(myVector->data.F32[j] - clippedMean) > - (stats->clipSigma * clippedStdev)) { + psTrace(__func__, 10, "Clipped %d: %f +/- %f\n", j, + myVector->data.F32[j], errors->data.F32[j]); + numClipped++; + } + } + } else { + for (long j = 0; j < myVector->n; j++) { + if (fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { tmpMask->data.U8[j] = 0xff; + psTrace(__func__, 10, "Clipped %d: %f\n", j, myVector->data.F32[j]); + numClipped++; } } @@ -1103,6 +857,6 @@ // b) compute new mean and stdev - p_psVectorSampleMean(myVector, errors, tmpMask, 0xff, statsTmp); - p_psVectorSampleStdev(myVector, errors, tmpMask, 0xff, statsTmp); + vectorSampleMean(myVector, errors, tmpMask, maskVal, statsTmp); + vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp); psTrace(__func__, 6, "The new sample mean is %f\n", statsTmp->sampleMean); psTrace(__func__, 6, "The new sample stdev is %f\n", statsTmp->sampleStdev); @@ -1113,7 +867,7 @@ // Exit loop. iter = stats->clipIter; - psError(PS_ERR_UNKNOWN, true, "p_psVectorSampleMean() or p_psVectorSampleStdev() returned a NAN.\n"); + psError(PS_ERR_UNKNOWN, true, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n"); psFree(tmpMask); - return(-1); + return false; } else { clippedMean = statsTmp->sampleMean; @@ -1122,10 +876,10 @@ } - // 7. The last calcuated value of x is the cliped mean. + // 7. The last calcuated value of x is the clipped mean. if (stats->options & PS_STAT_CLIPPED_MEAN) { stats->clippedMean = clippedMean; psTrace(__func__, 6, "The final clipped mean is %f\n", clippedMean); } - // 8. The last calcuated value of stdev is the cliped stdev. + // 8. The last calcuated value of stdev is the clipped stdev. if (stats->options & PS_STAT_CLIPPED_STDEV) { stats->clippedStdev = clippedStdev; @@ -1134,15 +888,16 @@ psFree(tmpMask); - psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, rc); - return(rc); + psFree(statsTmp); + psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); + return true; } -static psF32 QuadraticInverse( - psF32 a, - psF32 b, - psF32 c, - psF32 y, - psF32 xLo, - psF32 xHi) +static psF32 QuadraticInverse(psF32 a, + psF32 b, + psF32 c, + psF32 y, + psF32 xLo, + psF32 xHi + ) { psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a)); @@ -1151,20 +906,21 @@ psF64 x2 = -b/(2.0*a) - tmp; - if (0) { - psF64 y1 = (a * x1 * x1) + (b * x1) + c; - psF64 y2 = (a * x2 * x2) + (b * x2) + c; - printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c); - printf("QuadraticInverse: y is %f\n", y); - printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2); - printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2); - } - - if ((xLo <= x1) && (x1 <= xHi)) { - return(x1); - } else if ((xLo <= x2) && (x2 <= xHi)) { - return(x2); - } else { - return(0.5 * (xLo + xHi)); - } + #if 0 + + psF64 y1 = (a * x1 * x1) + (b * x1) + c; + psF64 y2 = (a * x2 * x2) + (b * x2) + c; + printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c); + printf("QuadraticInverse: y is %f\n", y); + printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2); + printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2); + #endif + + if (xLo <= x1 && x1 <= xHi) { + return x1; + } + if (xLo <= x2 && x2 <= xHi) { + return x2; + } + return 0.5 * (xLo + xHi); } @@ -1180,9 +936,9 @@ *****************************************************************************/ -psF32 fitQuadraticSearchForYThenReturnX( - psVector *xVec, - psVector *yVec, - psS32 binNum, - psF32 yVal) +static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, + psVector *yVec, + psS32 binNum, + psF32 yVal + ) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -1235,9 +991,10 @@ if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) || ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) { - psError(PS_ERR_UNKNOWN, true, "This routine must be called with monotonically increasing or decreasing data points.\n"); + psError(PS_ERR_UNKNOWN, true, + "This routine must be called with monotonically increasing or decreasing data points.\n"); psFree(x); psFree(y); psTrace(__func__, 5, "---- %s() end ----\n", __func__); - return(NAN); + return NAN; } @@ -1251,5 +1008,5 @@ psFree(y); psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__); - return(NAN); + return NAN; } psTrace(__func__, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]); @@ -1262,5 +1019,6 @@ psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); - tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]); + tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, + x->data.F64[0], x->data.F64[2]); psFree(myPoly); @@ -1295,5 +1053,5 @@ psTrace(__func__, 5, "---- %s(%f) end ----\n", __func__, tmpFloat); - return(tmpFloat); + return tmpFloat; } @@ -1301,8 +1059,8 @@ NOTE: We assume unnormalized gaussians. *****************************************************************************/ -psF32 psMinimizeLMChi2Gauss1D( - psVector *deriv, - const psVector *params, - const psVector *coords) +static psF32 minimizeLMChi2Gauss1D(psVector *deriv, + const psVector *params, + const psVector *coords + ) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -1318,23 +1076,20 @@ psF32 stdev = params->data.F32[1]; - if (deriv == NULL) { - deriv = psVectorAlloc(2, PS_TYPE_F32); - deriv->n = 2; - } else { + psF32 gauss = psGaussian(x, mean, stdev, false); + if (deriv) { PS_ASSERT_VECTOR_SIZE(deriv, 2, NAN); PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN); - } - - psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false); - deriv->data.F32[0] = tmp / (stdev * stdev); - tmp = (x - mean) * (x - mean) * psGaussian(x, mean, stdev, 0); - deriv->data.F32[1] = tmp / (stdev * stdev * stdev); + psF32 tmp = (x - mean) * gauss; + deriv->data.F32[0] = tmp / (stdev * stdev); + deriv->data.F32[1] = (x - mean) * tmp / (stdev * stdev * stdev); + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(psGaussian(x, mean, stdev, false)); + return gauss; } /****************************************************************************** -p_psVectorRobustStats(myVector, maskVector, maskVal, stats): This is the new +vectorRobustStats(myVector, maskVector, maskVal, stats): This is the new version of the robust stats routine. @@ -1354,9 +1109,9 @@ *****************************************************************************/ #define INITIAL_NUM_BINS 1000.0 -psS32 p_psVectorRobustStats(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) +static bool vectorRobustStats(const psVector* myVector, + const psVector* errors, + psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -1364,30 +1119,64 @@ PS_VECTOR_PRINT_F32(myVector); } - psHistogram *robustHistogram = NULL; - psHistogram *cumulativeRobustHistogram = NULL; - psS32 numBins = 0; - psScalar tmpScalar; + + psScalar tmpScalar; // Temporary scalar variable, for p_psVectorBinDisect tmpScalar.type.type = PS_TYPE_F32; - psS32 totalDataPoints = 0; - psS32 rc = 0; - psS32 rcBool = false; - psVector *tmpMaskVec = psVectorAlloc(myVector->n, PS_TYPE_U8); - tmpMaskVec->n = tmpMaskVec->nalloc; - if (maskVector != NULL) { - for (psS32 i = 0 ; i < myVector->n ; i++) { - tmpMaskVec->data.U8[i] = maskVector->data.U8[i]; - } + + // We need a mask for these statistics, so if the one we're given is NULL, make our own. + psVector *mask = NULL; // The actual mask we will use + if (maskVector) { + mask = psMemIncrRefCounter(maskVector); } else { - for (psS32 i = 0 ; i < myVector->n ; i++) { - tmpMaskVec->data.U8[i] = 0; - } - } - psS32 iterate = 1; - psF32 sigma; - psStats *tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); - - while (iterate > 0) { - psTrace(__func__, 6, "-------------------- Iterating on Bin size. Iteration number %d --------------------\n", iterate); - psF32 binSize = 0.0; + mask = psVectorAlloc(myVector->n, PS_TYPE_U8); + mask->n = myVector->n; + psVectorInit(mask, 0); + } + + psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max + psHistogram *histogram = NULL; // Histogram of the data + psHistogram *cumulative = NULL; // Cumulative histogram of the data + float min = NAN, max = NAN; // Mimimum and maximum values + float sigma = NAN; // The robust standard deviation + long totalDataPoints = 0; // Total number of (unmasked) data points + + // Iterate to get the best bin size + for (int iterate = 1; iterate > 0; iterate++) { + psTrace(__func__, 6, + "-------------------- Iterating on Bin size. Iteration number %d --------------------\n", + iterate); + + // Get the minimum and maximum values + int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of valid values + min = statsMinMax->min; + max = statsMinMax->max; + if (numValid == 0 || isnan(min) || isnan(max)) { + // Data range calculation failed + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(statsMinMax); + psFree(mask); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; + } + psTrace(__func__, 6, "Data min/max is (%.2f, %.2f)\n", min, max); + + // If all data points have the same value, then we set the appropiate members of stats and return. + if (fabs(max - min) <= FLT_EPSILON) { + psTrace(__func__, 5, "All data points have the same value: %f.\n", min); + if (stats->options & PS_STAT_ROBUST_MEDIAN) { + stats->robustMedian = min; + } + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + stats->robustUQ = min; + stats->robustLQ = min; + } + stats->robustN50 = numValid; + psFree(statsMinMax); + psFree(mask); + + psTrace(__func__, 4, "---- %s(0) end ----\n", __func__); + return false; + } + + float binSize = 0.0; // Size of bins for histogram if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) { // Set initial bin size to the specified value. @@ -1395,85 +1184,39 @@ psTrace(__func__, 6, "Setting initial robust bin size to %.2f\n", binSize); } else { - // Determine the bin size of the robust histogram. This is done - // by computing the total range of data values and dividing by 1000.0. - - rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax); - rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax); - if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - psFree(tmpStatsMinMax); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - psTrace(__func__, 6, "Data min/max is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max); - binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / INITIAL_NUM_BINS; + // Determine the bin size of the robust histogram, using the pre-defined number of bins + binSize = (max - min) / INITIAL_NUM_BINS; } psTrace(__func__, 6, "Initial robust bin size is %.2f\n", binSize); - // - // If all data points have the same value, then we set the appropiate - // members of stats and return. - // - if (fabs(tmpStatsMinMax->max - tmpStatsMinMax->min) <= FLT_EPSILON) { - psTrace(__func__, 5, "All data points have the same value.\n", binSize); - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - stats->robustMedian = tmpStatsMinMax->min; - } - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - stats->robustUQ = tmpStatsMinMax->min; - stats->robustLQ = tmpStatsMinMax->min; - } - // XXX: Set these to the number of unmasked data points? - stats->robustN50 = 0; - psFree(tmpStatsMinMax); - psFree(tmpMaskVec); - - psTrace(__func__, 4, "---- %s(0) end ----\n", __func__); - return(0); - } - - // - // ADD: Step 0. - // Construct the histogram with the specified bin size. - // - // NOTE: we can not specify the bin size precisely since the argument - // to psHistogramAlloc() is the number of bins, not the binSize. - // If we get here, we know that binSize != 0.0. - // - numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize); + // ADD step 0: Construct the histogram with the specified bin size. NOTE: we can not specify the bin + // size precisely since the argument to psHistogramAlloc() is the number of bins, not the binSize. If + // we get here, we know that binSize != 0.0. + long numBins = (max - min) / binSize; // Number of bins psTrace(__func__, 6, "Numbins is %d\n", numBins); - psTrace(__func__, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max); - robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); - cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); - - // Populate the histogram array. - robustHistogram = psVectorHistogram(robustHistogram, myVector, errors, tmpMaskVec, maskVal); + psTrace(__func__, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max); + // Generate the histogram + histogram = psHistogramAlloc(min, max, numBins); + histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal); if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(robustHistogram->bounds); - PS_VECTOR_PRINT_F32(robustHistogram->nums); - } - // - // ADD: Step 1. - // Construct the cumulative histogram from the specific histogram - // - cumulativeRobustHistogram->nums->data.F32[0] = robustHistogram->nums->data.F32[0]; - for (psS32 i = 1 ; i < robustHistogram->nums->n ; i++) { - cumulativeRobustHistogram->nums->data.F32[i] = cumulativeRobustHistogram->nums->data.F32[i-1] + - robustHistogram->nums->data.F32[i]; + PS_VECTOR_PRINT_F32(histogram->bounds); + PS_VECTOR_PRINT_F32(histogram->nums); + } + + // ADD step 1: Convert the specific histogram to a cumulative histogram + cumulative = psHistogramAlloc(min, max, numBins); + cumulative->nums->data.F32[0] = histogram->nums->data.F32[0]; + for (long i = 1; i < histogram->nums->n; i++) { + cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i]; } if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(cumulativeRobustHistogram->bounds); - PS_VECTOR_PRINT_F32(cumulativeRobustHistogram->nums); - } - - // - // ADD: Step 2. - // Find the bin which contains the 50% data point. - // - totalDataPoints = cumulativeRobustHistogram->nums->data.F32[numBins - 1]; + PS_VECTOR_PRINT_F32(cumulative->bounds); + PS_VECTOR_PRINT_F32(cumulative->nums); + } + + // ADD step 2: Find the bin which contains the 50% data point. + totalDataPoints = cumulative->nums->data.F32[numBins - 1]; psTrace(__func__, 6, "Total data points is %d\n", totalDataPoints); - psS32 binMedian; - if (totalDataPoints/2.0 < cumulativeRobustHistogram->nums->data.F32[0]) { + long binMedian; + if (totalDataPoints/2.0 < cumulative->nums->data.F32[0]) { // Special case: the median is in the first bin. Perhaps we should recode this. // XXX: Must try a special case where the median in in the last bin. @@ -1481,257 +1224,214 @@ } else { tmpScalar.data.F32 = totalDataPoints/2.0; - binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); + binMedian = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); if (binMedian < 0) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - } - psTrace(__func__, 6, "The median bin is %d (%.2f to %.2f)\n", binMedian, cumulativeRobustHistogram->bounds->data.F32[binMedian], cumulativeRobustHistogram->bounds->data.F32[binMedian+1]); - - // - // ADD: Step 3. - // Interpolate to the exact 50% position: this is the robust histogram median. - // - stats->robustMedian = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&cumulativeRobustHistogram->bounds, - *(psVector* *)&cumulativeRobustHistogram->nums, - binMedian, - totalDataPoints/2.0); + psError(PS_ERR_UNKNOWN, false, + "Failed to calculate the 50 precent data point (%d).\n", binMedian); + psFree(statsMinMax); + psFree(histogram); + psFree(cumulative); + psFree(mask); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; + } + } + psTrace(__func__, 6, "The median bin is %d (%.2f to %.2f)\n", binMedian, + cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]); + + // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median. + stats->robustMedian = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, + binMedian, totalDataPoints/2.0); if (isnan(stats->robustMedian)) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a quadratic and calculate the 50-percent position.\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); + psError(PS_ERR_UNKNOWN, false, + "Failed to fit a quadratic and calculate the 50-percent position.\n"); + psFree(statsMinMax); + psFree(histogram); + psFree(cumulative); + psFree(mask); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; } psTrace(__func__, 6, "Current robust median is %f\n", stats->robustMedian); - // - // ADD: Step 4. - // Find the bins which contains the 15.8655% and 84.1345% data points. - // - psS32 binLo = 0; - psS32 binHi = 0; + // ADD step 4: Find the bins which contains the 15.8655% and 84.1345% data points. + long binLo = 0; // Bin containing the 15.8655% mark tmpScalar.data.F32 = totalDataPoints * 0.158655f; - if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { - binLo = 0; - } else { + if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) { // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined. - binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); - if (binLo > cumulativeRobustHistogram->nums->n-1) { - binLo = cumulativeRobustHistogram->nums->n-1; - } - } + binLo = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); + if (binLo > cumulative->nums->n-1) { + binLo = cumulative->nums->n-1; + } + } + + long binHi = 0; // Bin containing the 84.1345% mark tmpScalar.data.F32 = totalDataPoints * 0.841345f; - if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { - binHi = 0; - } else { + if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) { // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined. - binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); - if (binHi > cumulativeRobustHistogram->nums->n-1) { - binHi = cumulativeRobustHistogram->nums->n-1; - } - } - psTrace(__func__, 6, "The 15.8655-percent and 84.1345-percent data point bins are (%d, %d).\n", + binHi = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); + if (binHi > cumulative->nums->n-1) { + binHi = cumulative->nums->n-1; + } + } + psTrace(__func__, 6, "The 15.8655%% and 84.1345%% data point bins are (%d, %d).\n", binLo, binHi); - psTrace(__func__, 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binLo)); - psTrace(__func__, 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binHi)); + psTrace(__func__, 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulative, binLo)); + psTrace(__func__, 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulative, binHi)); if ((binLo < 0) || (binHi < 0)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655 and 84.1345 percent data point\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - - // - // ADD: Step 4b. - // Interpolate Sigma to find these two positions exactly: these are the 1sigma positions. - // - // XXX: I am overriding the ADD for now. The following code implements - // the ADD exactly and is having problems fitting a 2nd-order polynomial - // to data y-values suchs as (1, 1, 100). Therefore, I commented out the - // code and am doing simply linear interpolation. - // - psF32 binLoF32; - psF32 binHiF32; - if (0) { - binLoF32 = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&cumulativeRobustHistogram->bounds, - *(psVector* *)&cumulativeRobustHistogram->nums, - binLo, - totalDataPoints * 0.158655f); - binHiF32 = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&cumulativeRobustHistogram->bounds, - *(psVector* *)&cumulativeRobustHistogram->nums, - binHi, - totalDataPoints * 0.841345f); - psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", - binLoF32, binHiF32); - } - - // + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n"); + psFree(statsMinMax); + psFree(histogram); + psFree(cumulative); + psFree(mask); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; + } + + // ADD step 4b: Interpolate Sigma to find these two positions exactly: these are the 1sigma positions. + #if 0 + // XXX: I am overriding the ADD for now. The following code implements the ADD exactly and is having + // problems fitting a 2nd-order polynomial to data y-values suchs as (1, 1, 100). Therefore, I + // commented out the code and am doing simply linear interpolation. + + // Value for the 15.8655% mark + float binLoF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binLo, + totalDataPoints * 0.158655f); + // Value for the 84.1345% mark + float binHiF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binHi, + totalDataPoints * 0.841345f); + psTrace(__func__, 6, "The exact 15.8655%% and 84.1345%% data point positions are: (%f, %f)\n", + binLoF32, binHiF32); + #else // This code basically interpolates to find the positions exactly. - // - if (1) { - psTrace(__func__, 6, "binLo is %d. Nums at that bin and the next are (%.2f, %.2f)\n", - binLo, cumulativeRobustHistogram->nums->data.F32[binLo], - cumulativeRobustHistogram->nums->data.F32[binLo+1]); - psTrace(__func__, 6, "binHi is %d. Nums at that bin and the next are (%.2f, %.2f)\n", - binHi, cumulativeRobustHistogram->nums->data.F32[binHi], - cumulativeRobustHistogram->nums->data.F32[binHi+1]); - - psF32 deltaNums; - psF32 deltaBounds; - psF32 prevPixels; - psF32 percentNums; - psF32 base; - deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binLo+1] - cumulativeRobustHistogram->bounds->data.F32[binLo]; - if (binLo == 0) { - deltaNums = cumulativeRobustHistogram->nums->data.F32[0]; - prevPixels = 0; - } else { - deltaNums = cumulativeRobustHistogram->nums->data.F32[binLo] - cumulativeRobustHistogram->nums->data.F32[binLo-1]; - prevPixels = cumulativeRobustHistogram->nums->data.F32[binLo-1]; - } - percentNums = (totalDataPoints * 0.158655f) - prevPixels; - base = cumulativeRobustHistogram->bounds->data.F32[binLo]; - binLoF32 = base + (deltaBounds / deltaNums) * percentNums; - psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", - base, deltaBounds, deltaNums, prevPixels, percentNums); - - deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binHi+1] - cumulativeRobustHistogram->bounds->data.F32[binHi]; - if (binHi == 0) { - deltaNums = cumulativeRobustHistogram->nums->data.F32[0]; - prevPixels = 0; - } else { - deltaNums = cumulativeRobustHistogram->nums->data.F32[binHi] - cumulativeRobustHistogram->nums->data.F32[binHi-1]; - prevPixels = cumulativeRobustHistogram->nums->data.F32[binHi-1]; - } - percentNums = (totalDataPoints * 0.841345f) - prevPixels; - base = cumulativeRobustHistogram->bounds->data.F32[binHi]; - binHiF32 = base + (deltaBounds / deltaNums) * percentNums; - psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", - base, deltaBounds, deltaNums, prevPixels, percentNums); - psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32); - } - - // - // ADD: Step 5. - // Determine SIGMA as 1/2 of the distance between these positions. - // + psTrace(__func__, 6, "binLo is %d. Nums at that bin and the next are (%.2f, %.2f)\n", + binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo+1]); + psTrace(__func__, 6, "binHi is %d. Nums at that bin and the next are (%.2f, %.2f)\n", + binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]); + + float deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo]; + float deltaNums; + float prevPixels; + if (binLo == 0) { + deltaNums = cumulative->nums->data.F32[0]; + prevPixels = 0; + } else { + deltaNums = cumulative->nums->data.F32[binLo] - cumulative->nums->data.F32[binLo-1]; + prevPixels = cumulative->nums->data.F32[binLo-1]; + } + float percentNums = (totalDataPoints * 0.158655f) - prevPixels; + float base = cumulative->bounds->data.F32[binLo]; + float binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark + psTrace(__func__, 6, + "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", + base, deltaBounds, deltaNums, prevPixels, percentNums); + + deltaBounds = cumulative->bounds->data.F32[binHi+1] - cumulative->bounds->data.F32[binHi]; + if (binHi == 0) { + deltaNums = cumulative->nums->data.F32[0]; + prevPixels = 0; + } else { + deltaNums = cumulative->nums->data.F32[binHi] - cumulative->nums->data.F32[binHi-1]; + prevPixels = cumulative->nums->data.F32[binHi-1]; + } + percentNums = (totalDataPoints * 0.841345f) - prevPixels; + base = cumulative->bounds->data.F32[binHi]; + float binHiF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 84.1345% mark + psTrace(__func__, 6, + "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", + base, deltaBounds, deltaNums, prevPixels, percentNums); + psTrace(__func__, 6, + "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", + binLoF32, binHiF32); + #endif + + // ADD step 5: Determine SIGMA as 1/2 of the distance between these positions. sigma = (binHiF32 - binLoF32) / 2.0; psTrace(__func__, 6, "The current sigma is %f.\n", sigma); stats->robustStdev = sigma; - // - // ADD: Step 6. - // If the measured SIGMA is less than 2 times the bin size, exclude - // points which are more than 25 bins from the median, - // recalculate the bin size, and perform the algorithm again. - // + // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more + // than 25 bins from the median, recalculate the bin size, and perform the algorithm again. if (sigma < (2 * binSize)) { psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize); - psS32 maskLo = PS_MAX(0, (binMedian - 25)); - psS32 maskHi = PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25)); - psF32 medianLo = robustHistogram->bounds->data.F32[maskLo]; - psF32 medianHi = robustHistogram->bounds->data.F32[maskHi]; - psTrace(__func__, 6, "Masking data more than 25 bins from the median (%.2f)\n", stats->robustMedian); - psTrace(__func__, 6, "The median is at bin number %d. We mask bins outside the bin range (%d:%d)\n", + long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region + long maskHi = PS_MIN(histogram->bounds->n - 1, (binMedian + 25)); // High index for masking + psF32 medianLo = histogram->bounds->data.F32[maskLo]; // Value at low index + psF32 medianHi = histogram->bounds->data.F32[maskHi]; // Value at high index + psTrace(__func__, 6, "Masking data more than 25 bins from the median\n"); + psTrace(__func__, 6, + "The median is at bin number %d. We mask bins outside the bin range (%d:%d)\n", binMedian, maskLo, maskHi); psTrace(__func__, 6, "Masking data outside (%f %f)\n", medianLo, medianHi); - for (psS32 i = 0 ; i < myVector->n ; i++) { + for (long i = 0 ; i < myVector->n ; i++) { if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) { - tmpMaskVec->data.U8[i] = 1; + mask->data.U8[i] = 0xff; psTrace(__func__, 6, "Masking element %d is %f\n", i, myVector->data.F32[i]); } } - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - iterate++; - } else { + // Free the histograms; they will be recreated on the next iteration, with new bounds + psFree(histogram); + psFree(cumulative); + } else { + // We've got the bin size correct now psTrace(__func__, 6, "*************: No more iteration. sigma is %f\n", sigma); - iterate = 0; - } - } - - // - // ADD: Step 7. - // Find the bins which contains the 25% and 75% data points. - // - psS32 binLo25 = 0; - psS32 binHi25 = 0; + iterate = -1; + } + } + psFree(histogram); + + // ADD step 7: Find the bins which contains the 25% and 75% data points. + long binLo25 = 0; // Bin containing the 25% value + long binHi25 = 0; // Bin containing the 75% value tmpScalar.data.F32 = totalDataPoints * 0.25f; - if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) { // Special case where its in first bin. The last bin should be handled automatically. binLo25 = 0; } else { - binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); + binLo25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); } tmpScalar.data.F32 = totalDataPoints * 0.75f; - if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) { // Special case where its in first bin. The last bin should be handled automatically. binHi25 = 0; } else { - binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); + binHi25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); } if ((binLo25 < 0) || (binHi25 < 0)) { psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25 and 75 percent data points\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpMaskVec); + psFree(statsMinMax); + psFree(cumulative); + psFree(mask); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; + } + psTrace(__func__, 6, "The 25-percent and 75-precent data point bins are (%d, %d).\n", binLo25, binHi25); + + // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile + // positions. + psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, + binLo25, totalDataPoints * 0.25f); + psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, + binHi25, totalDataPoints * 0.75f); + psFree(cumulative); + + if (isnan(binLo25F32) || isnan(binHi25F32)) { + psError(PS_ERR_UNKNOWN, false, + "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n"); + psFree(statsMinMax); + psFree(mask); psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - psTrace(__func__, 6, "The 25-percent and 75-precent data point bins are (%d, %d).\n", binLo25, binHi25); - - // - // ADD: Step 8. - // Interpolate to find these two positions exactly: these are the upper - // and lower quartile positions. - // - psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&cumulativeRobustHistogram->bounds, - *(psVector* *)&cumulativeRobustHistogram->nums, - binLo25, - totalDataPoints * 0.25f); - - psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&cumulativeRobustHistogram->bounds, - *(psVector* *)&cumulativeRobustHistogram->nums, - binHi25, - totalDataPoints * 0.75f); - - if (isnan(binLo25F32) || isnan(binHi25F32)) { - psError(PS_ERR_UNKNOWN, false, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); + return false; } stats->robustLQ = binLo25F32; stats->robustUQ = binHi25F32; - psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32); - psS32 N50 = 0; - for (psS32 i = 0 ; i < myVector->n ; i++) { - // XXX: use maskVal here? - if ((0 == tmpMaskVec->data.U8[i]) && - (binLo25F32 <= myVector->data.F32[i]) && - (binHi25F32 >= myVector->data.F32[i])) { + psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", + binLo25F32, binHi25F32); + long N50 = 0; + for (long i = 0 ; i < myVector->n ; i++) { + if (!(mask->data.U8[i] & maskVal) && + (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) { N50++; } @@ -1740,13 +1440,11 @@ psTrace(__func__, 6, "The robustN50 is %d.\n", N50); - // ************************************************************************ - if ((stats->options & PS_STAT_FITTED_MEAN) || - (stats->options & PS_STAT_FITTED_STDEV)) { - // + + // Fitted statistics + if (stats->options & PS_STAT_FITTED_MEAN || stats->options & PS_STAT_FITTED_STDEV) { // New algorithm for calculated mean and stdev. - // // XXX: Where is this documented? - psF32 dN = (psF32) (0.17 * N50); + psF32 dN = 0.17 * N50; if (dN < 1.0) { dN = 1.0; @@ -1756,116 +1454,102 @@ psF32 newBinSize = sigma / dN; - // // Determine the min/max of the vector (which prior outliers masked out) - // - rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); - rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); - if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { + int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values + min = statsMinMax->min; + max = statsMinMax->max; + if (numValid == 0 || isnan(min) || isnan(max)) { psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - - // + psFree(statsMinMax); + psFree(mask); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; + } + // Calculate the number of bins. - // - numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize); - psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max); + long numBins = (max - min) / newBinSize; + psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", min, max); psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize); psTrace(__func__, 6, "The numBins is %d\n", numBins); - psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); - newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1); + psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) + histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal); if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(newHistogram->nums); - } - - // - // FITTED STATISTICS HERE - // - // Smooth the resulting histogram with a Gaussian with sigma_s = 1 - // bin. - // - psF32 sigma_s = newBinSize; - - psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s); + PS_VECTOR_PRINT_F32(histogram->nums); + } + + // Smooth the resulting histogram with a Gaussian with sigma_s = 1 bin. + psVector *smoothed = vectorSmoothHistGaussian(histogram, newBinSize); // Smoothed histogram if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(newHistogramSmoothed); - } - - // - // Find the bin with the peak value in the range 2 sigma of the - // robust histogram median. - // - - psS32 binMin = 0; - psS32 binMax = 0; + PS_VECTOR_PRINT_F32(smoothed); + } + + // Find the bin with the peak value in the range 2 sigma of the robust histogram median. + long binMin = 0; // Low bin to check + long binMax = 0; // High bin to check tmpScalar.data.F32 = stats->robustMedian - (2.0 * sigma); - if (tmpScalar.data.F32 <= newHistogram->bounds->data.F32[0]) { + if (tmpScalar.data.F32 <= histogram->bounds->data.F32[0]) { binMin = 0; } else { - binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar); + binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar); } tmpScalar.data.F32 = stats->robustMedian + (2.0 + sigma); - if (tmpScalar.data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) { - binMax = newHistogram->bounds->n-1; - } else { - binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar); + if (tmpScalar.data.F32 >= histogram->bounds->data.F32[histogram->bounds->n-1]) { + binMax = histogram->bounds->n-1; + } else { + binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar); } if ((binMin < 0) || (binMax < 0)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n"); - psFree(tmpMaskVec); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpStatsMinMax); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - - psS32 binNum = binMin; - psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum]; + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +/- 2.0 * sigma bins\n"); + psFree(statsMinMax); + psFree(histogram); + psFree(mask); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; + } + + long binNum = binMin; // Bin number containing the peak + psF32 binMaxNums = smoothed->data.F32[binNum]; // The peak value for (psS32 i = binMin+1 ; i <= binMax ; i++) { - if (newHistogramSmoothed->data.F32[i] > binMaxNums) { + if (smoothed->data.F32[i] > binMaxNums) { binNum = i; - binMaxNums = newHistogramSmoothed->data.F32[i]; + binMaxNums = smoothed->data.F32[i]; } } psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums); - // - // Fit a Gaussian to the bins in the range 20 sigma of the robust - // histogram median. - // + // Fit a Gaussian to the bins in the range 20 sigma of the robust histogram median. tmpScalar.data.F32 = stats->robustMedian - (20.0 * sigma); - if (tmpScalar.data.F32 < tmpStatsMinMax->min) { + if (tmpScalar.data.F32 < min) { binMin = 0; } else { - binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar); - // XXX: check for errors here. + binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar); } tmpScalar.data.F32 = stats->robustMedian + (20.0 * sigma); - if (tmpScalar.data.F32 > tmpStatsMinMax->max) { - binMax = newHistogramSmoothed->n - 1; - } else { - binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar); - // XXX: check for errors here. - } - psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); - psArray *x = psArrayAlloc((1 + (binMax - binMin))); - psS32 j = 0; + if (tmpScalar.data.F32 > max) { + binMax = smoothed->n - 1; + } else { + binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar); + } + if (binMin < 0 || binMax < 0) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +/- 20.0 * sigma bins\n"); + psFree(statsMinMax); + psFree(histogram); + psFree(mask); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; + } + + // Generate the variables that will be used in the Gaussian fitting + psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates + psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates y->n = y->nalloc; - - for (psS32 i = binMin ; i <= binMax ; i++) { - y->data.F32[j] = newHistogramSmoothed->data.F32[i]; - x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32); - x->n++; - ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i); - ((psVector *) x->data[j])->n++; - j++; + x->n = x->nalloc; + for (long i = binMin, j = 0; i <= binMax ; i++, j++) { + y->data.F32[j] = smoothed->data.F32[i]; + psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value + ordinate->n = 1; + ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i); + x->data[j] = ordinate; } if (psTraceGetLevel(__func__) >= 8) { @@ -1873,30 +1557,27 @@ PS_VECTOR_PRINT_F32(y); } - - rc = p_psVectorMin(y, NULL, 0, tmpStatsMinMax); - rc|= p_psVectorMax(y, NULL, 0, tmpStatsMinMax); - if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { + psFree(smoothed); + psFree(histogram); + + numValid = vectorMinMax(y, NULL, 0, statsMinMax); // Number of valid values + min = statsMinMax->min; + max = statsMinMax->max; + if (numValid == 0 || isnan(min) || isnan(max)) { psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - psFree(tmpMaskVec); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpStatsMinMax); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - - // + psFree(mask); + psFree(histogram); + psFree(statsMinMax); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; + } + // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) - // XXX: Use the normalize routines for this. - // - for (psS32 i = 0 ; i < y->n ; i++) { - y->data.F32[i]= (y->data.F32[i] - tmpStatsMinMax->min) / (tmpStatsMinMax->max - tmpStatsMinMax->min); - } - - // - psMinimization *min = psMinimizationAlloc(100, 0.01); - psVector *params = psVectorAlloc(2, PS_TYPE_F32); + p_psNormalizeVectorRange(y, 0.0, 1.0); + + // Fit a Gaussian to the data. + psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information + psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian params->n = params->nalloc; - // Initial guess for the mean ([0]) and standard dev. + // Initial guess for the mean (index 0) and standard dev (index 1). params->data.F32[0] = stats->robustMedian; params->data.F32[1] = sigma; @@ -1905,22 +1586,14 @@ PS_VECTOR_PRINT_F32(y); } - - // - // Fit a Gaussian to the data. - // - rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D); - if (rcBool != true) { + if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) { psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(newHistogram); + psFree(statsMinMax); psFree(x); psFree(y); - psFree(min); + psFree(minimizer); psFree(params); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); + psFree(mask); + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return false; } if (psTraceGetLevel(__func__) >= 8) { @@ -1928,34 +1601,26 @@ } - // - // The fitted mean mean_r is derived directly from the fitted - // Gaussian mean. - // + // The fitted mean is the Gaussian mean. stats->fittedMean = params->data.F32[0]; psTrace(__func__, 6, "The fitted mean is %f.\n", params->data.F32[0]); - // - // The fitted standard deviation, SIGMA_r is determined by - // subtracting the smoothing scale in quadrature: - // SIGMA_r^2 = SIGMA^2 - sigma_s^2 - // - stats->fittedStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s)); + // The fitted standard deviation, SIGMA_r is determined by subtracting the smoothing scale in + // quadrature: SIGMA_r^2 = SIGMA^2 - sigma_s^2 + stats->fittedStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(newBinSize)); psTrace(__func__, 6, "The fitted stdev is %f.\n", stats->fittedStdev); - psFree(newHistogramSmoothed); - psFree(newHistogram); + // Clean up after fitting psFree(x); psFree(y); - psFree(min); + psFree(minimizer); psFree(params); } - psFree(tmpStatsMinMax); - psFree(cumulativeRobustHistogram); - psFree(tmpMaskVec); - psFree(robustHistogram); + // Clean up + psFree(statsMinMax); + psFree(mask); psTrace(__func__, 4, "---- %s(0) end ----\n", __func__); - return(0); + return true; } @@ -1968,4 +1633,5 @@ static void histogramFree(psHistogram* myHist); +// We keep statsFree so that we can identify statistics pointers from the memblock static void statsFree(psStats *newStruct) { @@ -2038,11 +1704,7 @@ PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(upper, lower, NULL); - psS32 i = 0; // Loop index variable - psHistogram* newHist = NULL; // The new histogram structure - psF32 binSize = 0.0; // The histogram bin size - - // 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)); + // Allocate memory for the new histogram structure. If there are N bins, then there are N+1 bounds to + // those bins. + psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure psMemSetDeallocator(newHist, (psFreeFunc) histogramFree); psVector* newBounds = psVectorAlloc(n + 1, PS_TYPE_F32); @@ -2051,9 +1713,8 @@ // Calculate the bounds for each bin. - binSize = (upper - lower) / (psF32)n; - // XXX: Is the following necessary? It prevents the max data point - // from being in a non-existant bin. + psF32 binSize = (upper - lower) / (psF32)n; // The histogram bin size + // XXX: 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++) { + for (long i = 0; i < n + 1; i++) { newBounds->data.F32[i] = lower + (binSize * (psF32)i); } @@ -2061,5 +1722,5 @@ // Allocate the bins, and initialize them to zero. newHist->nums = psVectorAlloc(n, PS_TYPE_F32); - for (i = 0; i < newHist->nums->nalloc; i++) { + for (long i = 0; i < newHist->nums->nalloc; i++) { newHist->nums->data.F32[i] = 0.0; newHist->nums->n++; @@ -2072,5 +1733,5 @@ psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return (newHist); + return newHist; } @@ -2091,14 +1752,11 @@ PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(bounds->n, 2, NULL); - psHistogram* newHist = NULL; // The new histogram structure - psS32 i; // Loop index variable - // Allocate memory for the new histogram structure. - newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); + psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure psMemSetDeallocator(newHist, (psFreeFunc) histogramFree); psVector* newBounds = psVectorAlloc(bounds->n, PS_TYPE_F32); newHist->bounds = newBounds; newBounds->n = newHist->bounds->nalloc; - for (i = 0; i < bounds->n; i++) { + for (long i = 0; i < bounds->n; i++) { newBounds->data.F32[i] = bounds->data.F32[i]; } @@ -2107,5 +1765,5 @@ // then there are N-1 bins. newHist->nums = psVectorAlloc((bounds->n) - 1, PS_TYPE_F32); - for (i = 0; i < newHist->nums->nalloc; i++) { + for (long i = 0; i < newHist->nums->nalloc; i++) { newHist->nums->data.F32[i] = 0.0; newHist->nums->n++; @@ -2140,37 +1798,33 @@ the data point as a boxcar PDF and update a range of points surrounding the histogram bin which contains the point. The width of that boxcar is defined -as 2.35 * error. Inputs: - binNum: the bin number of the data point in the histogram - out: the histogram structure - data: the data point value - error: the error in that data point +as 2.35 * error. XXX: Must test this. *****************************************************************************/ -psS32 UpdateHistogramBins(psS32 binNum, - psHistogram* out, - psF32 data, - psF32 error) +static bool UpdateHistogramBins(long binNum, // Bin number of the data point + psHistogram* out, // The histogram to be updated + psF32 data, // The data point value + psF32 error // The error in the data point + ) { psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(out, -1); - PS_ASSERT_PTR_NON_NULL(out->bounds, -1); - PS_ASSERT_PTR_NON_NULL(out->nums, -1); - PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, ((out->nums->n)-1), -2); - PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(error, 0.0, -3); - - PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4); - - psF32 boxcarWidth = 2.35 * error; + PS_ASSERT_PTR_NON_NULL(out, false); + PS_ASSERT_PTR_NON_NULL(out->bounds, false); + PS_ASSERT_PTR_NON_NULL(out->nums, false); + PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, ((out->nums->n)-1), false); + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(error, 0.0, false); + PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0], + out->bounds->data.F32[(out->bounds->n)-1], false); + + psF32 boxcarWidth = 2.35 * error; // Width of the boxcar psF32 boxcarCenter = (out->bounds->data.F32[binNum] + - out->bounds->data.F32[binNum+1]) / 2.0; - psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); - psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); - psS32 bin; - psS32 boxcarLeftBinNum = 0; - psS32 boxcarRightBinNum = 0; + out->bounds->data.F32[binNum+1]) / 2.0; // Centre of the boxcar + psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); // Left endpoint of the boxcar for the PDF + psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); // Right endpoint of the boxcar for the PDF + psS32 boxcarLeftBinNum = 0; // Bin number for left endpoint + psS32 boxcarRightBinNum = 0; // Bin number for right endpoint // Determine the left endpoint of the boxcar for the PDF. - for (bin=binNum ; bin >= 0 ; bin--) { + for (long bin = binNum; bin >= 0; bin--) { if (out->nums->data.F32[bin] <= boxcarLeft) { boxcarLeftBinNum = bin; @@ -2180,5 +1834,5 @@ // Determine the right endpoint of the boxcar for the PDF. - for (bin=binNum ; bin < out->nums->n ; bin++) { + for (long bin = binNum; bin < out->nums->n; bin++) { if (out->nums->data.F32[bin] >= boxcarRight) { boxcarRightBinNum = bin; @@ -2187,40 +1841,29 @@ } - // // If the boxcar fits entirely inside this bin, then simply add 1.0 to the // bin and return. - // if (boxcarLeftBinNum == boxcarRightBinNum) { - out->nums->data.F32[binNum]+= 1.0; - psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); - return(0); - } - - // - // If we get here, multiple bins must be updated. We handle the left - // endpoint, and right endpoint differently. - // - out->nums->data.F32[boxcarLeftBinNum]+= + out->nums->data.F32[binNum] += 1.0; + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return true; + } + + // If we get here, multiple bins must be updated. We handle the left-most endpoint, and right-most + // endpoints differently. + out->nums->data.F32[boxcarLeftBinNum] += (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth; - // // Loop through the center bins, if any. - // - for (bin = boxcarLeftBinNum + 1 ; bin < (boxcarRightBinNum - 1) ; bin++) { - out->nums->data.F32[bin]+= + for (long bin = boxcarLeftBinNum + 1; bin < (boxcarRightBinNum - 1); bin++) { + out->nums->data.F32[bin] += (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth; } - // // Handle the right endpoint differently. - // out->nums->data.F32[boxcarRightBinNum]+= (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth; - // - // Return 0 on success. - // - psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); - return(0); + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return true; } @@ -2255,40 +1898,37 @@ PS_ASSERT_INT_NONNEGATIVE(out->nums->n, NULL); PS_ASSERT_VECTOR_NON_NULL(values, out); - if (mask != NULL) { + if (mask) { PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, NULL); PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); } - if (errors != NULL) { + if (errors) { PS_ASSERT_VECTORS_SIZE_EQUAL(values, errors, NULL); PS_ASSERT_VECTOR_TYPE(errors, values->type.type, NULL); } - psS32 i = 0; // Loop index variable - psF32 binSize = 0.0; // Histogram bin size - psS32 binNum = 0; // A temporary bin number - psS32 numBins = 0; // The total number of bins + long binNum = 0; // A temporary bin number + long numBins = out->nums->n; // The total number of bins psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; - psVector* inF32 = NULL; - psVector* errorsF32 = NULL; - psS32 mustFreeVectorIn = 1; - psS32 mustFreeVectorErrors = 1; // Convert input and errors vectors to F32 if necessary. - inF32 = p_psConvertToF32((psVector *) values); - if (inF32 == NULL) { - inF32 = (psVector *) values; - mustFreeVectorIn = 0; - } - errorsF32 = p_psConvertToF32((psVector *) errors); - if (errorsF32 == NULL) { - errorsF32 = (psVector *) errors; - mustFreeVectorErrors = 0; - } - - numBins = out->nums->n; - for (i = 0; i < inF32->n; i++) { + psVector* inF32 = NULL; // F32 version of input vector + if (values->type.type == PS_TYPE_F32) { + inF32 = psMemIncrRefCounter((psPtr)values); + } else { + values = psVectorCopy(NULL, values, PS_TYPE_F32); + } + psVector* errorsF32 = NULL; // F32 version of errors vector + if (errors) { + if (errors->type.type == PS_TYPE_F32) { + errorsF32 = psMemIncrRefCounter((psPtr)errors); + } else { + errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32); + } + } + + for (long i = 0; i < inF32->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 || (mask && (!(mask->data.U8[i] & maskVal)))) { if (inF32->data.F32[i] < out->bounds->data.F32[0]) { // If this pixel is below minimum value, count it, then skip. @@ -2301,11 +1941,10 @@ // number is trivial. if (out->uniform == true) { - binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; - binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); - if (errorsF32 != NULL) { - psS32 rc = UpdateHistogramBins(binNum, out, inF32->data.F32[i], - errorsF32->data.F32[i]); - if (rc < 0) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram bins with the errors vector.\n"); + float binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; // Histogram bin size + binNum = (inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize; + if (errorsF32) { + if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errorsF32->data.F32[i])) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram " + "bins with the errors vector.\n"); } } else { @@ -2329,12 +1968,10 @@ } else { if (errorsF32 != NULL) { - psS32 rc = UpdateHistogramBins(binNum, out, - inF32->data.F32[i], - errors->data.F32[i]); - if (rc < 0) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram bins with the errors vector.\n"); + if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram " + "bins with the errors vector.\n"); } } else { - (out->nums->data.F32[binNum])+= 1.0; + out->nums->data.F32[binNum] += 1.0; } } @@ -2344,102 +1981,9 @@ } - if (mustFreeVectorIn == 1) { - psFree(inF32); - } - if (mustFreeVectorErrors == 1) { - psFree(errorsF32); - } + psFree(inF32); + psFree(errorsF32); psTrace(__func__, 3, "---- %s() end ----\n", __func__); return (out); -} - -/****************************************************************************** -p_psConvertToF32(in): this is the cheap way to support a variety of vector -data types: we simply convert the input vector to F32 at the beginning, and -write all of our functions in F32. If the vast majority of all vector stat -operations are F32 (or any other single type), then this is probably the -best way to go. Otherwise, when the algorithms stablize, we will then macro -everything and put type support in the various stat functions. - -XXX: Should the default data type be F64? Since we are buying Opterons... - -XXX: Use psLib functions instead. - *****************************************************************************/ -psVector* p_psConvertToF32(psVector* in) -{ - psTrace(__func__, 4,"---- %s() begin ----\n", __func__); - if (in == NULL) { - return(NULL); - } - psS32 i = 0; - psVector* tmp = NULL; - - if (in->type.type == PS_TYPE_S8) { - tmp = psVectorAlloc(in->n, PS_TYPE_F32); - for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (psF32)in->data.S8[i]; - tmp->n++; - } - } else if (in->type.type == PS_TYPE_S16) { - tmp = psVectorAlloc(in->n, PS_TYPE_F32); - for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (psF32) in->data.S16[i]; - tmp->n++; - } - } else 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] = (psF32)in->data.S32[i]; - tmp->n++; - } - } else if (in->type.type == PS_TYPE_S64) { - tmp = psVectorAlloc(in->n, PS_TYPE_F32); - for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (psF32)in->data.S64[i]; - tmp->n++; - } - } else if (in->type.type == PS_TYPE_U8) { - tmp = psVectorAlloc(in->n, PS_TYPE_F32); - for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (psF32)in->data.U8[i]; - tmp->n++; - } - } else if (in->type.type == PS_TYPE_U16) { - tmp = psVectorAlloc(in->n, PS_TYPE_F32); - for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (psF32)in->data.U16[i]; - tmp->n++; - } - } 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] = (psF32)in->data.U32[i]; - tmp->n++; - } - } else if (in->type.type == PS_TYPE_U64) { - tmp = psVectorAlloc(in->n, PS_TYPE_F32); - for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (psF32)in->data.U64[i]; - tmp->n++; - } - } 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] = (psF32)in->data.F64[i]; - tmp->n++; - } - } else if (in->type.type == PS_TYPE_F32) { - // do nothing - } else { - char* strType; - PS_TYPE_NAME(strType, in->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED, - strType); - } - - psTrace(__func__, 4,"---- %s() end ----\n", __func__); - return (tmp); } @@ -2469,36 +2013,35 @@ PS_ASSERT_PTR_NON_NULL(stats, NULL); PS_ASSERT_VECTOR_NON_NULL(in, NULL); - if (mask != NULL) { + if (mask) { PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, stats); PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, stats); } - if (errors != NULL) { + if (errors) { PS_ASSERT_VECTORS_SIZE_EQUAL(errors, in, stats); PS_ASSERT_VECTOR_TYPE(errors, in->type.type, stats); } - // XXX: Assert that "in" is F64, F32, U16, or S8 - - psVector* inF32 = NULL; - psVector* errorsF32 = NULL; - psS32 mustFreeVectorIn = 1; - psS32 mustFreeVectorErrors = 1; - - inF32 = p_psConvertToF32((psVector *) in); - if (inF32 == NULL) { - // printf("\n\ninF32 is NULL\n\n"); - // inF32 = (psVector *) in; - // mustFreeVectorIn = 0; - inF32 = psVectorCopy(inF32, in, PS_TYPE_F32); - inF32->n = inF32->nalloc; - mustFreeVectorIn = 1; - } - // else printf("\ninF32 has n=%ld, nalloc=%ld\n", inF32->n, inF32->nalloc); - errorsF32 = p_psConvertToF32((psVector *) errors); - if (errorsF32 == NULL) { - errorsF32 = (psVector *) errors; - mustFreeVectorErrors = 0; - } - // inF32->n = inF32->nalloc; - // errorsF32->n = errorsF32->nalloc; + + // Convert types, as necessary + psVector *inF32 = NULL; // Input vector of values, F32 version + if (in->type.type == PS_TYPE_F32) { + inF32 = psMemIncrRefCounter((psPtr)in); + } else { + inF32 = psVectorCopy(NULL, in, PS_TYPE_F32); + } + psVector *errorsF32 = NULL; // Input vector of errors, F32 version + if (errors) { + if (errors->type.type == PS_TYPE_F32) { + errorsF32 = psMemIncrRefCounter((psPtr)errors); + } else { + errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32); + } + } + psVector *maskU8 = NULL; // Input mask vector, U8 version + if (mask) { + // We want a copy of the mask, since we may change values + maskU8 = psVectorCopy(NULL, mask, PS_TYPE_U8); + } + + if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) { PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->max, stats->min, stats); @@ -2511,6 +2054,6 @@ // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_MEAN) { - if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n"); + if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n"); stats->sampleMean = NAN; } @@ -2518,8 +2061,10 @@ // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_MEDIAN) { - if (false == p_psVectorSampleMedian(inF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian() returned an error.\n"); + if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) { + if (!vectorSampleMedian(inF32, maskU8, maskVal, stats)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMedian() returned an error.\n"); stats->sampleMedian = NAN; + stats->sampleUQ = NAN; + stats->sampleLQ = NAN; } } @@ -2527,30 +2072,20 @@ // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_STDEV) { - if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n"); + if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n"); stats->sampleMean = NAN; } else { - p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats); + vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats); } } // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_QUARTILE) { - if (false == p_psVectorSampleQuartiles(inF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleQuartiles() returned an error.\n"); - stats->sampleLQ = NAN; - stats->sampleUQ = NAN; - } - } - - // ************************************************************************ - if ((stats->options & PS_STAT_ROBUST_MEDIAN) || - (stats->options & PS_STAT_ROBUST_STDEV) || - (stats->options & PS_STAT_ROBUST_QUARTILE) || - (stats->options & PS_STAT_FITTED_MEAN) || - (stats->options & PS_STAT_FITTED_STDEV)) { - if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) { + if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE | + PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { + if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED); psFree(stats); + psFree(inF32); + psFree(errorsF32); psTrace(__func__, 3,"---- %s(NULL) end ----\n", __func__); return(NULL); @@ -2560,6 +2095,5 @@ // ************************************************************************ if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { - psS32 rc = p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats); - if (rc < 0) { + if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n"); stats->clippedMean = NAN; @@ -2569,25 +2103,13 @@ // ************************************************************************ - if (stats->options & PS_STAT_MAX) { - if (0 != p_psVectorMax(inF32, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector maximum"); - stats->max = NAN; - } - } - - // ************************************************************************ - if (stats->options & PS_STAT_MIN) { - if (0 != p_psVectorMin(inF32, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum"); - stats->min = NAN; - } - } - - if (mustFreeVectorIn == 1) { - psFree(inF32); - } - if (mustFreeVectorErrors == 1) { - psFree(errorsF32); - } + if (stats->options & (PS_STAT_MAX | PS_STAT_MIN)) { + if (vectorMinMax(inF32, maskU8, maskVal, stats) == 0) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum and maximum"); + } + } + + psFree(inF32); + psFree(errorsF32); + psFree(maskU8); psTrace(__func__, 3,"---- %s() end ----\n", __func__); return (stats);