Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 1276) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 1277) @@ -9,11 +9,11 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-07-15 23:52:34 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-07-22 23:40:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ /*****************************************************************************/ -/* INCLUDE FILES */ +/* INCLUDE FILES */ /*****************************************************************************/ #include @@ -41,5 +41,5 @@ #include "psMatrix.h" /*****************************************************************************/ -/* DEFINE STATEMENTS */ +/* DEFINE STATEMENTS */ /*****************************************************************************/ #define MAX_LMM_ITERATIONS 100 @@ -85,5 +85,5 @@ /*****************************************************************************/ -/* TYPE DEFINITIONS */ +/* TYPE DEFINITIONS */ /*****************************************************************************/ typedef struct @@ -113,5 +113,5 @@ /*****************************************************************************/ -/* GLOBAL VARIABLES */ +/* GLOBAL VARIABLES */ /*****************************************************************************/ @@ -119,5 +119,5 @@ /*****************************************************************************/ -/* FILE STATIC VARIABLES */ +/* FILE STATIC VARIABLES */ /*****************************************************************************/ @@ -125,5 +125,5 @@ /*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - LOCAL */ +/* FUNCTION IMPLEMENTATION - LOCAL */ /*****************************************************************************/ @@ -558,5 +558,5 @@ /*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - PUBLIC */ +/* FUNCTION IMPLEMENTATION - PUBLIC */ /*****************************************************************************/ @@ -567,5 +567,5 @@ This routine must minimize an arbitrary function; it determines the set of -parameters of that function such that the. +parameters of that function such that the ... *****************************************************************************/ psVector * @@ -681,5 +681,6 @@ /****************************************************************************** - This routine must minimize an arbitrary function. + This routine must determine the parameters of an arbitrary function + such that they best fit the supplied data points. *****************************************************************************/ psVector * @@ -795,10 +796,17 @@ do { iter++; + for (i=0;in;i++) { + printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i)); + } // Perform an iteration of the GSL solver. status = gsl_multifit_fdfsolver_iterate(s); - // printf ("psMinimize() status = %s\n", gsl_strerror(status)); + printf("gsl_multifit_fdfsolver_iterate() status is %s\n", gsl_strerror(status)); + for (i=0;in;i++) { + printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i)); + } + // If there was a problem, abort. if (status) { - psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n"); + psAbort(__func__, "gsl_multifit_fdfsolver_iterate(%s)\n", gsl_strerror(status)); } @@ -811,4 +819,5 @@ // as specified in the ADD. *chiSq = gsl_blas_dnrm2(s->f); + printf("psMinimize.c: chiSq is %.3f\n", *chiSq); if (fabs(*chiSq - chiSqOld) < 1.0) { status = GSL_SUCCESS; @@ -855,5 +864,5 @@ This routine must fit a polynomial of degree myPoly to the data points (x, y) and return the coefficients of that polynomial, as well as the - error for each data poiny (yErr). + error for each data point (yErr). NOTE: yErr is currently ignored. @@ -877,4 +886,9 @@ psVector *xSums = NULL; + // printf("psGetArrayPolynomial()\n"); + // for (i=0;in;i++) { + // printf("(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]); + // } + PS_CHECK_NULL_1DPOLY(myPoly); PS_CHECK_NULL_VECTOR(x); @@ -929,5 +943,11 @@ for(k=0;k<(polyOrder);k++) { myPoly->coeff[k] = coeffs->data.F64[k]; - } + // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]); + } + + + // for (i=0;in;i++) { + // printf("HMMM: psEvalPolynomial1D(%f) is %f\n", x->data.F64[i], psEvalPolynomial1D(x->data.F64[i], myPoly)); + // } psFree(A); @@ -939,4 +959,4 @@ psFree(xSums); - return(NULL); -} + return(myPoly); +} Index: /trunk/psLib/src/dataManip/psStats.c =================================================================== --- /trunk/psLib/src/dataManip/psStats.c (revision 1276) +++ /trunk/psLib/src/dataManip/psStats.c (revision 1277) @@ -5,6 +5,9 @@ #include #include - +/*****************************************************************************/ +/* INCLUDE FILES */ +/*****************************************************************************/ #include "psMemory.h" +#include "psImage.h" #include "psVector.h" #include "psTrace.h" @@ -12,23 +15,1486 @@ #include "psAbort.h" #include "psStats.h" - -#include "float.h" -#define DEFAULT_ROBUST_SIZE_THRESHOLD 30000 // Vectors that are large than this -// will use robust statistical methods. -#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. -#define CLIPPED_NUM_ITER_LB 1 -#define CLIPPED_NUM_ITER_UB 10 -#define CLIPPED_SIGMA_LB 1.0 -#define CLIPPED_SIGMA_UB 10.0 -#define true 1 -#define false 0 -#define MYMAXFLOAT HUGE - -#ifndef DOXYGEN +#include "psSort.h" +#include "psMinimize.h" +#include "psFunctions.h" +======= + + /*****************************************************************************/ + /* DEFINE STATEMENTS */ + /*****************************************************************************/ + #define DEFAULT_ROBUST_SIZE_THRESHOLD 30000 // Vectors that are large than this + // will use robust statistical methods. + #define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. + #define CLIPPED_NUM_ITER_LB 1 + #define CLIPPED_NUM_ITER_UB 10 + #define CLIPPED_SIGMA_LB 1.0 + #define CLIPPED_SIGMA_UB 10.0 + #define true 1 + #define false 0 + #define MYMAXFLOAT HUGE + #define MAX_ITERATIONS 20 + + #ifndef DOXYGEN + void p_psVectorRobustStats(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats); +#endif + +/** Preprocessor macro to generate error on an incorrect type */ +#define PS_CHECK_VECTOR_TYPE(NAME, TYPE) \ +if (NAME->type.type != TYPE) { \ + psError(__func__,"Invalid operation: %s has incorrect type.", #NAME); \ +} + +/** Preprocessor macro to generate error on a NULL vector */ +#define PS_CHECK_NULL_VECTOR(NAME) \ +if (NAME == NULL || NAME->data.V == NULL) { \ + psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \ +} + +/** Preprocessor macro to generate error for zero length vector */ +#define PS_CHECK_EMPTY_VECTOR(NAME) \ +if (NAME->n < 1) { \ + psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \ +} + +/** Preprocessor macro to generate error on differing size vectors */ +#define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2) \ +if (VEC1->n != VEC2->n) { \ + psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \ +} + +/*****************************************************************************/ +/* TYPE DEFINITIONS */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +/***************************************************************************** +p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal +function that simply prints a vector to STDOUT. Used primarily for +debugging. + *****************************************************************************/ + +void p_psVectorPrint(psVector *myVector, + psVector *maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable. + + for (i=0;in;i++) { + if (maskVector != NULL) + printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]); + else + printf("Element %d is %f\n", i, myVector->data.F32[i]); + } +} + +/****************************************************************************** + ****************************************************************************** + ****************************************************************************** + MISC PRIVATE STATISTICAL FUNCTIONS + + NOTE: it is assumed that any call to these statistical functions will + have been preceded by a call to the psVectorStats() function. Various + sanity tests will only be performed in psVectorStats(). + Is the mask vector the same length as the data vector? + Is the mask vector of type PS_TYPE_U8? + Is the stats data structure NULL? + Is the in data structure NULL? + Is the in data structure of type PS_TYPE_F32? + ****************************************************************************** + ****************************************************************************** + *****************************************************************************/ + + +/****************************************************************************** +p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the +mean of the input vector. +Inputs +Returns + NULL +ASSUMPTION: the mean is always calculated exactly. Robust means are never +calculated in this routine. + *****************************************************************************/ +void p_psVectorSampleMean(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + float mean = 0.0; // The mean + int count = 0; // # of points in this mean? + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + // If PS_STAT_USE_RANGE is requested, then we enter a slightly different + // loop. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + // Check if the data is with the specified range + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + mean+= myVector->data.F32[i]; + count++; + } + } + mean/= (float) count; + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + mean+= myVector->data.F32[i]; + count++; + } + } + mean/= (float) count; + } + } else { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + mean+= myVector->data.F32[i]; + count++; + } + } + mean/= (float) count; + } else { + for (i=0;in;i++) { + mean+= myVector->data.F32[i]; + } + mean/= (float) myVector->n; + } + } + + stats->sampleMean = mean; +} + +/****************************************************************************** +p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the +max of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorMax(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + float max = -MYMAXFLOAT; // The calculated maximum + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + // If PS_STAT_USE_RANGE is requested, then we enter a different loop. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + if ((myVector->data.F32[i] > max) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + max = myVector->data.F32[i]; + } + } + } + } else { + for (i=0;in;i++) { + if ((myVector->data.F32[i] > max) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + max = myVector->data.F32[i]; + } + } + } + } else { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + if (myVector->data.F32[i] > max) { + max = myVector->data.F32[i]; + } + } + } + } else { + for (i=0;in;i++) { + if (myVector->data.F32[i] > max) { + max = myVector->data.F32[i]; + } + } + } + } + + stats->max = max; +} + +/****************************************************************************** +p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the +minimum of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorMin(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + float min = MYMAXFLOAT; // The calculated maximum + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + if ((myVector->data.F32[i] < min) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + min = myVector->data.F32[i]; + } + } + } + } else { + for (i=0;in;i++) { + if ((myVector->data.F32[i] < min) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + min = myVector->data.F32[i]; + } + } + } + } else { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + if (myVector->data.F32[i] < min) { + min = myVector->data.F32[i]; + } + } + } + } else { + for (i=0;in;i++) { + if (myVector->data.F32[i] < min) { + min = myVector->data.F32[i]; + } + } + } + } + + stats->min = min; +} + +/****************************************************************************** +p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the +number of non-masked pixels in the vector that fall within the min/max +range, if given. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +int p_psVectorNValues(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + int numData = 0; // The number of data points + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + numData++; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + numData++; + } + } + } + } else { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + numData++; + } + } + } else { + numData = myVector->n; + } + } + return(numData); +} + + + +/****************************************************************************** +p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the +median of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorSampleMedian(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + psVector *unsortedVector = NULL; // Temporary vector + psVector *sortedVector = NULL; // Temporary vector + int i = 0; // Loop index variable + int count = 0; // # of points in this mean? + int nValues = 0; // # of points in vector + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + psStats *stats2 = NULL; // Temporary stats structure + + + // Determine if the number of data points exceed a threshold which will + // cause to generate robust stats, as opposed to exact stats. + + if (myVector->n > stats->sampleLimit) { + psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); + + // Calculate the robust quartiles. + stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); + + // Store the robust quartiles into the sample quartile members. + stats->sampleMedian = stats2->robustMedian; + + // Free temporary data buffers. + psFree(stats2); + + // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. + stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; + + return; + } + + // Determine how many data points fit inside this min/max range + // and are not masked, IF the maskVector is not NULL> + nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); + + // Allocate temporary vectors for the data. + unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + unsortedVector->n = unsortedVector->nalloc; + + sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + sortedVector->n = sortedVector->nalloc; + + // Determine if we must only use data points within a min/max range. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + + // Store all non-masked data points within the min/max range + // into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + unsortedVector->data.F32[count++] = maskVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } + } else { + // Store all non-masked data points into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + unsortedVector->data.F32[i] = myVector->data.F32[i]; + } + } + } + // Sort the temporary vectors. + psVectorSort(sortedVector, unsortedVector); + + // Calculate the median exactly. + if (0 == (nValues % 2)) { + stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] + + sortedVector->data.F32[nValues/2]); + } else { + stats->sampleMedian = sortedVector->data.F32[nValues/2]; + } + + // Free the temporary data structures. + psFree(unsortedVector); + psFree(sortedVector); +} + +/****************************************************************************** + This routine smoothes the data in the input robustHistogram with a + Gaussian of width sigma. + *****************************************************************************/ +psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram, + float sigma) +{ + int i = 0; // Loop index variable + int j = 0; // Loop index variable + float denom = 0.0; // Temporary variable + float expo = 0.0; // Temporary variable + float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients + psVector *smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32); + + for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) { + if (fabs(sigma) >= FLT_EPSILON) { + // If sigma does not equal zero, then we use Gaussian smoothing. + #ifdef DARWIN + denom = (float) sqrt(2.0 * M_PI * sigma * sigma); + #else + + denom = sqrtf(2.0 * M_PI * sigma * sigma); + #endif + + expo = - (float) ((i-GAUSS_WIDTH) * (i-GAUSS_WIDTH)); + expo/= (2.0 * sigma * sigma); + gaussianCoefs[i] = exp(expo/denom); + + // NOTE: Gaussian smoothing just isn't working with low sigma + // values. The problem is that the Gaussian coefficients are + // all zero, except for the middle coefficient, which is exactly + // one. Therefore, I'm using boxcar smoothing. + gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float) GAUSS_WIDTH)); + // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); + } else { + /* If sigma equals zero (all pixels have the same value) + * the above code will divide by zero. Therefore, we don't need + * to smooth the data. + */ + for(i=0;inums->n;i++) { + smooth->data.F32[i] = (float) robustHistogram->nums->data.S32[i]; + } + return(smooth); + } + } + + for(i=0;inums->n;i++) { + smooth->data.F32[i] = 0.0; + for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) { + if (((j+i) >= 0) && ((j+i) < smooth->n)) { + smooth->data.F32[i]+= (gaussianCoefs[j+GAUSS_WIDTH] * + (float) robustHistogram->nums->data.S32[j+i]); + } + } + } + return(smooth); +} + +/****************************************************************************** +p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates +the upper and/or lower quartiles of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorSampleQuartiles(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + psVector *unsortedVector = NULL; // Temporary vector + psVector *sortedVector = NULL; // Temporary vector + int i = 0; // Loop index variable + int count = 0; // # of points in this mean? + int nValues = 0; // # data points + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + // Determine if the number of data points exceed a threshold which will + // cause to generate robust stats, as opposed to exact stats. + if (myVector->n > stats->sampleLimit) { + psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); + psStats *stats2 = NULL; + // Calculate the robust quartiles. + stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE); + p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); + + // Store the robust quartiles into the sample quartile members. + stats->sampleUQ = stats2->robustUQ; + stats->sampleLQ = stats2->robustLQ; + + // Free temporary data buffers. + psFree(stats2); + + // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. + stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; + + return; + } + + // Determine how many data points fit inside this min/max range + // and are not maxed, IF the maskVector is not NULL> + nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); + + // Allocate temporary vectors for the data. + unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + unsortedVector->n = unsortedVector->nalloc; + sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + sortedVector->n = sortedVector->nalloc; + + // Determine if we must only use data points within a min/max range. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + // Store all non-masked data points within the min/max range + // into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } + } else { + // Store all non-masked data points into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + unsortedVector->data.F32[i] = myVector->data.F32[i]; + } + } + } + + // Sort the temporary vectors. + psVectorSort(sortedVector, unsortedVector); + + // Calculate the quartile points exactly. + stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)]; + stats->sampleLQ = sortedVector->data.F32[nValues / 4]; + + // Free the temporary data structures. + psFree(unsortedVector); + psFree(sortedVector); + // NOTE: This is the +} + + +/****************************************************************************** +p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the +stdev of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + +NOTE: the mean is always calculated exactly. Robust means are never +calculated in this routine. + *****************************************************************************/ +void p_psVectorSampleStdev(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + int countInt = 0; // # of data points being used + float countFloat = 0.0; // # of data points being used + float mean = 0.0; // The mean + float diff = 0.0; // Used in calculating stdev + float sumSquares = 0.0; // temporary variable + float sumDiffs = 0.0; // temporary variable + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + // This procedure requires the mean. If it has not been already + // calculated, then call p_psVectorSampleMean() + if (0 != isnan(stats->sampleMean)) { + p_psVectorSampleMean(myVector, maskVector, maskVal, stats); + } + mean = stats->sampleMean; + + if (stats->options & PS_STAT_USE_RANGE) { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + } + countInt = myVector->n; + } + } else { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + } + } else { + for (i=0;in;i++) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + countInt = myVector->n; + } + } + countFloat = (float) countInt; + + #ifdef DARWIN + + stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs * + sumDiffs/countFloat))/ (countFloat-1)); + #else + + stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs * + sumDiffs/countFloat))/ (countFloat-1)); + #endif +} + +/****************************************************************************** +p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the +clipped stats (mean or stdev) of the input vector. + +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorClippedStats(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + int j = 0; // Loop index variable + float clippedMean = 0.0; // self-explanatory + float clippedStdev = 0.0; // self-explanatory + float oldStanMean = 0.0; // Temporary variable + float oldStanStdev = 0.0; // Temporary variable + psVector *tmpMask = NULL; // Temporary vector + + // Endure that stats->clipIter is within the proper range. + if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) && + (stats->clipIter <= CLIPPED_NUM_ITER_UB))) { + psAbort(__func__, "Unallowed value for clipIter (%d).\n", + stats->clipIter); + } + + // Endure that stats->clipSigma is within the proper range. + if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && + (stats->clipSigma <= CLIPPED_SIGMA_UB))) { + psAbort(__func__, "Unallowed value for clipSigma (%f).\n", + stats->clipSigma); + } + + // 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 = myVector->n; + + // If we were called with a mask vector, then initialize the temporary + // mask vector with those values. + if (maskVector != NULL) { + for (i=0;in;i++) { + tmpMask->data.U8[i] = maskVector->data.U8[i]; + } + } + + // 1. Compute the sample median. + // NOTE: This seems odd. Verify with IfA that we want to calculate the + // median here, not the mean. + p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); + + // 2. Compute the sample standard deviation. + p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); + + // 3. Use the sample median as the first estimator of the mean X. + clippedMean = stats->sampleMean; + + // 4. Use the sample stdev as the first estimator of the mean stdev. + clippedStdev = stats->sampleStdev; + + // Must save the old sampleMean and sampleStdev since the following code + // block overwrites them. + oldStanMean = stats->sampleMean; + oldStanStdev = stats->sampleStdev; + + // 5. Repeat N times: + for (i=0;iclipIter;i++) { + for (j=0;jn;j++) { + // a) Exclude all values x_i for which |x_i - x| > K * stdev + if (fabs(myVector->data.F32[j] - clippedMean) > + (stats->clipSigma * clippedStdev)) { + tmpMask->data.U8[i] = 0xff; + } + // b) compute new mean and stdev + p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); + p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); + + // c) Use the new mean for x + clippedMean = stats->sampleMean; + + // d) Use the new stdev for stdev + clippedStdev = stats->sampleStdev; + } + + } + stats->sampleMean = oldStanMean; + stats->sampleStdev= oldStanStdev; + + // 7. The last calcuated value of x is the cliped mean. + if (stats->options & PS_STAT_CLIPPED_MEAN) { + stats->clippedMean = clippedMean; + } + + // 8. The last calcuated value of stdev is the cliped stdev. + if (stats->options & PS_STAT_CLIPPED_STDEV) { + stats->clippedStdev = clippedStdev; + } + + psFree(tmpMask); +} + +/***************************************************************************** + *****************************************************************************/ +void p_psNormalizeVector(psVector *myData) +{ + float min = (float) HUGE; + float max = (float) -HUGE; + float range = 0.0; + int i = 0; + + for (i=0;in;i++) { + if (myData->data.F32[i] < min) { + min = myData->data.F32[i]; + } + if (myData->data.F32[i] > max) { + max = myData->data.F32[i]; + } + } + + range = max - min; + for (i=0;in;i++) { + myData->data.F32[i] = (myData->data.F32[i] - min) / range; + } +} + + +float p_psGaussian(const psVector *restrict myData, + const psVector *restrict myParams) +{ + float x = myData->data.F32[0]; + float mean = myParams->data.F32[0]; + float stdev = myParams->data.F32[1]; + float tmp = exp(-((x-mean) * (x-mean)) / (2.0 * stdev * stdev)); + tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev))); + + printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); + return(tmp); +} + +float p_psGaussianDeriv(const psVector *restrict myData, + const psVector *restrict myParams, + int whichParam) +{ + float x = myData->data.F32[0]; + float mean = myParams->data.F32[0]; + float stdev = myParams->data.F32[1]; + float tmp = 0.0; + + if (whichParam == 0) { + // Return the derivative w.r.t. the mean. + tmp = (x - mean) * p_psGaussian(myData, myParams); + tmp/= (stdev * stdev); + } else if (whichParam == 1) { + // Return the derivative w.r.t. the stdev. + tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams); + tmp/= (stdev * stdev * stdev); + } + printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp); + + return(tmp); +} + + +float p_psQuadratic(const psVector *restrict myParams, + const psVector *restrict myCoords) +{ + float x = myCoords->data.F32[0]; + float A = myParams->data.F32[0]; + float B = myParams->data.F32[1]; + float C = myParams->data.F32[2]; + float tmp = 0.0; + + tmp = (A * x * x) + (B * x) + C; + return(tmp); +} + +float p_psQuadraticDeriv(const psVector *restrict myParams, + const psVector *restrict myCoords, + int whichParamDeriv) +{ + float x = myCoords->data.F32[0]; + float tmp = 0.0; + + if (whichParamDeriv == 0) { + tmp = x * x; + } else if (whichParamDeriv == 1) { + tmp = x; + } else if (whichParamDeriv == 2) { + tmp = 1.0; + } + + return(tmp); +} + +/****************************************************************************** +p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes +as input a 1-D polynomial of arbitrary order (though we are using 2nd-order +polynomials here) and a range of x-values for which it is defined: +[rangeLow, rangeHigh]. It determines the x-value of that polynomial such +that f(x) == midpoint. This functions uses a binary-search algorithm on the +range and assumes that the polnomial is monotonically increasing within that +range. + *****************************************************************************/ +float p_ps1DPolyMedian(psPolynomial1D *myPoly, + float rangeLow, + float rangeHigh, + float getThisValue) +{ + int numIterations=0; + float midpoint = 0.0; + float oldMidpoint = 1.0; + float f = 0.0; + + // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue); + + while (numIterations < MAX_ITERATIONS) { + midpoint = (rangeHigh + rangeLow) / 2.0; + if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) { + return(midpoint); + } + oldMidpoint = midpoint; + + f = psEvalPolynomial1D(midpoint, myPoly); + // printf("p_ps1DPolyMedian() iteration %d. f(%f) is %f\n", numIterations, midpoint, f); + if (fabs(f - getThisValue) <= FLT_EPSILON) { + return(midpoint); + } + + if (f > getThisValue) { + rangeHigh = midpoint; + } else { + rangeLow = midpoint; + } + numIterations++; + } + return(midpoint); +} + +/****************************************************************************** +p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes +as input a 1-D polynomial of arbitrary order (though we are using 2nd-order +polynomials here) and a range of x-values for which it is defined: +[rangeLow, rangeHigh]. It determines the x-value of that polynomial such +that f(x) == midpoint. This functions uses a binary-search algorithm on the +range and assumes that the polnomial is monotonically increasing within that +range. + *****************************************************************************/ +float p_psFitQuadratic(psHistogram *robustHistogram, + int binNum, + float fitFloat) +{ + /* + if ((binNum > 0) && + (binNum < (robustHistogram->nums->n+1))) { + x->data.F64[0] = (double) 0.5 * + (robustHistogram->bounds->data.F32[binNum-1] + + robustHistogram->bounds->data.F32[binNum]); + x->data.F64[1] = (double) 0.5 * + (robustHistogram->bounds->data.F32[binNum] + + robustHistogram->bounds->data.F32[binNum+1]); + x->data.F64[2] = (double) 0.5 * + (robustHistogram->bounds->data.F32[binNum+1] + + robustHistogram->bounds->data.F32[binNum+2]); + + y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1]; + y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum]; + y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1]; + + if (!((y->data.F64[0] <= fitFloat) && + (fitFloat <= y->data.F64[2]))) { + psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); + } + + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat); + } else { + return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] + + robustHistogram->bounds->data.F32[binNum])); + } + */ + return(0.0); +} + +#define PS_PRINT_VECTOR(NAME) \ +for (int my_i=0;my_in;my_i++) { \ + printf("%s->data.F32[%d] is %f\n", #NAME, my_i, NAME->data.F32[my_i]); \ +} \ +printf("\n"); \ + +/****************************************************************************** +p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure +calculates a variety of robust stat measures: +PS_STAT_ROBUST_MEAN +PS_STAT_ROBUST_MEDIAN +PS_STAT_ROBUST_MODE +PS_STAT_ROBUST_STDEV +PS_STAT_ROBUST_QUARTILE +I have included all that computation in a single function, as opposed to +breaking it across several functions for one primary reason: +they all +require the same basic initial processing steps (calculate the histogram, +etc.) +however there is no currently defined means, in the SDRS, to +specify this computation as a separate step. If the above robust stat +measures were calcualted in separate functiosn, then much of the initial +processing would be duplicated. +Inputs +myVector +maskVector +maskVal +stats +Returns +NULL +*****************************************************************************/ void p_psVectorRobustStats(const psVector *restrict myVector, const psVector *restrict maskVector, unsigned int maskVal, - psStats *stats); -#endif + psStats *stats) +{ + psHistogram *robustHistogram = NULL; + psVector *robustHistogramVector = NULL; + float binSize = 0.0; // Size of the histogram bins + int LQBinNum = -1; // Bin num for lower quartile + int UQBinNum = -1; // Bin num for upper quartile + int i = 0; // Loop index variable + int maxBinNum = 0; + float maxBinCount = 0.0; + float dL = 0.0; + int numBins = 0; + psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN); + // psImage *domain; + // psVector *errors; + // psVector *data; + // psVector *initialGuess; + // psVector *theParams; + // float chiSq=0.0; + // float max = -HUGE; + // float max_pos; + float myMean = 0.0; + float myStdev = 0.0; + float countFloat = 0.0; + float diff = 0.0; + float sumSquares = 0.0; + float sumDiffs = 0.0; + psVector *x = psVectorAlloc(3, PS_TYPE_F64); + psVector *y = psVectorAlloc(3, PS_TYPE_F64); + psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); + psPolynomial1D *myPoly = psPolynomial1DAlloc(2); + psVector *cumulativeRobustSumsFullRange = NULL; + psVector *cumulativeRobustSumsDlRange = NULL; + float sumRobust = 0.0; + float sumN50 = 0.0; + float sumNfit = 0.0; + float cumulativeMedian = 0.0; + + // Compute the initial bin size of the robust histogram. This is done + // by computing the clipped standard deviation of the vector, and dividing + // that by 10.0; + p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats); + binSize = tmpStats->clippedStdev / 10.0f; + + // If stats->clippedStdev == 0.0, then all data elements have the same + // value. Therefore, we can set the appropiate results and return. + if (fabs(binSize) <= FLT_EPSILON) { + if (stats->options & PS_STAT_ROBUST_MEAN) { + stats->robustMean = stats->clippedMean; + } + if (stats->options & PS_STAT_ROBUST_MEDIAN) { + stats->robustMedian = stats->clippedMean; + } + if (stats->options & PS_STAT_ROBUST_MODE) { + stats->robustMode = stats->clippedMean; + } + if (stats->options & PS_STAT_ROBUST_STDEV) { + stats->robustStdev = 0.0; + } + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + stats->robustUQ = stats->clippedMean; + stats->robustLQ = stats->clippedMean; + } + // XXX: Set these to the number of unmasked data points? + stats->robustNfit = 0.0; + stats->robustN50 = 0.0; + psFree(tmpStats); + return; + } + + // Determine minimum and maximum values in the data vector. + if (isnan(stats->min)) { + p_psVectorMin(myVector, maskVector, maskVal, stats); + } + if (isnan(stats->max)) { + p_psVectorMax(myVector, maskVector, maskVal, stats); + } + + // Create the histogram structure. NOTE: we can not specify the bin size + // precisely since the argument to psHistogramAlloc() is the number of + // bins, not the binSize. Also, if we get here, we know that + // binSize != 0.0. + numBins = (int) ((stats->max - stats->min) / binSize); + robustHistogram = psHistogramAlloc(stats->min, + stats->max, + numBins); + + // Populate the histogram array. + psHistogramVector(robustHistogram, myVector, maskVector, maskVal); + + // Smooth the histogram. + // XXX: is that the right stdev? + robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram, + tmpStats->clippedStdev/4.0f); + + // The following was necessary to fit a gaussian to the data, since + // gaussian functions produce data between 0.0 and 1.0. + // p_psNormalizeVector(robustHistogramVector); + + /************************************************************************** + Determine the lower/upper quartiles. + **************************************************************************/ + // We define a vector called "cumulativeRobustSums..." where the value at + // index position i is equal to the sum of bins 0:i. This will be used + // now and later in determining the lower/upper quartiles. + cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0]; + for (i=1;in;i++) { + cumulativeRobustSumsFullRange->data.F32[i] = + cumulativeRobustSumsFullRange->data.F32[i-1] + + robustHistogramVector->data.F32[i]; + } + sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1]; + + // Determine the bin number containing the lower quartile point. + LQBinNum = -1; + for (i=0;in;i++) { + if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust/4.0)) { + LQBinNum = i; + break; + } + } + + // Determine the bin number containing the upper quartile point. + UQBinNum = -1; + for (i=cumulativeRobustSumsFullRange->n-1;i>=0;i--) { + if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0*sumRobust/4.0)) { + UQBinNum = i; + break; + } + } + + if ((LQBinNum == -1) || + (UQBinNum == -1)) { + psAbort(__func__, "Could not determine the robust lower/upper quartiles."); + } + /************************************************************************** + Determine the mode in the range LQ:UQ. + **************************************************************************/ + // Determine the bin with the peak value in the range LQ to UQ. + maxBinNum = LQBinNum; + maxBinCount = robustHistogramVector->data.F32[LQBinNum]; + sumN50 = (float) robustHistogram->nums->data.S32[LQBinNum]; + for (i=LQBinNum+1;i<=UQBinNum;i++) { + if (robustHistogramVector->data.F32[i] > maxBinCount) { + maxBinNum = i; + maxBinCount = robustHistogramVector->data.F32[i]; + } + sumN50+= (float) robustHistogram->nums->data.S32[i]; + } + // XXX: is dL defined as the value at the LQ/UQ, or the bin number? + dL = (UQBinNum - LQBinNum) / 4; + + printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum); + + /************************************************************************** + Determine the mean/stdev for the bins in the range mode-dL to mode+dL + **************************************************************************/ + cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + for (i=0;in;i++) { + cumulativeRobustSumsDlRange->data.F32[i] = 0.0; + } + sumNfit = 0.0; + cumulativeMedian = 0.0; + for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) { + if ((0 <= i) && (i < robustHistogramVector->n)) { + cumulativeRobustSumsDlRange->data.F32[i] = + cumulativeRobustSumsDlRange->data.F32[i-1] + + robustHistogramVector->data.F32[i]; + cumulativeMedian+= robustHistogramVector->data.F32[i]; + sumNfit+= (float) robustHistogramVector->data.S32[i]; + } + } + + // Calculate the mean of the smoothed robust histogram in the range + // mode-dL to mode+dL. We use the midpoint of each bin as the mean for + // that bin (this is a non-exact approximation). + myMean = 0.0; + for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) { + if ((0 <= i) && (i < robustHistogramVector->n)) { + myMean+= (robustHistogramVector->data.F32[i]) * 0.5 * + (robustHistogram->bounds->data.F32[i+1] + + robustHistogram->bounds->data.F32[i]); + countFloat+= robustHistogramVector->data.F32[i]; + } + } + myMean/= countFloat; + + // Calculate the stdev of the smoothed robust histogram in the range + // mode-dL to mode+dL. We use the midpoint of each bin as the mean for + // that bin. + for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) { + if ((0 <= i) && (i < robustHistogramVector->n)) { + diff = (0.5 * (robustHistogram->bounds->data.F32[i+1] + + robustHistogram->bounds->data.F32[i])) - myMean; + sumSquares+= diff * diff * robustHistogramVector->data.F32[i]; + sumDiffs+= diff * robustHistogramVector->data.F32[i]; + } + } + myStdev = sqrt((sumSquares-(sumDiffs * sumDiffs/countFloat))/ (countFloat-1)); + + /************************************************************************** + Set the appropriate members in the output stats struct. + **************************************************************************/ + if (stats->options & PS_STAT_ROBUST_MEAN) { + stats->robustMean = myMean; + } + + if (stats->options & PS_STAT_ROBUST_MODE) { + stats->robustMode = 0.5 * + (robustHistogram->bounds->data.F32[maxBinNum] + + robustHistogram->bounds->data.F32[maxBinNum+1]); + } + + if (stats->options & PS_STAT_ROBUST_STDEV) { + stats->robustStdev = myStdev; + } + + // To determine the median (and later, the lower/upper quartile), we fit + // a quadratic to the three bins surrounding the bin containing the median. + // The quadratic y=f(x) with x being the midpoint of each bin, and y being + // the cumulative number of data points in all bins up to, and including, + // this bin. We then solve the quadratic for + + if (stats->options & PS_STAT_ROBUST_MEDIAN) { + if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n-1))) { + x->data.F64[0] = (double) 0.5 * + (robustHistogram->bounds->data.F32[maxBinNum-1] + + robustHistogram->bounds->data.F32[maxBinNum]); + x->data.F64[1] = (double) 0.5 * + (robustHistogram->bounds->data.F32[maxBinNum] + + robustHistogram->bounds->data.F32[maxBinNum+1]); + x->data.F64[2] = (double) 0.5 * + (robustHistogram->bounds->data.F32[maxBinNum+1] + + robustHistogram->bounds->data.F32[maxBinNum+2]); + + y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum-1]; + y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum]; + y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum+1]; + + // Ensure that cumulativeMedian/2 is actually within the range of the bins + // we are using. + cumulativeMedian*= 0.5; + if (!((y->data.F64[0] <= cumulativeMedian) && + (cumulativeMedian <= y->data.F64[2]))) { + printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]); + psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n"); + } + // XXX: yErr is not currently used by psGetArrayPolynomial(). We + // may have to set this meaningfully later. + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + // Determine the coefficients of the polynomial. + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + // Call p_ps1DPolyMedian(), which does a binary search on the + // polynomial, looking for the value x such that + // f(x) = cumulativeMedian. + stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], + x->data.F64[2], cumulativeMedian); + } else { + // If the mode is the first/last histogram bin, then simply use + // the midpoint of that bin. + stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum+1] + + robustHistogram->bounds->data.F32[maxBinNum]); + } + } + + // The lower/upper quartile calculations are very similar to the median + // calculations. We fit a quadratic to the array containing the + // cumulative data points in each bin, then search for x such that + // f(x) equals the lower/upper quartile exactly. + // + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1]; + + if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n-1))) { + x->data.F64[0] = (double) 0.5 * + (robustHistogram->bounds->data.F32[LQBinNum-1] + + robustHistogram->bounds->data.F32[LQBinNum]); + x->data.F64[1] = (double) 0.5 * + (robustHistogram->bounds->data.F32[LQBinNum] + + robustHistogram->bounds->data.F32[LQBinNum+1]); + x->data.F64[2] = (double) 0.5 * + (robustHistogram->bounds->data.F32[LQBinNum+1] + + robustHistogram->bounds->data.F32[LQBinNum+2]); + + y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum-1]; + y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum]; + y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum+1]; + + if (!((y->data.F64[0] <= (countFloat/4.0)) && + ((countFloat/4.0) <= y->data.F64[2]))) { + psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n"); + } + + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + stats->robustLQ = p_ps1DPolyMedian(myPoly, + x->data.F64[0], + x->data.F64[2], + countFloat/4.0); + + } else { + // If the LQ is the first/last histogram bin, then simply use + // the midpoint of that bin. + stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum+1] + + robustHistogram->bounds->data.F32[LQBinNum]); + } + + if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n-1))) { + x->data.F64[0] = (double) 0.5 * + (robustHistogram->bounds->data.F32[UQBinNum-1] + + robustHistogram->bounds->data.F32[UQBinNum]); + x->data.F64[1] = (double) 0.5 * + (robustHistogram->bounds->data.F32[UQBinNum] + + robustHistogram->bounds->data.F32[UQBinNum+1]); + x->data.F64[2] = (double) 0.5 * + (robustHistogram->bounds->data.F32[UQBinNum+1] + + robustHistogram->bounds->data.F32[UQBinNum+2]); + + y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum-1]; + y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum]; + y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum+1]; + + if (!((y->data.F64[0] <= (3.0 * countFloat/4.0)) && + ((3.0 * countFloat/4.0) <= y->data.F64[2]))) { + psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n"); + } + + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + stats->robustUQ = p_ps1DPolyMedian(myPoly, + x->data.F64[0], + x->data.F64[2], + 3.0*countFloat/4.0); + } else { + // If the UQ is the first/last histogram bin, then simply use + // the midpoint of that bin. + stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum+1] + + robustHistogram->bounds->data.F32[UQBinNum]); + } + } + stats->robustNfit = sumNfit; + stats->robustN50 = sumN50; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(tmpStats); + psFree(robustHistogram); + psFree(myPoly); + psFree(cumulativeRobustSumsFullRange); + psFree(cumulativeRobustSumsDlRange); +} + + + +/* +void p_ps_FitTheGaussian() +{ + + // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL + // NOTE: This code uses the psMinimize.c functions to perform + // the fit. It doesn't quite work, 100% of the time. For the time being + // I am commenting this code out, and replacing it with code which + // calculates the mean directly on the robustHistogram. + + domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32); + errors = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + data = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + initialGuess = psVectorAlloc(2, PS_TYPE_F32); + + max = -HUGE; + max_pos = -1; + for (i=0;in;i++) { + domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] + + robustHistogram->bounds->data.F32[i]); + data->data.F32[i] = (float) robustHistogramVector->data.F32[i]; + errors->data.F32[i] = 1.0; + //printf("DATA (%.2f, %.2f)\n", domain->data.F32[i][0], data->data.F32[i]); + if (data->data.F32[i] > max) { + max = data->data.F32[i]; + max_pos = domain->data.F32[i][0]; + } + } + + initialGuess->data.F32[0] = max_pos; + initialGuess->data.F32[1] = 1.0; + printf("Initial (mean. stdev) is (%.3f, %.3f)\n", initialGuess->data.F32[0], initialGuess->data.F32[1]); + + printf("Calling psMinimizeChi2()\n"); + theParams = psMinimizeChi2(p_psGaussian, + p_psGaussianDeriv, + domain, + data, + errors, + initialGuess, + NULL, + &chiSq); + printf("Called psMinimizeChi2()\n"); + printf("p_psVectorRobustStats() is (%f, %f)\n", theParams->data.F32[0], theParams->data.F32[1]); +} +*/ + + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ static void histogramFree(psHistogram *myHist); @@ -92,4 +1558,7 @@ return(NULL); } + if (n < 0) { + psAbort(__func__, "psHistogramAlloc() called with bin size %d.\n", n); + } // NOTE: Verify that this is the correct action. @@ -107,4 +1576,7 @@ // Calculate the bounds for each bin. binSize = (upper - lower) / (float) n; + // NOTE: Is the following necessary? It prevents the max data point + // from being in a non-existant bin. + binSize+= FLT_EPSILON; for (i=0;ibounds->data.F32[i] = lower + (binSize * (float) i); @@ -127,7 +1599,7 @@ /****************************************************************************** -psHistogramAlloc(lower, upper, n): allocate a non-uniform histogram structure -with the specifed bounds. The number of elements in the bounds vector is n. -Therefore, the number of bins is n-1. +psHistogramAllocGenric(bounds): allocate a non-uniform histogram structure +with the specifed bounds. + Input: bounds @@ -238,5 +1710,4 @@ numBins = out->nums->n; - for (i=0;in;i++) { // Check if this pixel is masked, and if so, skip it. @@ -256,8 +1727,16 @@ if (out->uniform == true) { binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; - binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) / binSize); + + // NOTE: This next if-statement really shouldn't be necessary. + // However, do to numerical lack of precision, we occasionally + // produce a binNum outside the range of bins. + if (binNum >= out->nums->n) { + binNum = out->nums->n-1; + } + (out->nums->data.S32[binNum])++; + // If this is a non-uniform histogram, determining the correct // bin number requires a bit more work. @@ -277,925 +1756,4 @@ return(out); } - -/***************************************************************************** -p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal -function that simply prints a vector to STDOUT. Used primarily for -debugging. - *****************************************************************************/ - -void p_psVectorPrint(psVector *myVector, - psVector *maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable. - - for (i=0;in;i++) { - if (maskVector != NULL) - printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]); - else - printf("Element %d is %f\n", i, myVector->data.F32[i]); - } -} - - -/****************************************************************************** - ****************************************************************************** - ****************************************************************************** - MISC PRIVATE STATISTICAL FUNCTIONS - - NOTE: it is assumed that any call to these statistical functions will - have been preceded by a call to the psVectorStats() function. Various - sanity tests will only be performed in psVectorStats(). - Is the mask vector the same length as the data vector? - Is the mask vector of type PS_TYPE_U8? - Is the stats data structure NULL? - Is the in data structure NULL? - Is the in data structure of type PS_TYPE_F32? - ****************************************************************************** - ****************************************************************************** - *****************************************************************************/ - - -/****************************************************************************** -p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the -mean of the input vector. -Inputs -Returns - NULL -ASSUMPTION: the mean is always calculated exactly. Robust means are never -calculated in this routine. - *****************************************************************************/ -void p_psVectorSampleMean(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - float mean = 0.0; // The mean - int count = 0; // # of points in this mean? - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - // If PS_STAT_USE_RANGE is requested, then we enter a slightly different - // loop. - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - // Check if the data is with the specified range - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - mean+= myVector->data.F32[i]; - count++; - } - } - mean/= (float) count; - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - mean+= myVector->data.F32[i]; - count++; - } - } - mean/= (float) count; - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - mean+= myVector->data.F32[i]; - count++; - } - } - mean/= (float) count; - } else { - for (i=0;in;i++) { - mean+= myVector->data.F32[i]; - } - mean/= (float) myVector->n; - } - } - - stats->sampleMean = mean; -} - -/****************************************************************************** -p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the -max of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorMax(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - float max = -MYMAXFLOAT; // The calculated maximum - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - // If PS_STAT_USE_RANGE is requested, then we enter a different loop. - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if ((myVector->data.F32[i] > max) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - max = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if ((myVector->data.F32[i] > max) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - max = myVector->data.F32[i]; - } - } - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if (myVector->data.F32[i] > max) { - max = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if (myVector->data.F32[i] > max) { - max = myVector->data.F32[i]; - } - } - } - } - - stats->max = max; -} - -/****************************************************************************** -p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the -minimum of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorMin(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - float min = MYMAXFLOAT; // The calculated maximum - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if ((myVector->data.F32[i] < min) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - min = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if ((myVector->data.F32[i] < min) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - min = myVector->data.F32[i]; - } - } - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if (myVector->data.F32[i] < min) { - min = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if (myVector->data.F32[i] < min) { - min = myVector->data.F32[i]; - } - } - } - } - - stats->min = min; -} - -/****************************************************************************** -p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the -number of non-masked pixels in the vector that fall within the min/max -range, if given. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -int p_psVectorNValues(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - int numData = 0; // The number of data points - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - numData++; - } - } - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - numData++; - } - } - } - } else { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - numData++; - } - } - } else { - numData = myVector->n; - } - } - return(numData); -} - - - -/****************************************************************************** -p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the -median of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorSampleMedian(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - psVector *unsortedVector = NULL; // Temporary vector - psVector *sortedVector = NULL; // Temporary vector - int i = 0; // Loop index variable - int count = 0; // # of points in this mean? - int nValues = 0; // # of points in vector - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - psStats *stats2 = NULL; // Temporary stats structure - - - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleMedian = stats2->robustMedian; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } - - // Determine how many data points fit inside this min/max range - // and are not masked, IF the maskVector is not NULL> - nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); - - // Allocate temporary vectors for the data. - unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = unsortedVector->nalloc; - - sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - sortedVector->n = sortedVector->nalloc; - - // Determine if we must only use data points within a min/max range. - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - - // Store all non-masked data points within the min/max range - // into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = maskVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } - } else { - // Store all non-masked data points into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - unsortedVector->data.F32[i] = myVector->data.F32[i]; - } - } - } - // Sort the temporary vectors. - psVectorSort(sortedVector, unsortedVector); - - // Calculate the median exactly. - if (0 == (nValues % 2)) { - stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] + - sortedVector->data.F32[nValues/2]); - } else { - stats->sampleMedian = sortedVector->data.F32[nValues/2]; - } - - // Free the temporary data structures. - psFree(unsortedVector); - psFree(sortedVector); -} - -/****************************************************************************** - This routine smoothes the data in the input robustHistogram with a - Gaussian of width sigma. - *****************************************************************************/ -void p_psVectorsmoothHistGaussian(psHistogram *robustHistogram, - float sigma) -{ - int i = 0; // Loop index variable - int j = 0; // Loop index variable - float tmpf = 0.0; // Temporary variable - float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients - - for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) { - if (fabs(sigma) <= FLT_EPSILON) { - // If sigma does not equal zero, then we use Gaussian smoothing. - #ifdef DARWIN - tmpf = (float) sqrt(2.0f * M_PI * sigma * sigma); - #else - - tmpf = sqrtf(2.0f * M_PI * sigma * sigma); - #endif - - gaussianCoefs[i] = (float) exp( (-((float) (i-GAUSS_WIDTH)) * - ((float) (i-GAUSS_WIDTH))) / - (2.0f * sigma * sigma)) / tmpf; - } else { - /* If sigma equals zero (all pixels have the same value) - * the above code will divide by zero. Therefore, we don't need - * to smooth the data. - */ - return; - } - } - - for(i=0;inums->n;i++) { - for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) { - if (((j+i) >= 0) && ((j+i) < robustHistogram->nums->n)) { - robustHistogram->nums->data.S32[j+i]+= - (gaussianCoefs[j+GAUSS_WIDTH] * - (float) robustHistogram->nums->data.S32[j+i]); - } - } - } -} - -/****************************************************************************** -p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates -the upper and/or lower quartiles of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorSampleQuartiles(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - psVector *unsortedVector = NULL; // Temporary vector - psVector *sortedVector = NULL; // Temporary vector - int i = 0; // Loop index variable - int count = 0; // # of points in this mean? - int nValues = 0; // # data points - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - psStats *stats2 = NULL; - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleUQ = stats2->robustUQ; - stats->sampleLQ = stats2->robustLQ; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } - - // Determine how many data points fit inside this min/max range - // and are not maxed, IF the maskVector is not NULL> - nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); - - // Allocate temporary vectors for the data. - unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = unsortedVector->nalloc; - sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - sortedVector->n = sortedVector->nalloc; - - // Determine if we must only use data points within a min/max range. - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - // Store all non-masked data points within the min/max range - // into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } - } else { - // Store all non-masked data points into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - unsortedVector->data.F32[i] = myVector->data.F32[i]; - } - } - } - - // Sort the temporary vectors. - psVectorSort(sortedVector, unsortedVector); - - // Calculate the quartile points exactly. - stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)]; - stats->sampleLQ = sortedVector->data.F32[nValues / 4]; - - // Free the temporary data structures. - psFree(unsortedVector); - psFree(sortedVector); - // NOTE: This is the -} - - -/****************************************************************************** -p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the -stdev of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - -NOTE: the mean is always calculated exactly. Robust means are never -calculated in this routine. - *****************************************************************************/ -void p_psVectorSampleStdev(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - int countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - // This procedure requires the mean. If it has not been already - // calculated, then call p_psVectorSampleMean() - if (0 != isnan(stats->sampleMean)) { - p_psVectorSampleMean(myVector, maskVector, maskVal, stats); - } - mean = stats->sampleMean; - - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - countInt = myVector->n; - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - } else { - for (i=0;in;i++) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - countInt = myVector->n; - } - } - countFloat = (float) countInt; - - #ifdef DARWIN - - stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs * - sumDiffs/countFloat))/ (countFloat-1)); - #else - - stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs * - sumDiffs/countFloat))/ (countFloat-1)); - #endif -} - -/****************************************************************************** -p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the -clipped stats (mean or stdev) of the input vector. - -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorClippedStats(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - int j = 0; // Loop index variable - float clippedMean = 0.0; // self-explanatory - float clippedStdev = 0.0; // self-explanatory - float oldStanMean = 0.0; // Temporary variable - float oldStanStdev = 0.0; // Temporary variable - psVector *tmpMask = NULL; // Temporary vector - - // Endure that stats->clipIter is within the proper range. - if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) && - (stats->clipIter <= CLIPPED_NUM_ITER_UB))) { - psAbort(__func__, "Unallowed value for clipIter (%d).\n", - stats->clipIter); - } - - // Endure that stats->clipSigma is within the proper range. - if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && - (stats->clipSigma <= CLIPPED_SIGMA_UB))) { - psAbort(__func__, "Unallowed value for clipSigma (%f).\n", - stats->clipSigma); - } - - // 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 = myVector->n; - - // If we were called with a mask vector, then initialize the temporary - // mask vector with those values. - if (maskVector != NULL) { - for (i=0;in;i++) { - tmpMask->data.U8[i] = maskVector->data.U8[i]; - } - } - - // 1. Compute the sample median. - // NOTE: This seems odd. Verify with IfA that we want to calculate the - // median here, not the mean. - p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); - - // 2. Compute the sample standard deviation. - p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); - - // 3. Use the sample median as the first estimator of the mean X. - clippedMean = stats->sampleMean; - - // 4. Use the sample stdev as the first estimator of the mean stdev. - clippedStdev = stats->sampleStdev; - - // Must save the old sampleMean and sampleStdev since the following code - // block overwrites them. - oldStanMean = stats->sampleMean; - oldStanStdev = stats->sampleStdev; - - // 5. Repeat N times: - for (i=0;iclipIter;i++) { - for (j=0;jn;j++) { - // a) Exclude all values x_i for which |x_i - x| > K * stdev - if (fabs(myVector->data.F32[j] - clippedMean) > - (stats->clipSigma * clippedStdev)) { - tmpMask->data.U8[i] = 0xff; - } - // b) compute new mean and stdev - p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); - p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); - - // c) Use the new mean for x - clippedMean = stats->sampleMean; - - // d) Use the new stdev for stdev - clippedStdev = stats->sampleStdev; - } - } - stats->sampleMean = oldStanMean; - stats->sampleStdev= oldStanStdev; - - // 7. The last calcuated value of x is the cliped mean. - if (stats->options & PS_STAT_CLIPPED_MEAN) { - stats->clippedMean = clippedMean; - } - - // 8. The last calcuated value of stdev is the cliped stdev. - if (stats->options & PS_STAT_CLIPPED_STDEV) { - stats->clippedStdev = clippedStdev; - } - - psFree(tmpMask); -} - - -/****************************************************************************** -p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure -calculates a variety of robust stat measures: - PS_STAT_ROBUST_MEAN - PS_STAT_ROBUST_MEDIAN - PS_STAT_ROBUST_MODE - PS_STAT_ROBUST_STDEV - PS_STAT_ROBUST_QUARTILE - I have included all that computation in a single function, as opposed to - breaking it across several functions for one primary reason: they all - require the same basic initial processing steps (calculate the histogram, - etc.) however there is no currently defined means, in the SDRS, to - specify this computation as a separate step. If the above robust stat - measures were calcualted in separate functiosn, then much of the initial - processing would be duplicated. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorRobustStats(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - psHistogram *robustHistogram = NULL; - float binSize = 0.0; // Size of the histogram bins - float sigmaE = 0.0; - int LQBinNum = -1; // Bin num for lower quartile - int UQBinNum = -1; // Bin num for upper quartile - int i = 0; // Loop index variable - int maxBinNum = 0; - int maxBinCount = 0; - float dL = 0.0; - int numBins = 0; - psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN); - - // NOTE: The SDRS states that the sample quartiles must be used to - // determine the initial bin sizes. However, the sample quartiles are - // calculated based on a full sort of the data set, regardless of the - // size of the data set. We should consult with IfA to ensure that this - // is really required. - /* - if (isnan(stats->sampleUQ) || - isnan(stats->sampleLQ)) { - stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE; - p_psVectorSampleQuartiles(myVector, - maskVector, - maskVal, - stats); - } - */ - // Compute the initial bin size of the robust histogram. - p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats); - binSize = stats->clippedStdev / 10.0f; - - // If stats->clippedStdev == 0.0, then all data elements have the same - // value. Therefore, we can set the appropiate results and return. - if (fabs(binSize) <= FLT_EPSILON) { - if (stats->options & PS_STAT_ROBUST_MEAN) { - stats->robustMean = stats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - stats->robustMedian = stats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = stats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_STDEV) { - stats->robustStdev = 0.0; - } - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - stats->robustUQ = stats->clippedMean; - stats->robustLQ = stats->clippedMean; - } - psFree(tmpStats); - psFree(robustHistogram); - return; - } - - // Detemine minimum and maximum values in the data vector. - if (isnan(stats->min)) { - p_psVectorMin(myVector, maskVector, maskVal, stats); - } - if (isnan(stats->max)) { - p_psVectorMax(myVector, maskVector, maskVal, stats); - } - - // Create the histogram structure (yes, 2 is necessary, not 1). Also, - // if we get here, we know that binSize != 0.0. - numBins = 2 + (int) ((stats->max - stats->min) / binSize); - - robustHistogram = psHistogramAlloc(stats->min, - stats->max, - numBins); - // Populate the histogram array. - psHistogramVector(robustHistogram, myVector, maskVector, maskVal); - - // Smooth the histogram. - p_psVectorsmoothHistGaussian(robustHistogram, sigmaE/4.0f); - - LQBinNum = -1; - UQBinNum = -1; - for (i=0;inums->n;i++) { - if ((robustHistogram->nums->data.S32[i] <= stats->sampleLQ) && - (stats->sampleLQ <= robustHistogram->nums->data.S32[i])) { - LQBinNum = i; - } - - if ((robustHistogram->nums->data.S32[i] <= stats->sampleUQ) && - (stats->sampleUQ <= robustHistogram->nums->data.S32[i])) { - UQBinNum = i; - } - } - - // Determine the bin with the peak value in the range LQ to UQ. - maxBinNum = LQBinNum; - maxBinCount = robustHistogram->nums->data.S32[maxBinNum]; - for (i=LQBinNum;i<=UQBinNum;i++) { - if (robustHistogram->nums->data.S32[i] > maxBinCount) { - maxBinNum = i; - maxBinCount = robustHistogram->nums->data.S32[i]; - } - } - - dL = (stats->robustUQ - stats->robustLQ) / 8.0; - - // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL - // NOTE: This step is dependent on the functions in psMinimize.c being - // implemented. Currently, they are not. - - if (stats->options & PS_STAT_ROBUST_MEAN) { - stats->robustMean = 0.0; - } - - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - stats->robustMedian = 0.0; - } - if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = maxBinNum; - } - if (stats->options & PS_STAT_ROBUST_STDEV) { - stats->robustStdev = 0.0; - } - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - stats->robustUQ = 0.0; - stats->robustLQ = 0.0; - } - stats->robustNfit = 0.0; - stats->robustN50 = 0.0; - psFree(tmpStats); - psFree(robustHistogram); -} - /****************************************************************************** @@ -1229,22 +1787,10 @@ } - // Ensure that the data is of type PS_TYPE_F32. Eventually, more data - // types will be implemented. - if (in->type.type != PS_TYPE_F32) { - psAbort(__func__, - "Only data type PS_TYPE_F32 is currently supported (0x%x).", - in->type.type); - } - - // Ensure that the mask vector is of the proper size and type. + PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32); if (mask != NULL) { - if (in->n != mask->n) { - psAbort(__func__, - "Vector data and vector mask are of different sizes."); - } - if (mask->type.type != PS_TYPE_U8) { - psAbort(__func__, - "Vector mask must be type PS_TYPE_U8"); - } + PS_CHECK_NULL_VECTOR(mask); + PS_CHECK_EMPTY_VECTOR(mask); + PS_CHECK_VECTOR_SIZE_EQUAL(mask, in); + PS_CHECK_VECTOR_TYPE(mask, PS_TYPE_U8); } @@ -1281,4 +1827,6 @@ (stats->options & PS_STAT_ROBUST_QUARTILE)) { p_psVectorRobustStats(in, mask, maskVal, stats); + + printf("HMMMM stats->robustMode is %f\n", stats->robustMode); } Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 1276) +++ /trunk/psLib/src/math/psStats.c (revision 1277) @@ -5,6 +5,9 @@ #include #include - +/*****************************************************************************/ +/* INCLUDE FILES */ +/*****************************************************************************/ #include "psMemory.h" +#include "psImage.h" #include "psVector.h" #include "psTrace.h" @@ -12,23 +15,1486 @@ #include "psAbort.h" #include "psStats.h" - -#include "float.h" -#define DEFAULT_ROBUST_SIZE_THRESHOLD 30000 // Vectors that are large than this -// will use robust statistical methods. -#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. -#define CLIPPED_NUM_ITER_LB 1 -#define CLIPPED_NUM_ITER_UB 10 -#define CLIPPED_SIGMA_LB 1.0 -#define CLIPPED_SIGMA_UB 10.0 -#define true 1 -#define false 0 -#define MYMAXFLOAT HUGE - -#ifndef DOXYGEN +#include "psSort.h" +#include "psMinimize.h" +#include "psFunctions.h" +======= + + /*****************************************************************************/ + /* DEFINE STATEMENTS */ + /*****************************************************************************/ + #define DEFAULT_ROBUST_SIZE_THRESHOLD 30000 // Vectors that are large than this + // will use robust statistical methods. + #define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. + #define CLIPPED_NUM_ITER_LB 1 + #define CLIPPED_NUM_ITER_UB 10 + #define CLIPPED_SIGMA_LB 1.0 + #define CLIPPED_SIGMA_UB 10.0 + #define true 1 + #define false 0 + #define MYMAXFLOAT HUGE + #define MAX_ITERATIONS 20 + + #ifndef DOXYGEN + void p_psVectorRobustStats(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats); +#endif + +/** Preprocessor macro to generate error on an incorrect type */ +#define PS_CHECK_VECTOR_TYPE(NAME, TYPE) \ +if (NAME->type.type != TYPE) { \ + psError(__func__,"Invalid operation: %s has incorrect type.", #NAME); \ +} + +/** Preprocessor macro to generate error on a NULL vector */ +#define PS_CHECK_NULL_VECTOR(NAME) \ +if (NAME == NULL || NAME->data.V == NULL) { \ + psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \ +} + +/** Preprocessor macro to generate error for zero length vector */ +#define PS_CHECK_EMPTY_VECTOR(NAME) \ +if (NAME->n < 1) { \ + psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \ +} + +/** Preprocessor macro to generate error on differing size vectors */ +#define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2) \ +if (VEC1->n != VEC2->n) { \ + psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \ +} + +/*****************************************************************************/ +/* TYPE DEFINITIONS */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +/***************************************************************************** +p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal +function that simply prints a vector to STDOUT. Used primarily for +debugging. + *****************************************************************************/ + +void p_psVectorPrint(psVector *myVector, + psVector *maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable. + + for (i=0;in;i++) { + if (maskVector != NULL) + printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]); + else + printf("Element %d is %f\n", i, myVector->data.F32[i]); + } +} + +/****************************************************************************** + ****************************************************************************** + ****************************************************************************** + MISC PRIVATE STATISTICAL FUNCTIONS + + NOTE: it is assumed that any call to these statistical functions will + have been preceded by a call to the psVectorStats() function. Various + sanity tests will only be performed in psVectorStats(). + Is the mask vector the same length as the data vector? + Is the mask vector of type PS_TYPE_U8? + Is the stats data structure NULL? + Is the in data structure NULL? + Is the in data structure of type PS_TYPE_F32? + ****************************************************************************** + ****************************************************************************** + *****************************************************************************/ + + +/****************************************************************************** +p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the +mean of the input vector. +Inputs +Returns + NULL +ASSUMPTION: the mean is always calculated exactly. Robust means are never +calculated in this routine. + *****************************************************************************/ +void p_psVectorSampleMean(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + float mean = 0.0; // The mean + int count = 0; // # of points in this mean? + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + // If PS_STAT_USE_RANGE is requested, then we enter a slightly different + // loop. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + // Check if the data is with the specified range + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + mean+= myVector->data.F32[i]; + count++; + } + } + mean/= (float) count; + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + mean+= myVector->data.F32[i]; + count++; + } + } + mean/= (float) count; + } + } else { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + mean+= myVector->data.F32[i]; + count++; + } + } + mean/= (float) count; + } else { + for (i=0;in;i++) { + mean+= myVector->data.F32[i]; + } + mean/= (float) myVector->n; + } + } + + stats->sampleMean = mean; +} + +/****************************************************************************** +p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the +max of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorMax(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + float max = -MYMAXFLOAT; // The calculated maximum + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + // If PS_STAT_USE_RANGE is requested, then we enter a different loop. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + if ((myVector->data.F32[i] > max) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + max = myVector->data.F32[i]; + } + } + } + } else { + for (i=0;in;i++) { + if ((myVector->data.F32[i] > max) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + max = myVector->data.F32[i]; + } + } + } + } else { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + if (myVector->data.F32[i] > max) { + max = myVector->data.F32[i]; + } + } + } + } else { + for (i=0;in;i++) { + if (myVector->data.F32[i] > max) { + max = myVector->data.F32[i]; + } + } + } + } + + stats->max = max; +} + +/****************************************************************************** +p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the +minimum of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorMin(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + float min = MYMAXFLOAT; // The calculated maximum + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + if ((myVector->data.F32[i] < min) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + min = myVector->data.F32[i]; + } + } + } + } else { + for (i=0;in;i++) { + if ((myVector->data.F32[i] < min) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + min = myVector->data.F32[i]; + } + } + } + } else { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + if (myVector->data.F32[i] < min) { + min = myVector->data.F32[i]; + } + } + } + } else { + for (i=0;in;i++) { + if (myVector->data.F32[i] < min) { + min = myVector->data.F32[i]; + } + } + } + } + + stats->min = min; +} + +/****************************************************************************** +p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the +number of non-masked pixels in the vector that fall within the min/max +range, if given. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +int p_psVectorNValues(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + int numData = 0; // The number of data points + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + numData++; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + numData++; + } + } + } + } else { + rangeMin = stats->min; + rangeMax = stats->max; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + numData++; + } + } + } else { + numData = myVector->n; + } + } + return(numData); +} + + + +/****************************************************************************** +p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the +median of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorSampleMedian(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + psVector *unsortedVector = NULL; // Temporary vector + psVector *sortedVector = NULL; // Temporary vector + int i = 0; // Loop index variable + int count = 0; // # of points in this mean? + int nValues = 0; // # of points in vector + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + psStats *stats2 = NULL; // Temporary stats structure + + + // Determine if the number of data points exceed a threshold which will + // cause to generate robust stats, as opposed to exact stats. + + if (myVector->n > stats->sampleLimit) { + psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); + + // Calculate the robust quartiles. + stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); + + // Store the robust quartiles into the sample quartile members. + stats->sampleMedian = stats2->robustMedian; + + // Free temporary data buffers. + psFree(stats2); + + // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. + stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; + + return; + } + + // Determine how many data points fit inside this min/max range + // and are not masked, IF the maskVector is not NULL> + nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); + + // Allocate temporary vectors for the data. + unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + unsortedVector->n = unsortedVector->nalloc; + + sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + sortedVector->n = sortedVector->nalloc; + + // Determine if we must only use data points within a min/max range. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + + // Store all non-masked data points within the min/max range + // into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + unsortedVector->data.F32[count++] = maskVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } + } else { + // Store all non-masked data points into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + unsortedVector->data.F32[i] = myVector->data.F32[i]; + } + } + } + // Sort the temporary vectors. + psVectorSort(sortedVector, unsortedVector); + + // Calculate the median exactly. + if (0 == (nValues % 2)) { + stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] + + sortedVector->data.F32[nValues/2]); + } else { + stats->sampleMedian = sortedVector->data.F32[nValues/2]; + } + + // Free the temporary data structures. + psFree(unsortedVector); + psFree(sortedVector); +} + +/****************************************************************************** + This routine smoothes the data in the input robustHistogram with a + Gaussian of width sigma. + *****************************************************************************/ +psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram, + float sigma) +{ + int i = 0; // Loop index variable + int j = 0; // Loop index variable + float denom = 0.0; // Temporary variable + float expo = 0.0; // Temporary variable + float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients + psVector *smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32); + + for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) { + if (fabs(sigma) >= FLT_EPSILON) { + // If sigma does not equal zero, then we use Gaussian smoothing. + #ifdef DARWIN + denom = (float) sqrt(2.0 * M_PI * sigma * sigma); + #else + + denom = sqrtf(2.0 * M_PI * sigma * sigma); + #endif + + expo = - (float) ((i-GAUSS_WIDTH) * (i-GAUSS_WIDTH)); + expo/= (2.0 * sigma * sigma); + gaussianCoefs[i] = exp(expo/denom); + + // NOTE: Gaussian smoothing just isn't working with low sigma + // values. The problem is that the Gaussian coefficients are + // all zero, except for the middle coefficient, which is exactly + // one. Therefore, I'm using boxcar smoothing. + gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float) GAUSS_WIDTH)); + // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); + } else { + /* If sigma equals zero (all pixels have the same value) + * the above code will divide by zero. Therefore, we don't need + * to smooth the data. + */ + for(i=0;inums->n;i++) { + smooth->data.F32[i] = (float) robustHistogram->nums->data.S32[i]; + } + return(smooth); + } + } + + for(i=0;inums->n;i++) { + smooth->data.F32[i] = 0.0; + for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) { + if (((j+i) >= 0) && ((j+i) < smooth->n)) { + smooth->data.F32[i]+= (gaussianCoefs[j+GAUSS_WIDTH] * + (float) robustHistogram->nums->data.S32[j+i]); + } + } + } + return(smooth); +} + +/****************************************************************************** +p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates +the upper and/or lower quartiles of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorSampleQuartiles(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + psVector *unsortedVector = NULL; // Temporary vector + psVector *sortedVector = NULL; // Temporary vector + int i = 0; // Loop index variable + int count = 0; // # of points in this mean? + int nValues = 0; // # data points + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + // Determine if the number of data points exceed a threshold which will + // cause to generate robust stats, as opposed to exact stats. + if (myVector->n > stats->sampleLimit) { + psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); + psStats *stats2 = NULL; + // Calculate the robust quartiles. + stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE); + p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); + + // Store the robust quartiles into the sample quartile members. + stats->sampleUQ = stats2->robustUQ; + stats->sampleLQ = stats2->robustLQ; + + // Free temporary data buffers. + psFree(stats2); + + // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. + stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; + + return; + } + + // Determine how many data points fit inside this min/max range + // and are not maxed, IF the maskVector is not NULL> + nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); + + // Allocate temporary vectors for the data. + unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + unsortedVector->n = unsortedVector->nalloc; + sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); + sortedVector->n = sortedVector->nalloc; + + // Determine if we must only use data points within a min/max range. + if (stats->options & PS_STAT_USE_RANGE) { + rangeMin = stats->min; + rangeMax = stats->max; + // Store all non-masked data points within the min/max range + // into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } + } else { + // Store all non-masked data points into the temporary vectors. + count = 0; + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + unsortedVector->data.F32[count++] = myVector->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + unsortedVector->data.F32[i] = myVector->data.F32[i]; + } + } + } + + // Sort the temporary vectors. + psVectorSort(sortedVector, unsortedVector); + + // Calculate the quartile points exactly. + stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)]; + stats->sampleLQ = sortedVector->data.F32[nValues / 4]; + + // Free the temporary data structures. + psFree(unsortedVector); + psFree(sortedVector); + // NOTE: This is the +} + + +/****************************************************************************** +p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the +stdev of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + +NOTE: the mean is always calculated exactly. Robust means are never +calculated in this routine. + *****************************************************************************/ +void p_psVectorSampleStdev(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + int countInt = 0; // # of data points being used + float countFloat = 0.0; // # of data points being used + float mean = 0.0; // The mean + float diff = 0.0; // Used in calculating stdev + float sumSquares = 0.0; // temporary variable + float sumDiffs = 0.0; // temporary variable + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + // This procedure requires the mean. If it has not been already + // calculated, then call p_psVectorSampleMean() + if (0 != isnan(stats->sampleMean)) { + p_psVectorSampleMean(myVector, maskVector, maskVal, stats); + } + mean = stats->sampleMean; + + if (stats->options & PS_STAT_USE_RANGE) { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + } + countInt = myVector->n; + } + } else { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + } + } else { + for (i=0;in;i++) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + countInt = myVector->n; + } + } + countFloat = (float) countInt; + + #ifdef DARWIN + + stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs * + sumDiffs/countFloat))/ (countFloat-1)); + #else + + stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs * + sumDiffs/countFloat))/ (countFloat-1)); + #endif +} + +/****************************************************************************** +p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the +clipped stats (mean or stdev) of the input vector. + +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorClippedStats(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + int j = 0; // Loop index variable + float clippedMean = 0.0; // self-explanatory + float clippedStdev = 0.0; // self-explanatory + float oldStanMean = 0.0; // Temporary variable + float oldStanStdev = 0.0; // Temporary variable + psVector *tmpMask = NULL; // Temporary vector + + // Endure that stats->clipIter is within the proper range. + if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) && + (stats->clipIter <= CLIPPED_NUM_ITER_UB))) { + psAbort(__func__, "Unallowed value for clipIter (%d).\n", + stats->clipIter); + } + + // Endure that stats->clipSigma is within the proper range. + if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && + (stats->clipSigma <= CLIPPED_SIGMA_UB))) { + psAbort(__func__, "Unallowed value for clipSigma (%f).\n", + stats->clipSigma); + } + + // 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 = myVector->n; + + // If we were called with a mask vector, then initialize the temporary + // mask vector with those values. + if (maskVector != NULL) { + for (i=0;in;i++) { + tmpMask->data.U8[i] = maskVector->data.U8[i]; + } + } + + // 1. Compute the sample median. + // NOTE: This seems odd. Verify with IfA that we want to calculate the + // median here, not the mean. + p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); + + // 2. Compute the sample standard deviation. + p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); + + // 3. Use the sample median as the first estimator of the mean X. + clippedMean = stats->sampleMean; + + // 4. Use the sample stdev as the first estimator of the mean stdev. + clippedStdev = stats->sampleStdev; + + // Must save the old sampleMean and sampleStdev since the following code + // block overwrites them. + oldStanMean = stats->sampleMean; + oldStanStdev = stats->sampleStdev; + + // 5. Repeat N times: + for (i=0;iclipIter;i++) { + for (j=0;jn;j++) { + // a) Exclude all values x_i for which |x_i - x| > K * stdev + if (fabs(myVector->data.F32[j] - clippedMean) > + (stats->clipSigma * clippedStdev)) { + tmpMask->data.U8[i] = 0xff; + } + // b) compute new mean and stdev + p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); + p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); + + // c) Use the new mean for x + clippedMean = stats->sampleMean; + + // d) Use the new stdev for stdev + clippedStdev = stats->sampleStdev; + } + + } + stats->sampleMean = oldStanMean; + stats->sampleStdev= oldStanStdev; + + // 7. The last calcuated value of x is the cliped mean. + if (stats->options & PS_STAT_CLIPPED_MEAN) { + stats->clippedMean = clippedMean; + } + + // 8. The last calcuated value of stdev is the cliped stdev. + if (stats->options & PS_STAT_CLIPPED_STDEV) { + stats->clippedStdev = clippedStdev; + } + + psFree(tmpMask); +} + +/***************************************************************************** + *****************************************************************************/ +void p_psNormalizeVector(psVector *myData) +{ + float min = (float) HUGE; + float max = (float) -HUGE; + float range = 0.0; + int i = 0; + + for (i=0;in;i++) { + if (myData->data.F32[i] < min) { + min = myData->data.F32[i]; + } + if (myData->data.F32[i] > max) { + max = myData->data.F32[i]; + } + } + + range = max - min; + for (i=0;in;i++) { + myData->data.F32[i] = (myData->data.F32[i] - min) / range; + } +} + + +float p_psGaussian(const psVector *restrict myData, + const psVector *restrict myParams) +{ + float x = myData->data.F32[0]; + float mean = myParams->data.F32[0]; + float stdev = myParams->data.F32[1]; + float tmp = exp(-((x-mean) * (x-mean)) / (2.0 * stdev * stdev)); + tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev))); + + printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); + return(tmp); +} + +float p_psGaussianDeriv(const psVector *restrict myData, + const psVector *restrict myParams, + int whichParam) +{ + float x = myData->data.F32[0]; + float mean = myParams->data.F32[0]; + float stdev = myParams->data.F32[1]; + float tmp = 0.0; + + if (whichParam == 0) { + // Return the derivative w.r.t. the mean. + tmp = (x - mean) * p_psGaussian(myData, myParams); + tmp/= (stdev * stdev); + } else if (whichParam == 1) { + // Return the derivative w.r.t. the stdev. + tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams); + tmp/= (stdev * stdev * stdev); + } + printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp); + + return(tmp); +} + + +float p_psQuadratic(const psVector *restrict myParams, + const psVector *restrict myCoords) +{ + float x = myCoords->data.F32[0]; + float A = myParams->data.F32[0]; + float B = myParams->data.F32[1]; + float C = myParams->data.F32[2]; + float tmp = 0.0; + + tmp = (A * x * x) + (B * x) + C; + return(tmp); +} + +float p_psQuadraticDeriv(const psVector *restrict myParams, + const psVector *restrict myCoords, + int whichParamDeriv) +{ + float x = myCoords->data.F32[0]; + float tmp = 0.0; + + if (whichParamDeriv == 0) { + tmp = x * x; + } else if (whichParamDeriv == 1) { + tmp = x; + } else if (whichParamDeriv == 2) { + tmp = 1.0; + } + + return(tmp); +} + +/****************************************************************************** +p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes +as input a 1-D polynomial of arbitrary order (though we are using 2nd-order +polynomials here) and a range of x-values for which it is defined: +[rangeLow, rangeHigh]. It determines the x-value of that polynomial such +that f(x) == midpoint. This functions uses a binary-search algorithm on the +range and assumes that the polnomial is monotonically increasing within that +range. + *****************************************************************************/ +float p_ps1DPolyMedian(psPolynomial1D *myPoly, + float rangeLow, + float rangeHigh, + float getThisValue) +{ + int numIterations=0; + float midpoint = 0.0; + float oldMidpoint = 1.0; + float f = 0.0; + + // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue); + + while (numIterations < MAX_ITERATIONS) { + midpoint = (rangeHigh + rangeLow) / 2.0; + if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) { + return(midpoint); + } + oldMidpoint = midpoint; + + f = psEvalPolynomial1D(midpoint, myPoly); + // printf("p_ps1DPolyMedian() iteration %d. f(%f) is %f\n", numIterations, midpoint, f); + if (fabs(f - getThisValue) <= FLT_EPSILON) { + return(midpoint); + } + + if (f > getThisValue) { + rangeHigh = midpoint; + } else { + rangeLow = midpoint; + } + numIterations++; + } + return(midpoint); +} + +/****************************************************************************** +p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes +as input a 1-D polynomial of arbitrary order (though we are using 2nd-order +polynomials here) and a range of x-values for which it is defined: +[rangeLow, rangeHigh]. It determines the x-value of that polynomial such +that f(x) == midpoint. This functions uses a binary-search algorithm on the +range and assumes that the polnomial is monotonically increasing within that +range. + *****************************************************************************/ +float p_psFitQuadratic(psHistogram *robustHistogram, + int binNum, + float fitFloat) +{ + /* + if ((binNum > 0) && + (binNum < (robustHistogram->nums->n+1))) { + x->data.F64[0] = (double) 0.5 * + (robustHistogram->bounds->data.F32[binNum-1] + + robustHistogram->bounds->data.F32[binNum]); + x->data.F64[1] = (double) 0.5 * + (robustHistogram->bounds->data.F32[binNum] + + robustHistogram->bounds->data.F32[binNum+1]); + x->data.F64[2] = (double) 0.5 * + (robustHistogram->bounds->data.F32[binNum+1] + + robustHistogram->bounds->data.F32[binNum+2]); + + y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1]; + y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum]; + y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1]; + + if (!((y->data.F64[0] <= fitFloat) && + (fitFloat <= y->data.F64[2]))) { + psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); + } + + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat); + } else { + return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] + + robustHistogram->bounds->data.F32[binNum])); + } + */ + return(0.0); +} + +#define PS_PRINT_VECTOR(NAME) \ +for (int my_i=0;my_in;my_i++) { \ + printf("%s->data.F32[%d] is %f\n", #NAME, my_i, NAME->data.F32[my_i]); \ +} \ +printf("\n"); \ + +/****************************************************************************** +p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure +calculates a variety of robust stat measures: +PS_STAT_ROBUST_MEAN +PS_STAT_ROBUST_MEDIAN +PS_STAT_ROBUST_MODE +PS_STAT_ROBUST_STDEV +PS_STAT_ROBUST_QUARTILE +I have included all that computation in a single function, as opposed to +breaking it across several functions for one primary reason: +they all +require the same basic initial processing steps (calculate the histogram, +etc.) +however there is no currently defined means, in the SDRS, to +specify this computation as a separate step. If the above robust stat +measures were calcualted in separate functiosn, then much of the initial +processing would be duplicated. +Inputs +myVector +maskVector +maskVal +stats +Returns +NULL +*****************************************************************************/ void p_psVectorRobustStats(const psVector *restrict myVector, const psVector *restrict maskVector, unsigned int maskVal, - psStats *stats); -#endif + psStats *stats) +{ + psHistogram *robustHistogram = NULL; + psVector *robustHistogramVector = NULL; + float binSize = 0.0; // Size of the histogram bins + int LQBinNum = -1; // Bin num for lower quartile + int UQBinNum = -1; // Bin num for upper quartile + int i = 0; // Loop index variable + int maxBinNum = 0; + float maxBinCount = 0.0; + float dL = 0.0; + int numBins = 0; + psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN); + // psImage *domain; + // psVector *errors; + // psVector *data; + // psVector *initialGuess; + // psVector *theParams; + // float chiSq=0.0; + // float max = -HUGE; + // float max_pos; + float myMean = 0.0; + float myStdev = 0.0; + float countFloat = 0.0; + float diff = 0.0; + float sumSquares = 0.0; + float sumDiffs = 0.0; + psVector *x = psVectorAlloc(3, PS_TYPE_F64); + psVector *y = psVectorAlloc(3, PS_TYPE_F64); + psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); + psPolynomial1D *myPoly = psPolynomial1DAlloc(2); + psVector *cumulativeRobustSumsFullRange = NULL; + psVector *cumulativeRobustSumsDlRange = NULL; + float sumRobust = 0.0; + float sumN50 = 0.0; + float sumNfit = 0.0; + float cumulativeMedian = 0.0; + + // Compute the initial bin size of the robust histogram. This is done + // by computing the clipped standard deviation of the vector, and dividing + // that by 10.0; + p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats); + binSize = tmpStats->clippedStdev / 10.0f; + + // If stats->clippedStdev == 0.0, then all data elements have the same + // value. Therefore, we can set the appropiate results and return. + if (fabs(binSize) <= FLT_EPSILON) { + if (stats->options & PS_STAT_ROBUST_MEAN) { + stats->robustMean = stats->clippedMean; + } + if (stats->options & PS_STAT_ROBUST_MEDIAN) { + stats->robustMedian = stats->clippedMean; + } + if (stats->options & PS_STAT_ROBUST_MODE) { + stats->robustMode = stats->clippedMean; + } + if (stats->options & PS_STAT_ROBUST_STDEV) { + stats->robustStdev = 0.0; + } + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + stats->robustUQ = stats->clippedMean; + stats->robustLQ = stats->clippedMean; + } + // XXX: Set these to the number of unmasked data points? + stats->robustNfit = 0.0; + stats->robustN50 = 0.0; + psFree(tmpStats); + return; + } + + // Determine minimum and maximum values in the data vector. + if (isnan(stats->min)) { + p_psVectorMin(myVector, maskVector, maskVal, stats); + } + if (isnan(stats->max)) { + p_psVectorMax(myVector, maskVector, maskVal, stats); + } + + // Create the histogram structure. NOTE: we can not specify the bin size + // precisely since the argument to psHistogramAlloc() is the number of + // bins, not the binSize. Also, if we get here, we know that + // binSize != 0.0. + numBins = (int) ((stats->max - stats->min) / binSize); + robustHistogram = psHistogramAlloc(stats->min, + stats->max, + numBins); + + // Populate the histogram array. + psHistogramVector(robustHistogram, myVector, maskVector, maskVal); + + // Smooth the histogram. + // XXX: is that the right stdev? + robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram, + tmpStats->clippedStdev/4.0f); + + // The following was necessary to fit a gaussian to the data, since + // gaussian functions produce data between 0.0 and 1.0. + // p_psNormalizeVector(robustHistogramVector); + + /************************************************************************** + Determine the lower/upper quartiles. + **************************************************************************/ + // We define a vector called "cumulativeRobustSums..." where the value at + // index position i is equal to the sum of bins 0:i. This will be used + // now and later in determining the lower/upper quartiles. + cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0]; + for (i=1;in;i++) { + cumulativeRobustSumsFullRange->data.F32[i] = + cumulativeRobustSumsFullRange->data.F32[i-1] + + robustHistogramVector->data.F32[i]; + } + sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1]; + + // Determine the bin number containing the lower quartile point. + LQBinNum = -1; + for (i=0;in;i++) { + if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust/4.0)) { + LQBinNum = i; + break; + } + } + + // Determine the bin number containing the upper quartile point. + UQBinNum = -1; + for (i=cumulativeRobustSumsFullRange->n-1;i>=0;i--) { + if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0*sumRobust/4.0)) { + UQBinNum = i; + break; + } + } + + if ((LQBinNum == -1) || + (UQBinNum == -1)) { + psAbort(__func__, "Could not determine the robust lower/upper quartiles."); + } + /************************************************************************** + Determine the mode in the range LQ:UQ. + **************************************************************************/ + // Determine the bin with the peak value in the range LQ to UQ. + maxBinNum = LQBinNum; + maxBinCount = robustHistogramVector->data.F32[LQBinNum]; + sumN50 = (float) robustHistogram->nums->data.S32[LQBinNum]; + for (i=LQBinNum+1;i<=UQBinNum;i++) { + if (robustHistogramVector->data.F32[i] > maxBinCount) { + maxBinNum = i; + maxBinCount = robustHistogramVector->data.F32[i]; + } + sumN50+= (float) robustHistogram->nums->data.S32[i]; + } + // XXX: is dL defined as the value at the LQ/UQ, or the bin number? + dL = (UQBinNum - LQBinNum) / 4; + + printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum); + + /************************************************************************** + Determine the mean/stdev for the bins in the range mode-dL to mode+dL + **************************************************************************/ + cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + for (i=0;in;i++) { + cumulativeRobustSumsDlRange->data.F32[i] = 0.0; + } + sumNfit = 0.0; + cumulativeMedian = 0.0; + for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) { + if ((0 <= i) && (i < robustHistogramVector->n)) { + cumulativeRobustSumsDlRange->data.F32[i] = + cumulativeRobustSumsDlRange->data.F32[i-1] + + robustHistogramVector->data.F32[i]; + cumulativeMedian+= robustHistogramVector->data.F32[i]; + sumNfit+= (float) robustHistogramVector->data.S32[i]; + } + } + + // Calculate the mean of the smoothed robust histogram in the range + // mode-dL to mode+dL. We use the midpoint of each bin as the mean for + // that bin (this is a non-exact approximation). + myMean = 0.0; + for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) { + if ((0 <= i) && (i < robustHistogramVector->n)) { + myMean+= (robustHistogramVector->data.F32[i]) * 0.5 * + (robustHistogram->bounds->data.F32[i+1] + + robustHistogram->bounds->data.F32[i]); + countFloat+= robustHistogramVector->data.F32[i]; + } + } + myMean/= countFloat; + + // Calculate the stdev of the smoothed robust histogram in the range + // mode-dL to mode+dL. We use the midpoint of each bin as the mean for + // that bin. + for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) { + if ((0 <= i) && (i < robustHistogramVector->n)) { + diff = (0.5 * (robustHistogram->bounds->data.F32[i+1] + + robustHistogram->bounds->data.F32[i])) - myMean; + sumSquares+= diff * diff * robustHistogramVector->data.F32[i]; + sumDiffs+= diff * robustHistogramVector->data.F32[i]; + } + } + myStdev = sqrt((sumSquares-(sumDiffs * sumDiffs/countFloat))/ (countFloat-1)); + + /************************************************************************** + Set the appropriate members in the output stats struct. + **************************************************************************/ + if (stats->options & PS_STAT_ROBUST_MEAN) { + stats->robustMean = myMean; + } + + if (stats->options & PS_STAT_ROBUST_MODE) { + stats->robustMode = 0.5 * + (robustHistogram->bounds->data.F32[maxBinNum] + + robustHistogram->bounds->data.F32[maxBinNum+1]); + } + + if (stats->options & PS_STAT_ROBUST_STDEV) { + stats->robustStdev = myStdev; + } + + // To determine the median (and later, the lower/upper quartile), we fit + // a quadratic to the three bins surrounding the bin containing the median. + // The quadratic y=f(x) with x being the midpoint of each bin, and y being + // the cumulative number of data points in all bins up to, and including, + // this bin. We then solve the quadratic for + + if (stats->options & PS_STAT_ROBUST_MEDIAN) { + if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n-1))) { + x->data.F64[0] = (double) 0.5 * + (robustHistogram->bounds->data.F32[maxBinNum-1] + + robustHistogram->bounds->data.F32[maxBinNum]); + x->data.F64[1] = (double) 0.5 * + (robustHistogram->bounds->data.F32[maxBinNum] + + robustHistogram->bounds->data.F32[maxBinNum+1]); + x->data.F64[2] = (double) 0.5 * + (robustHistogram->bounds->data.F32[maxBinNum+1] + + robustHistogram->bounds->data.F32[maxBinNum+2]); + + y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum-1]; + y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum]; + y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum+1]; + + // Ensure that cumulativeMedian/2 is actually within the range of the bins + // we are using. + cumulativeMedian*= 0.5; + if (!((y->data.F64[0] <= cumulativeMedian) && + (cumulativeMedian <= y->data.F64[2]))) { + printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]); + psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n"); + } + // XXX: yErr is not currently used by psGetArrayPolynomial(). We + // may have to set this meaningfully later. + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + // Determine the coefficients of the polynomial. + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + // Call p_ps1DPolyMedian(), which does a binary search on the + // polynomial, looking for the value x such that + // f(x) = cumulativeMedian. + stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], + x->data.F64[2], cumulativeMedian); + } else { + // If the mode is the first/last histogram bin, then simply use + // the midpoint of that bin. + stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum+1] + + robustHistogram->bounds->data.F32[maxBinNum]); + } + } + + // The lower/upper quartile calculations are very similar to the median + // calculations. We fit a quadratic to the array containing the + // cumulative data points in each bin, then search for x such that + // f(x) equals the lower/upper quartile exactly. + // + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1]; + + if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n-1))) { + x->data.F64[0] = (double) 0.5 * + (robustHistogram->bounds->data.F32[LQBinNum-1] + + robustHistogram->bounds->data.F32[LQBinNum]); + x->data.F64[1] = (double) 0.5 * + (robustHistogram->bounds->data.F32[LQBinNum] + + robustHistogram->bounds->data.F32[LQBinNum+1]); + x->data.F64[2] = (double) 0.5 * + (robustHistogram->bounds->data.F32[LQBinNum+1] + + robustHistogram->bounds->data.F32[LQBinNum+2]); + + y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum-1]; + y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum]; + y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum+1]; + + if (!((y->data.F64[0] <= (countFloat/4.0)) && + ((countFloat/4.0) <= y->data.F64[2]))) { + psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n"); + } + + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + stats->robustLQ = p_ps1DPolyMedian(myPoly, + x->data.F64[0], + x->data.F64[2], + countFloat/4.0); + + } else { + // If the LQ is the first/last histogram bin, then simply use + // the midpoint of that bin. + stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum+1] + + robustHistogram->bounds->data.F32[LQBinNum]); + } + + if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n-1))) { + x->data.F64[0] = (double) 0.5 * + (robustHistogram->bounds->data.F32[UQBinNum-1] + + robustHistogram->bounds->data.F32[UQBinNum]); + x->data.F64[1] = (double) 0.5 * + (robustHistogram->bounds->data.F32[UQBinNum] + + robustHistogram->bounds->data.F32[UQBinNum+1]); + x->data.F64[2] = (double) 0.5 * + (robustHistogram->bounds->data.F32[UQBinNum+1] + + robustHistogram->bounds->data.F32[UQBinNum+2]); + + y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum-1]; + y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum]; + y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum+1]; + + if (!((y->data.F64[0] <= (3.0 * countFloat/4.0)) && + ((3.0 * countFloat/4.0) <= y->data.F64[2]))) { + psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n"); + } + + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + stats->robustUQ = p_ps1DPolyMedian(myPoly, + x->data.F64[0], + x->data.F64[2], + 3.0*countFloat/4.0); + } else { + // If the UQ is the first/last histogram bin, then simply use + // the midpoint of that bin. + stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum+1] + + robustHistogram->bounds->data.F32[UQBinNum]); + } + } + stats->robustNfit = sumNfit; + stats->robustN50 = sumN50; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(tmpStats); + psFree(robustHistogram); + psFree(myPoly); + psFree(cumulativeRobustSumsFullRange); + psFree(cumulativeRobustSumsDlRange); +} + + + +/* +void p_ps_FitTheGaussian() +{ + + // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL + // NOTE: This code uses the psMinimize.c functions to perform + // the fit. It doesn't quite work, 100% of the time. For the time being + // I am commenting this code out, and replacing it with code which + // calculates the mean directly on the robustHistogram. + + domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32); + errors = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + data = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + initialGuess = psVectorAlloc(2, PS_TYPE_F32); + + max = -HUGE; + max_pos = -1; + for (i=0;in;i++) { + domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] + + robustHistogram->bounds->data.F32[i]); + data->data.F32[i] = (float) robustHistogramVector->data.F32[i]; + errors->data.F32[i] = 1.0; + //printf("DATA (%.2f, %.2f)\n", domain->data.F32[i][0], data->data.F32[i]); + if (data->data.F32[i] > max) { + max = data->data.F32[i]; + max_pos = domain->data.F32[i][0]; + } + } + + initialGuess->data.F32[0] = max_pos; + initialGuess->data.F32[1] = 1.0; + printf("Initial (mean. stdev) is (%.3f, %.3f)\n", initialGuess->data.F32[0], initialGuess->data.F32[1]); + + printf("Calling psMinimizeChi2()\n"); + theParams = psMinimizeChi2(p_psGaussian, + p_psGaussianDeriv, + domain, + data, + errors, + initialGuess, + NULL, + &chiSq); + printf("Called psMinimizeChi2()\n"); + printf("p_psVectorRobustStats() is (%f, %f)\n", theParams->data.F32[0], theParams->data.F32[1]); +} +*/ + + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ static void histogramFree(psHistogram *myHist); @@ -92,4 +1558,7 @@ return(NULL); } + if (n < 0) { + psAbort(__func__, "psHistogramAlloc() called with bin size %d.\n", n); + } // NOTE: Verify that this is the correct action. @@ -107,4 +1576,7 @@ // Calculate the bounds for each bin. binSize = (upper - lower) / (float) n; + // NOTE: Is the following necessary? It prevents the max data point + // from being in a non-existant bin. + binSize+= FLT_EPSILON; for (i=0;ibounds->data.F32[i] = lower + (binSize * (float) i); @@ -127,7 +1599,7 @@ /****************************************************************************** -psHistogramAlloc(lower, upper, n): allocate a non-uniform histogram structure -with the specifed bounds. The number of elements in the bounds vector is n. -Therefore, the number of bins is n-1. +psHistogramAllocGenric(bounds): allocate a non-uniform histogram structure +with the specifed bounds. + Input: bounds @@ -238,5 +1710,4 @@ numBins = out->nums->n; - for (i=0;in;i++) { // Check if this pixel is masked, and if so, skip it. @@ -256,8 +1727,16 @@ if (out->uniform == true) { binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; - binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) / binSize); + + // NOTE: This next if-statement really shouldn't be necessary. + // However, do to numerical lack of precision, we occasionally + // produce a binNum outside the range of bins. + if (binNum >= out->nums->n) { + binNum = out->nums->n-1; + } + (out->nums->data.S32[binNum])++; + // If this is a non-uniform histogram, determining the correct // bin number requires a bit more work. @@ -277,925 +1756,4 @@ return(out); } - -/***************************************************************************** -p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal -function that simply prints a vector to STDOUT. Used primarily for -debugging. - *****************************************************************************/ - -void p_psVectorPrint(psVector *myVector, - psVector *maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable. - - for (i=0;in;i++) { - if (maskVector != NULL) - printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]); - else - printf("Element %d is %f\n", i, myVector->data.F32[i]); - } -} - - -/****************************************************************************** - ****************************************************************************** - ****************************************************************************** - MISC PRIVATE STATISTICAL FUNCTIONS - - NOTE: it is assumed that any call to these statistical functions will - have been preceded by a call to the psVectorStats() function. Various - sanity tests will only be performed in psVectorStats(). - Is the mask vector the same length as the data vector? - Is the mask vector of type PS_TYPE_U8? - Is the stats data structure NULL? - Is the in data structure NULL? - Is the in data structure of type PS_TYPE_F32? - ****************************************************************************** - ****************************************************************************** - *****************************************************************************/ - - -/****************************************************************************** -p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the -mean of the input vector. -Inputs -Returns - NULL -ASSUMPTION: the mean is always calculated exactly. Robust means are never -calculated in this routine. - *****************************************************************************/ -void p_psVectorSampleMean(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - float mean = 0.0; // The mean - int count = 0; // # of points in this mean? - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - // If PS_STAT_USE_RANGE is requested, then we enter a slightly different - // loop. - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - // Check if the data is with the specified range - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - mean+= myVector->data.F32[i]; - count++; - } - } - mean/= (float) count; - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - mean+= myVector->data.F32[i]; - count++; - } - } - mean/= (float) count; - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - mean+= myVector->data.F32[i]; - count++; - } - } - mean/= (float) count; - } else { - for (i=0;in;i++) { - mean+= myVector->data.F32[i]; - } - mean/= (float) myVector->n; - } - } - - stats->sampleMean = mean; -} - -/****************************************************************************** -p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the -max of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorMax(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - float max = -MYMAXFLOAT; // The calculated maximum - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - // If PS_STAT_USE_RANGE is requested, then we enter a different loop. - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if ((myVector->data.F32[i] > max) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - max = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if ((myVector->data.F32[i] > max) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - max = myVector->data.F32[i]; - } - } - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if (myVector->data.F32[i] > max) { - max = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if (myVector->data.F32[i] > max) { - max = myVector->data.F32[i]; - } - } - } - } - - stats->max = max; -} - -/****************************************************************************** -p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the -minimum of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorMin(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - float min = MYMAXFLOAT; // The calculated maximum - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if ((myVector->data.F32[i] < min) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - min = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if ((myVector->data.F32[i] < min) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - min = myVector->data.F32[i]; - } - } - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - if (myVector->data.F32[i] < min) { - min = myVector->data.F32[i]; - } - } - } - } else { - for (i=0;in;i++) { - if (myVector->data.F32[i] < min) { - min = myVector->data.F32[i]; - } - } - } - } - - stats->min = min; -} - -/****************************************************************************** -p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the -number of non-masked pixels in the vector that fall within the min/max -range, if given. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -int p_psVectorNValues(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - int numData = 0; // The number of data points - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - numData++; - } - } - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - numData++; - } - } - } - } else { - rangeMin = stats->min; - rangeMax = stats->max; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - numData++; - } - } - } else { - numData = myVector->n; - } - } - return(numData); -} - - - -/****************************************************************************** -p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the -median of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorSampleMedian(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - psVector *unsortedVector = NULL; // Temporary vector - psVector *sortedVector = NULL; // Temporary vector - int i = 0; // Loop index variable - int count = 0; // # of points in this mean? - int nValues = 0; // # of points in vector - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - psStats *stats2 = NULL; // Temporary stats structure - - - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleMedian = stats2->robustMedian; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } - - // Determine how many data points fit inside this min/max range - // and are not masked, IF the maskVector is not NULL> - nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); - - // Allocate temporary vectors for the data. - unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = unsortedVector->nalloc; - - sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - sortedVector->n = sortedVector->nalloc; - - // Determine if we must only use data points within a min/max range. - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - - // Store all non-masked data points within the min/max range - // into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = maskVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } - } else { - // Store all non-masked data points into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - unsortedVector->data.F32[i] = myVector->data.F32[i]; - } - } - } - // Sort the temporary vectors. - psVectorSort(sortedVector, unsortedVector); - - // Calculate the median exactly. - if (0 == (nValues % 2)) { - stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] + - sortedVector->data.F32[nValues/2]); - } else { - stats->sampleMedian = sortedVector->data.F32[nValues/2]; - } - - // Free the temporary data structures. - psFree(unsortedVector); - psFree(sortedVector); -} - -/****************************************************************************** - This routine smoothes the data in the input robustHistogram with a - Gaussian of width sigma. - *****************************************************************************/ -void p_psVectorsmoothHistGaussian(psHistogram *robustHistogram, - float sigma) -{ - int i = 0; // Loop index variable - int j = 0; // Loop index variable - float tmpf = 0.0; // Temporary variable - float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients - - for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) { - if (fabs(sigma) <= FLT_EPSILON) { - // If sigma does not equal zero, then we use Gaussian smoothing. - #ifdef DARWIN - tmpf = (float) sqrt(2.0f * M_PI * sigma * sigma); - #else - - tmpf = sqrtf(2.0f * M_PI * sigma * sigma); - #endif - - gaussianCoefs[i] = (float) exp( (-((float) (i-GAUSS_WIDTH)) * - ((float) (i-GAUSS_WIDTH))) / - (2.0f * sigma * sigma)) / tmpf; - } else { - /* If sigma equals zero (all pixels have the same value) - * the above code will divide by zero. Therefore, we don't need - * to smooth the data. - */ - return; - } - } - - for(i=0;inums->n;i++) { - for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) { - if (((j+i) >= 0) && ((j+i) < robustHistogram->nums->n)) { - robustHistogram->nums->data.S32[j+i]+= - (gaussianCoefs[j+GAUSS_WIDTH] * - (float) robustHistogram->nums->data.S32[j+i]); - } - } - } -} - -/****************************************************************************** -p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates -the upper and/or lower quartiles of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorSampleQuartiles(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - psVector *unsortedVector = NULL; // Temporary vector - psVector *sortedVector = NULL; // Temporary vector - int i = 0; // Loop index variable - int count = 0; // # of points in this mean? - int nValues = 0; // # data points - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - psStats *stats2 = NULL; - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleUQ = stats2->robustUQ; - stats->sampleLQ = stats2->robustLQ; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } - - // Determine how many data points fit inside this min/max range - // and are not maxed, IF the maskVector is not NULL> - nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); - - // Allocate temporary vectors for the data. - unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = unsortedVector->nalloc; - sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - sortedVector->n = sortedVector->nalloc; - - // Determine if we must only use data points within a min/max range. - if (stats->options & PS_STAT_USE_RANGE) { - rangeMin = stats->min; - rangeMax = stats->max; - // Store all non-masked data points within the min/max range - // into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } - } else { - // Store all non-masked data points into the temporary vectors. - count = 0; - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - unsortedVector->data.F32[count++] = myVector->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - unsortedVector->data.F32[i] = myVector->data.F32[i]; - } - } - } - - // Sort the temporary vectors. - psVectorSort(sortedVector, unsortedVector); - - // Calculate the quartile points exactly. - stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)]; - stats->sampleLQ = sortedVector->data.F32[nValues / 4]; - - // Free the temporary data structures. - psFree(unsortedVector); - psFree(sortedVector); - // NOTE: This is the -} - - -/****************************************************************************** -p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the -stdev of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - -NOTE: the mean is always calculated exactly. Robust means are never -calculated in this routine. - *****************************************************************************/ -void p_psVectorSampleStdev(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - int countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - // This procedure requires the mean. If it has not been already - // calculated, then call p_psVectorSampleMean() - if (0 != isnan(stats->sampleMean)) { - p_psVectorSampleMean(myVector, maskVector, maskVal, stats); - } - mean = stats->sampleMean; - - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - countInt = myVector->n; - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - } else { - for (i=0;in;i++) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - countInt = myVector->n; - } - } - countFloat = (float) countInt; - - #ifdef DARWIN - - stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs * - sumDiffs/countFloat))/ (countFloat-1)); - #else - - stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs * - sumDiffs/countFloat))/ (countFloat-1)); - #endif -} - -/****************************************************************************** -p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the -clipped stats (mean or stdev) of the input vector. - -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorClippedStats(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - int j = 0; // Loop index variable - float clippedMean = 0.0; // self-explanatory - float clippedStdev = 0.0; // self-explanatory - float oldStanMean = 0.0; // Temporary variable - float oldStanStdev = 0.0; // Temporary variable - psVector *tmpMask = NULL; // Temporary vector - - // Endure that stats->clipIter is within the proper range. - if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) && - (stats->clipIter <= CLIPPED_NUM_ITER_UB))) { - psAbort(__func__, "Unallowed value for clipIter (%d).\n", - stats->clipIter); - } - - // Endure that stats->clipSigma is within the proper range. - if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && - (stats->clipSigma <= CLIPPED_SIGMA_UB))) { - psAbort(__func__, "Unallowed value for clipSigma (%f).\n", - stats->clipSigma); - } - - // 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 = myVector->n; - - // If we were called with a mask vector, then initialize the temporary - // mask vector with those values. - if (maskVector != NULL) { - for (i=0;in;i++) { - tmpMask->data.U8[i] = maskVector->data.U8[i]; - } - } - - // 1. Compute the sample median. - // NOTE: This seems odd. Verify with IfA that we want to calculate the - // median here, not the mean. - p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); - - // 2. Compute the sample standard deviation. - p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); - - // 3. Use the sample median as the first estimator of the mean X. - clippedMean = stats->sampleMean; - - // 4. Use the sample stdev as the first estimator of the mean stdev. - clippedStdev = stats->sampleStdev; - - // Must save the old sampleMean and sampleStdev since the following code - // block overwrites them. - oldStanMean = stats->sampleMean; - oldStanStdev = stats->sampleStdev; - - // 5. Repeat N times: - for (i=0;iclipIter;i++) { - for (j=0;jn;j++) { - // a) Exclude all values x_i for which |x_i - x| > K * stdev - if (fabs(myVector->data.F32[j] - clippedMean) > - (stats->clipSigma * clippedStdev)) { - tmpMask->data.U8[i] = 0xff; - } - // b) compute new mean and stdev - p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); - p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); - - // c) Use the new mean for x - clippedMean = stats->sampleMean; - - // d) Use the new stdev for stdev - clippedStdev = stats->sampleStdev; - } - } - stats->sampleMean = oldStanMean; - stats->sampleStdev= oldStanStdev; - - // 7. The last calcuated value of x is the cliped mean. - if (stats->options & PS_STAT_CLIPPED_MEAN) { - stats->clippedMean = clippedMean; - } - - // 8. The last calcuated value of stdev is the cliped stdev. - if (stats->options & PS_STAT_CLIPPED_STDEV) { - stats->clippedStdev = clippedStdev; - } - - psFree(tmpMask); -} - - -/****************************************************************************** -p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure -calculates a variety of robust stat measures: - PS_STAT_ROBUST_MEAN - PS_STAT_ROBUST_MEDIAN - PS_STAT_ROBUST_MODE - PS_STAT_ROBUST_STDEV - PS_STAT_ROBUST_QUARTILE - I have included all that computation in a single function, as opposed to - breaking it across several functions for one primary reason: they all - require the same basic initial processing steps (calculate the histogram, - etc.) however there is no currently defined means, in the SDRS, to - specify this computation as a separate step. If the above robust stat - measures were calcualted in separate functiosn, then much of the initial - processing would be duplicated. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorRobustStats(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - psHistogram *robustHistogram = NULL; - float binSize = 0.0; // Size of the histogram bins - float sigmaE = 0.0; - int LQBinNum = -1; // Bin num for lower quartile - int UQBinNum = -1; // Bin num for upper quartile - int i = 0; // Loop index variable - int maxBinNum = 0; - int maxBinCount = 0; - float dL = 0.0; - int numBins = 0; - psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN); - - // NOTE: The SDRS states that the sample quartiles must be used to - // determine the initial bin sizes. However, the sample quartiles are - // calculated based on a full sort of the data set, regardless of the - // size of the data set. We should consult with IfA to ensure that this - // is really required. - /* - if (isnan(stats->sampleUQ) || - isnan(stats->sampleLQ)) { - stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE; - p_psVectorSampleQuartiles(myVector, - maskVector, - maskVal, - stats); - } - */ - // Compute the initial bin size of the robust histogram. - p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats); - binSize = stats->clippedStdev / 10.0f; - - // If stats->clippedStdev == 0.0, then all data elements have the same - // value. Therefore, we can set the appropiate results and return. - if (fabs(binSize) <= FLT_EPSILON) { - if (stats->options & PS_STAT_ROBUST_MEAN) { - stats->robustMean = stats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - stats->robustMedian = stats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = stats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_STDEV) { - stats->robustStdev = 0.0; - } - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - stats->robustUQ = stats->clippedMean; - stats->robustLQ = stats->clippedMean; - } - psFree(tmpStats); - psFree(robustHistogram); - return; - } - - // Detemine minimum and maximum values in the data vector. - if (isnan(stats->min)) { - p_psVectorMin(myVector, maskVector, maskVal, stats); - } - if (isnan(stats->max)) { - p_psVectorMax(myVector, maskVector, maskVal, stats); - } - - // Create the histogram structure (yes, 2 is necessary, not 1). Also, - // if we get here, we know that binSize != 0.0. - numBins = 2 + (int) ((stats->max - stats->min) / binSize); - - robustHistogram = psHistogramAlloc(stats->min, - stats->max, - numBins); - // Populate the histogram array. - psHistogramVector(robustHistogram, myVector, maskVector, maskVal); - - // Smooth the histogram. - p_psVectorsmoothHistGaussian(robustHistogram, sigmaE/4.0f); - - LQBinNum = -1; - UQBinNum = -1; - for (i=0;inums->n;i++) { - if ((robustHistogram->nums->data.S32[i] <= stats->sampleLQ) && - (stats->sampleLQ <= robustHistogram->nums->data.S32[i])) { - LQBinNum = i; - } - - if ((robustHistogram->nums->data.S32[i] <= stats->sampleUQ) && - (stats->sampleUQ <= robustHistogram->nums->data.S32[i])) { - UQBinNum = i; - } - } - - // Determine the bin with the peak value in the range LQ to UQ. - maxBinNum = LQBinNum; - maxBinCount = robustHistogram->nums->data.S32[maxBinNum]; - for (i=LQBinNum;i<=UQBinNum;i++) { - if (robustHistogram->nums->data.S32[i] > maxBinCount) { - maxBinNum = i; - maxBinCount = robustHistogram->nums->data.S32[i]; - } - } - - dL = (stats->robustUQ - stats->robustLQ) / 8.0; - - // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL - // NOTE: This step is dependent on the functions in psMinimize.c being - // implemented. Currently, they are not. - - if (stats->options & PS_STAT_ROBUST_MEAN) { - stats->robustMean = 0.0; - } - - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - stats->robustMedian = 0.0; - } - if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = maxBinNum; - } - if (stats->options & PS_STAT_ROBUST_STDEV) { - stats->robustStdev = 0.0; - } - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - stats->robustUQ = 0.0; - stats->robustLQ = 0.0; - } - stats->robustNfit = 0.0; - stats->robustN50 = 0.0; - psFree(tmpStats); - psFree(robustHistogram); -} - /****************************************************************************** @@ -1229,22 +1787,10 @@ } - // Ensure that the data is of type PS_TYPE_F32. Eventually, more data - // types will be implemented. - if (in->type.type != PS_TYPE_F32) { - psAbort(__func__, - "Only data type PS_TYPE_F32 is currently supported (0x%x).", - in->type.type); - } - - // Ensure that the mask vector is of the proper size and type. + PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32); if (mask != NULL) { - if (in->n != mask->n) { - psAbort(__func__, - "Vector data and vector mask are of different sizes."); - } - if (mask->type.type != PS_TYPE_U8) { - psAbort(__func__, - "Vector mask must be type PS_TYPE_U8"); - } + PS_CHECK_NULL_VECTOR(mask); + PS_CHECK_EMPTY_VECTOR(mask); + PS_CHECK_VECTOR_SIZE_EQUAL(mask, in); + PS_CHECK_VECTOR_TYPE(mask, PS_TYPE_U8); } @@ -1281,4 +1827,6 @@ (stats->options & PS_STAT_ROBUST_QUARTILE)) { p_psVectorRobustStats(in, mask, maskVal, stats); + + printf("HMMMM stats->robustMode is %f\n", stats->robustMode); }