Index: /trunk/psLib/src/dataManip/psFunctions.c =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.c (revision 1019) +++ /trunk/psLib/src/dataManip/psFunctions.c (revision 1020) @@ -4,8 +4,9 @@ * routines. * - * @author George Gusciora, MHPCC + * This file will hold the functions for allocated, freeing, and evaluating + * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-06-07 00:32:53 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-14 19:32:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -103,4 +104,5 @@ psPolynomial1D *psPolynomial1DAlloc(int n) { + int i = 0; psPolynomial1D *newPoly = NULL; @@ -110,4 +112,9 @@ newPoly->coeffErr = (float *) psAlloc(n * sizeof(float)); newPoly->mask = (char *) psAlloc(n * sizeof(char)); + for (i=0;icoeff[i] = 0.0; + newPoly->coeffErr[i] = 0.0; + newPoly->mask[i] = 0.0; + } return(newPoly); @@ -117,4 +124,5 @@ { int x = 0; + int y = 0; psPolynomial2D *newPoly = NULL; @@ -131,5 +139,11 @@ newPoly->mask[x] = (char *) psAlloc(nY * sizeof(char)); } - + for (x=0;xcoeff[x][y] = 0.0; + newPoly->coeffErr[x][y] = 0.0; + newPoly->mask[x][y] = 0; + } + } return(newPoly); @@ -140,4 +154,5 @@ int x = 0; int y = 0; + int z = 0; psPolynomial3D *newPoly = NULL; @@ -160,4 +175,13 @@ } } + for (x=0;xcoeff[x][y][z] = 0.0; + newPoly->coeffErr[x][y][z] = 0.0; + newPoly->mask[x][y][z] = 0; + } + } + } return(newPoly); @@ -169,4 +193,5 @@ int x = 0; int y = 0; + int z = 0; psPolynomial4D *newPoly = NULL; @@ -195,4 +220,15 @@ } } + for (w=0;wcoeff[w][x][y][z] = 0.0; + newPoly->coeffErr[w][x][y][z] = 0.0; + newPoly->mask[w][x][y][z] = 0; + } + } + } + } return(newPoly); @@ -394,4 +430,5 @@ psDPolynomial1D *psDPolynomial1DAlloc(int n) { + int i = 0; psDPolynomial1D *newPoly = NULL; @@ -401,4 +438,9 @@ newPoly->coeffErr = (double *) psAlloc(n * sizeof(double)); newPoly->mask = (char *) psAlloc(n * sizeof(char)); + for (i=0;icoeff[i] = 0.0; + newPoly->coeffErr[i] = 0.0; + newPoly->mask[i] = 0.0; + } return(newPoly); @@ -408,4 +450,5 @@ { int x = 0; + int y = 0; psDPolynomial2D *newPoly = NULL; @@ -422,4 +465,11 @@ newPoly->mask[x] = (char *) psAlloc(nY * sizeof(char)); } + for (x=0;xcoeff[x][y] = 0.0; + newPoly->coeffErr[x][y] = 0.0; + newPoly->mask[x][y] = 0; + } + } return(newPoly); @@ -430,4 +480,5 @@ int x = 0; int y = 0; + int z = 0; psDPolynomial3D *newPoly = NULL; @@ -450,4 +501,13 @@ } } + for (x=0;xcoeff[x][y][z] = 0.0; + newPoly->coeffErr[x][y][z] = 0.0; + newPoly->mask[x][y][z] = 0; + } + } + } return(newPoly); @@ -459,4 +519,5 @@ int x = 0; int y = 0; + int z = 0; psDPolynomial4D *newPoly = NULL; @@ -485,4 +546,15 @@ } } + for (w=0;wcoeff[w][x][y][z] = 0.0; + newPoly->coeffErr[w][x][y][z] = 0.0; + newPoly->mask[w][x][y][z] = 0; + } + } + } + } return(newPoly); Index: /trunk/psLib/src/dataManip/psFunctions.h =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.h (revision 1019) +++ /trunk/psLib/src/dataManip/psFunctions.h (revision 1020) @@ -1,9 +1,24 @@ +/** @file psFunctions.h + * \brief Standard Mathematical Functions. + * \ingroup Stats + * + * This file will hold the prototypes for procedures which allocate, free, + * and evaluate various polynomials. Those polynomial structures are also + * defined here. + * + * @ingroup Stats + * + * @author Someone at IfA + * @author George Gusciora, MHPCC + * + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-14 19:32:42 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + #if !defined(PS_FUNCTIONS_H) #define PS_FUNCTIONS_H -/** \file psFunctions.h - * \brief Standard Mathematical Functions. - * \ingroup Stats - */ /** \addtogroup Stats @@ -23,17 +38,21 @@ -psVector *psGaussianDev(float mean, - float sigma, - int Npts); - -/*****************************************************************************/ +/** Produce a vector of random numbers from a Gaussian distribution with + the specified mean and sigma */ +psVector *psGaussianDev(float mean, ///< The mean of the Gaussian + float sigma, ///< The sigma of the Gaussian + int Npts); ///< The size of the vector + + + + /** One-dimensional polynomial */ typedef struct { - int n; ///< Number of terms - float *coeff; ///< Coefficients - float *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask + int n; ///< Number of terms + float *coeff; ///< Coefficients + float *coeffErr; ///< Error in coefficients + char *mask; ///< Coefficient mask } psPolynomial1D; Index: /trunk/psLib/src/dataManip/psStats.c =================================================================== --- /trunk/psLib/src/dataManip/psStats.c (revision 1019) +++ /trunk/psLib/src/dataManip/psStats.c (revision 1020) @@ -15,9 +15,7 @@ #include "float.h" -#define ROBUST_SIZE_THRESHOLD 30000 // Vectors that are large than this +#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 PI 3.141592653 +#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. #define CLIPPED_NUM_ITER_LB 1 #define CLIPPED_NUM_ITER_UB 10 @@ -26,5 +24,5 @@ #define true 1 #define false 0 -#define MYMAXFLOAT 1e99 +#define MYMAXFLOAT HUGE #ifndef DOXYGEN @@ -48,5 +46,5 @@ newStruct->sampleUQ = NAN; newStruct->sampleLQ = NAN; - newStruct->sampleLimit = ROBUST_SIZE_THRESHOLD; + newStruct->sampleLimit = DEFAULT_ROBUST_SIZE_THRESHOLD; newStruct->robustMean = NAN; newStruct->robustMedian = NAN; @@ -78,5 +76,14 @@ /****************************************************************************** - HISTOGRAM FUNCTIONS +psHistogramAlloc(lower, upper, n): allocate a uniform histogram structure +with the specifed upper and lower limits, and the specifed number of bins. +This routine will also set the bounds for each of the bins. + +Input: + lower + upper + n +Returns: + The histogram structure *****************************************************************************/ psHistogram *psHistogramAlloc(float lower, @@ -84,26 +91,31 @@ int n) { - int i = 0; - psHistogram *newHist = NULL; - float binSize = 0.0; - + int i = 0; // Loop index variable + psHistogram *newHist = NULL; // The new histogram structure + float binSize = 0.0; // The histogram bin size + + // GUS: Verify that this is the correct action. if (n == 0) { - psAbort(__func__, "psHistogram requested with 0 bins"); + return(NULL); } // GUS: Verify that this is the correct action. if (lower > upper) { - // psAbort(__func__, "psHistogram requested with lower>upper"); return(NULL); } + // Allocate memory for the new histogram structure. If there are N + // bins, then there are N+1 bounds to those bins. newHist = (psHistogram *) psAlloc(sizeof(psHistogram)); newHist->bounds = psVectorAlloc(n+1, PS_TYPE_F32); newHist->bounds->n = newHist->bounds->nalloc; + // Calculate the bounds for each bin. binSize = (upper - lower) / (float) n; for (i=0;ibounds->data.F32[i] = lower + (binSize * (float) i); } + + // Allocate the bins, and initialize them to zero. newHist->nums = psVectorAlloc(n, PS_TYPE_S32); newHist->nums->n = newHist->nums->nalloc; @@ -111,4 +123,6 @@ newHist->nums->data.S32[i] = 0; } + + // Initialize the other members. newHist->minNum = 0; newHist->maxNum = 0; @@ -118,10 +132,17 @@ } -// When this is called, the number of data elements in bounds -// is n. Therefore, the number of bins is n-1. +/****************************************************************************** +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. +Input: + bounds +Returns: + The histogram structure + *****************************************************************************/ psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds) { - psHistogram *newHist = NULL; - int i; + psHistogram *newHist = NULL; // The new histogram structure + int i; // Loop index variable // GUS: Verify that this is the correct action. @@ -131,4 +152,11 @@ } + // GUS: Verify that this is the correct action. + if (bounds->n <= 1) { + // psAbort(__func__, "psHistogram requested with NULL bounds"); + return(NULL); + } + + // Allocate memory for the new histogram structure. newHist = (psHistogram *) psAlloc(sizeof(psHistogram)); newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32); @@ -137,4 +165,7 @@ newHist->bounds->data.F32[i] = bounds->data.F32[i]; } + + // Allocate the bins, and initialize them to zero. If there are N bounds, + // then there are N-1 bins. newHist->nums = psVectorAlloc((bounds->n)-1, PS_TYPE_S32); newHist->nums->n = newHist->nums->nalloc; @@ -143,4 +174,5 @@ } + // Initialize the other members. newHist->minNum = 0; newHist->maxNum = 0; @@ -159,30 +191,30 @@ /***************************************************************************** - NOTE: Can we assume fixed size bins in the myHist data structure? +psHistogramVector(out, in, mask, maskVal): this procedure takes as input a +preallocated and initialized histogram structure. It fills the bins in that +histogram structure in accordance with the input data "in" and the, possibly +NULL, mask vector. - NOTE: This procedure assumes that the following has been already set - myHist->lower - myHist->upper - myHist->minVal - myHist->maxVal - This procedure sets - myHist->nums - myHist->minNum - myHist->maxNum +Inputs: + out + in + mask + maskVal +Returns: + The histogram structure "out". *****************************************************************************/ psHistogram *psHistogramVector(psHistogram *out, - psVector *in, - psVector *mask, + const psVector *restrict in, + const psVector *restrict mask, int maskVal) { - int i = 0; - int j = 0; - float binSize = 0.0; - int binNum = 0; - int numBins = 0; + int i = 0; // Loop index variable + int j = 0; // Loop index variable + float binSize = 0.0; // Histogram bin size + int binNum = 0; // A temporary bin number + int numBins = 0; // The total number of bins // GUS: Verify that this is the correct action. if (out == NULL) { - // psAbort(__func__, "psHistogram requested with NULL struct"); return(NULL); } @@ -190,56 +222,38 @@ // GUS: Verify that this is the correct action. if (in == NULL) { - // psAbort(__func__, "psHistogram requested with NULL struct"); return(out); } + if (in->type.type != PS_TYPE_F32) { + psAbort(__func__, + "Only data type PS_TYPE_F32 is currently supported (0x%x).", + in->type.type); + } + + if (mask != NULL) { + if (in->n != mask->n) { + psAbort(__func__, + "Vector data and vector mask are of different sizes."); + } + if (mask->type.type != PS_TYPE_U8) { + psAbort(__func__, "Vector mask must be type PS_TYPE_U8"); + } + } + // GUS: determine the correct action for a variety of other cases: + // in vector has 0 elements, and histogram structure has zero bins. + numBins = out->nums->n; - if (mask != NULL) { - for (i=0;in;i++) { - // Check if this pixel is masked, and if so, skip it. - if (!(mask->data.U8[i] & maskVal)) { - // Check if this pixel is below the minimum value, and if so - // count it, then skip it. - if (in->data.F32[i] < out->bounds->data.F32[0]) { - out->minNum++; - - // Check if this pixel is above the maximum value, and if so - // count it, then skip it. - } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { - out->maxNum++; - } else { - // If this is a uniform histogram, determining the correct - // number is trivial. - if (out->uniform == true) { - binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; - - binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) / - binSize); - (out->nums->data.S32[binNum])++; - // If this is a non-uniform histogram, determining the correct - // bin number requires a bit more work. - } else { - // GUS: This is slow. Put a smarter algorithm here to - // find the correct bin number (bin search, probably) - for (j=0;j<(out->bounds->n)-1;j++) { - if ((out->bounds->data.S32[j] <= in->data.F32[i]) && - (in->data.F32[i] <= out->bounds->data.S32[j+1])) { - (out->nums->data.S32[j])++; - } - } - } - } - } - } - } else { - for (i=0;in;i++) { + for (i=0;in;i++) { + // Check if this pixel is masked, and if so, skip it. + if ((mask == NULL) || + ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) { + // Check if this pixel is below the minimum value, and if so + // count it, then skip it. if (in->data.F32[i] < out->bounds->data.F32[0]) { - // Check if this pixel is below the minimum value, and if so - // count it, then skip it. out->minNum++; - } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { // Check if this pixel is above the maximum value, and if so // count it, then skip it. + } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { out->maxNum++; } else { @@ -251,6 +265,4 @@ binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) / binSize); - (out->nums->data.S32[binNum])++; - // If this is a non-uniform histogram, determining the correct // bin number requires a bit more work. @@ -267,15 +279,20 @@ } } - } return(out); } -void p_printVector(psVector *myVector, - psVector *maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; +/***************************************************************************** +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++) { @@ -291,5 +308,14 @@ ****************************************************************************** ****************************************************************************** - MISC STATISTICAL FUNCTIONS + 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? ****************************************************************************** ****************************************************************************** @@ -298,6 +324,11 @@ /****************************************************************************** - ASSUMPTION: the mean is always calculated exactly. Robust means are never - calculated in this routine. +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, @@ -306,10 +337,12 @@ psStats *stats) { - int i = 0; - float mean = 0.0; - int count = 0; - float rangeMin = 0.0; - float rangeMax = 0.0; - + 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; @@ -317,4 +350,5 @@ 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]) && @@ -355,4 +389,15 @@ } +/****************************************************************************** +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, @@ -360,9 +405,10 @@ psStats *stats) { - int i = 0; - float max = -MYMAXFLOAT; - float rangeMin = 0.0; - float rangeMax = 0.0; - + 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; @@ -408,4 +454,15 @@ } +/****************************************************************************** +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, @@ -413,8 +470,8 @@ psStats *stats) { - int i = 0; - float min = MYMAXFLOAT; - float rangeMin = 0.0; - float rangeMax = 0.0; + 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) { @@ -462,6 +519,14 @@ /****************************************************************************** - This routine calculates the number of non-masked pixels in the vector. - that fall within the min/max range, if given. +p_psVectorSampleMax(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, @@ -470,8 +535,8 @@ psStats *stats) { - int i = 0; - int numData = 0; - float rangeMin = 0.0; - float rangeMax = 0.0; + 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) { @@ -512,4 +577,15 @@ +/****************************************************************************** +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, @@ -517,12 +593,12 @@ psStats *stats) { - psVector *unsortedVector = NULL; - psVector *sortedVector = NULL; - int count = 0; - int i = 0; - int nValues = 0; - float rangeMin = 0.0; - float rangeMax = 0.0; - psStats *stats2 = NULL; + 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 @@ -531,4 +607,6 @@ 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); @@ -548,5 +626,5 @@ // Determine how many data points fit inside this min/max range - // and are not maxed, IF the maskVector is not NULL> + // and are not masked, IF the maskVector is not NULL> nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); @@ -620,10 +698,8 @@ float sigma) { - int i = 0; - int j = 0; - float tmpf = 0.0; - float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; - // The coefficients used in the histogram - // smoothing calculation. + 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++) { @@ -631,20 +707,19 @@ // If sigma does not equal zero, then we use Gaussian smoothing. #ifdef DARWIN - tmpf = (float) sqrt(2.0f * PI * sigma * sigma); + tmpf = (float) sqrt(2.0f * M_PI * sigma * sigma); #else - tmpf = sqrtf(2.0f * PI * sigma * sigma); + 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; + (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 instead - * use boxcar smoothing. + * the above code will divide by zero. Therefore, we don't need + * to smooth the data. */ - gaussianCoefs[i] = 1.0f / (1.0f + (2.0f * (float) GAUSS_WIDTH)); + return; } } @@ -662,8 +737,13 @@ /****************************************************************************** - p_psVectorSampleQuartiles() - This procedure calculates the upper and/or lower quartiles of the - data set. - +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, @@ -672,35 +752,33 @@ psStats *stats) { - psVector *unsortedVector = NULL; - psVector *sortedVector = NULL; - int count = 0; - int ind = 0; - int i = 0; - int nValues = 0; - float rangeMin = 0.0; - float rangeMax = 0.0; + 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. - /* GUS: When IfA provides the algorithm, insert is here. - if (myVector->n > stats->sampleLimit) { - 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. - psStatsFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } - */ + if (myVector->n > stats->sampleLimit) { + psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); + 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. + psStatsFree(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 @@ -757,8 +835,6 @@ // Calculate the quartile points exactly. - ind = 3 * (nValues / 4); - stats->sampleUQ = sortedVector->data.F32[ind]; - ind = (nValues / 4); - stats->sampleLQ = sortedVector->data.F32[ind]; + stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)]; + stats->sampleLQ = sortedVector->data.F32[nValues / 4]; // Free the temporary data structures. @@ -770,6 +846,6 @@ /****************************************************************************** - p_psVectorRobustStats(): this procedure calculates a variety of robust - stat measures: +p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure +calculates a variety of robust stat measures: PS_STAT_ROBUST_MEAN PS_STAT_ROBUST_MEDIAN @@ -784,4 +860,11 @@ 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, @@ -791,12 +874,13 @@ { psHistogram *robustHistogram = NULL; - float binSize = 0.0; + float binSize = 0.0; // Size of the histogram bins float sigmaE = 0.0; - int LQBinNum = -1; - int UQBinNum = -1; - int i = 0; // Loop index variable. + 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; // NOTE: The SDRS states that the sample quartiles must be used to @@ -827,14 +911,12 @@ } - // Create the histogram structure. + // Create the histogram structure (yes, 2 is necessary, not 1). + numBins = 2 + (int) ((stats->max - stats->min) / binSize); + robustHistogram = psHistogramAlloc(stats->min, stats->max, - 1 + (int) ((stats->max - stats->min) / binSize)); - - + numBins); // Populate the histogram array. - // GUS: fix this - // robustHistogram = psHistogramVector(robustHistogram, myVector); - // + psHistogramVector(robustHistogram, myVector, maskVector, maskVal); // Smooth the histogram. @@ -868,8 +950,6 @@ // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL - // What algorithm should I use to do this? - - // The resulting fit parameters are the robust mean, mean_r, and sigma - // What is the mean_r? + // GUS: This step is dependent on the functions in psMinimize.c being + // implemented. Currently, they are not. if (stats->options & PS_STAT_ROBUST_MEAN) { @@ -895,10 +975,17 @@ } -/***************************************************************************** - NOTE: This function assumes that p_psVectorMean() has already been called - and the correct value is stored in stats->sampleMean. +/****************************************************************************** +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. +NOTE: the mean is always calculated exactly. Robust means are never +calculated in this routine. *****************************************************************************/ void p_psVectorSampleStdev(const psVector *restrict myVector, @@ -907,14 +994,16 @@ psStats *stats) { - int i = 0; - int countInt = 0; - float countFloat = 0.0; - float mean = 0.0; - float diff = 0.0; - float sumSquares = 0.0; - float sumDiffs = 0.0; - float rangeMin = 0.0; - float rangeMax = 0.0; - + 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); @@ -979,4 +1068,16 @@ } +/****************************************************************************** +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, @@ -984,10 +1085,13 @@ psStats *stats) { - int i = 0; - int j = 0; - float clippedMean = 0.0; - float clippedStdev = 0.0; - psVector *tmpMask = NULL; - + 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))) { @@ -996,4 +1100,5 @@ } + // Endure that stats->clipSigma is within the proper range. if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && (stats->clipSigma <= CLIPPED_SIGMA_UB))) { @@ -1002,8 +1107,12 @@ } - + // 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++) { @@ -1026,9 +1135,11 @@ 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++) { - // printf("p_psVectorClippedStats(): Iteration %d (%f %f)\n", i, - // clippedMean, clippedStdev); for (j=0;jn;j++) { // a) Exclude all values x_i for which |x_i - x| > K * stdev @@ -1038,6 +1149,4 @@ } // b) compute new mean and stdev - // GUS: I should probably create a new struct here since the - // following calls will overwrite any old values in sampleMean. p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); @@ -1050,4 +1159,6 @@ } } + stats->sampleMean = oldStanMean; + stats->sampleStdev= oldStanStdev; // 7. The last calcuated value of x is the cliped mean. @@ -1066,9 +1177,19 @@ /****************************************************************************** +psVectorStats(myVector, maskVector, maskVal, stats): this is the public API +function which calls the above private stats functions based on what bits +were set in stats->options. - NOTE: The current strategy is to implement everything assuming that all - input data is of type PS_TYPE_F32. Once the basic code is in place, - we will macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64. +Inputs + myVector + maskVector + maskVal + stats +Returns + The stats structure. +NOTE: The current strategy is to implement everything assuming that all +input data is of type PS_TYPE_F32. Once the basic code is in place, we will +macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64. *****************************************************************************/ psStats *psVectorStats(psStats *stats, @@ -1077,4 +1198,5 @@ unsigned int maskVal) { + // GUS: Verify that this is the correct action. if (in == NULL) { return(stats); @@ -1084,4 +1206,6 @@ } + // 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__, @@ -1090,4 +1214,5 @@ } + // Ensure that the mask vector is of the proper size and type. if (mask != NULL) { if (in->n != mask->n) { @@ -1124,4 +1249,7 @@ } + // Since the various robust stats quantities share much computation, they + // are grouped together in a single private function: + // p_psVectorRobustStats() if ((stats->options & PS_STAT_ROBUST_MEAN) || (stats->options & PS_STAT_ROBUST_MEDIAN) || @@ -1147,8 +1275,4 @@ } - // OLD CODE: Should we check for an unknown option? - // default: - // psAbort(__func__, "Unknown options 0x%x.\n", stats->options); - return(stats); } Index: /trunk/psLib/src/dataManip/psStats.h =================================================================== --- /trunk/psLib/src/dataManip/psStats.h (revision 1019) +++ /trunk/psLib/src/dataManip/psStats.h (revision 1020) @@ -1,11 +1,20 @@ +/** @file psStats.h + * \brief basic statistical operations + * @ingroup Stats + * + * This file will hold the definition of the histogram and stats data + * structures. It also contains prototypes for procedures which operate + * on those data structures. + * + * @author George Gusciora, MHPCC + * + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-14 19:33:09 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ #if !defined(PS_STATS_H) #define PS_STATS_H -/****************************************************************************** - This file will histogram/statistical functions and data structures. - *****************************************************************************/ -/** \file psStats.h - * \brief basic statistical operations - * \ingroup Stats - */ + #include "psVector.h" @@ -36,5 +45,7 @@ -/** This is the generic statistics structure */ +/** This is the generic statistics structure. It contails the data members + for the various statistic values. It contains the options member to + specifiy which statistics must be calculated. */ typedef struct { @@ -74,8 +85,8 @@ ); -/** Constructor */ +/** A constructor for the stats structure.*/ psStats *psStatsAlloc(psStatsOptions options); ///< Statistics to measure -/** Destructor */ +/** A destructor for the stats structure.*/ void psStatsFree(psStats *restrict stats); ///< Stats structure to destroy @@ -83,4 +94,5 @@ Histogram functions and data structures. *****************************************************************************/ +/** The basic histogram structure which contains bounds and bins. */ typedef struct { @@ -111,6 +123,6 @@ /** Calculate a histogram \ingroup MathGroup **/ psHistogram *psHistogramVector (psHistogram *out, ///< Histogram data - psVector *in, ///< Vector to analyse - psVector *mask, ///< Mask dat for input vector + const psVector *restrict in, ///< Vector to analyse + const psVector *restrict mask, ///< Mask dat for input vector int maskVal); ///< Mask value Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 1019) +++ /trunk/psLib/src/math/psPolynomial.c (revision 1020) @@ -4,8 +4,9 @@ * routines. * - * @author George Gusciora, MHPCC + * This file will hold the functions for allocated, freeing, and evaluating + * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-06-07 00:32:53 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-14 19:32:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -103,4 +104,5 @@ psPolynomial1D *psPolynomial1DAlloc(int n) { + int i = 0; psPolynomial1D *newPoly = NULL; @@ -110,4 +112,9 @@ newPoly->coeffErr = (float *) psAlloc(n * sizeof(float)); newPoly->mask = (char *) psAlloc(n * sizeof(char)); + for (i=0;icoeff[i] = 0.0; + newPoly->coeffErr[i] = 0.0; + newPoly->mask[i] = 0.0; + } return(newPoly); @@ -117,4 +124,5 @@ { int x = 0; + int y = 0; psPolynomial2D *newPoly = NULL; @@ -131,5 +139,11 @@ newPoly->mask[x] = (char *) psAlloc(nY * sizeof(char)); } - + for (x=0;xcoeff[x][y] = 0.0; + newPoly->coeffErr[x][y] = 0.0; + newPoly->mask[x][y] = 0; + } + } return(newPoly); @@ -140,4 +154,5 @@ int x = 0; int y = 0; + int z = 0; psPolynomial3D *newPoly = NULL; @@ -160,4 +175,13 @@ } } + for (x=0;xcoeff[x][y][z] = 0.0; + newPoly->coeffErr[x][y][z] = 0.0; + newPoly->mask[x][y][z] = 0; + } + } + } return(newPoly); @@ -169,4 +193,5 @@ int x = 0; int y = 0; + int z = 0; psPolynomial4D *newPoly = NULL; @@ -195,4 +220,15 @@ } } + for (w=0;wcoeff[w][x][y][z] = 0.0; + newPoly->coeffErr[w][x][y][z] = 0.0; + newPoly->mask[w][x][y][z] = 0; + } + } + } + } return(newPoly); @@ -394,4 +430,5 @@ psDPolynomial1D *psDPolynomial1DAlloc(int n) { + int i = 0; psDPolynomial1D *newPoly = NULL; @@ -401,4 +438,9 @@ newPoly->coeffErr = (double *) psAlloc(n * sizeof(double)); newPoly->mask = (char *) psAlloc(n * sizeof(char)); + for (i=0;icoeff[i] = 0.0; + newPoly->coeffErr[i] = 0.0; + newPoly->mask[i] = 0.0; + } return(newPoly); @@ -408,4 +450,5 @@ { int x = 0; + int y = 0; psDPolynomial2D *newPoly = NULL; @@ -422,4 +465,11 @@ newPoly->mask[x] = (char *) psAlloc(nY * sizeof(char)); } + for (x=0;xcoeff[x][y] = 0.0; + newPoly->coeffErr[x][y] = 0.0; + newPoly->mask[x][y] = 0; + } + } return(newPoly); @@ -430,4 +480,5 @@ int x = 0; int y = 0; + int z = 0; psDPolynomial3D *newPoly = NULL; @@ -450,4 +501,13 @@ } } + for (x=0;xcoeff[x][y][z] = 0.0; + newPoly->coeffErr[x][y][z] = 0.0; + newPoly->mask[x][y][z] = 0; + } + } + } return(newPoly); @@ -459,4 +519,5 @@ int x = 0; int y = 0; + int z = 0; psDPolynomial4D *newPoly = NULL; @@ -485,4 +546,15 @@ } } + for (w=0;wcoeff[w][x][y][z] = 0.0; + newPoly->coeffErr[w][x][y][z] = 0.0; + newPoly->mask[w][x][y][z] = 0; + } + } + } + } return(newPoly); Index: /trunk/psLib/src/math/psPolynomial.h =================================================================== --- /trunk/psLib/src/math/psPolynomial.h (revision 1019) +++ /trunk/psLib/src/math/psPolynomial.h (revision 1020) @@ -1,9 +1,24 @@ +/** @file psFunctions.h + * \brief Standard Mathematical Functions. + * \ingroup Stats + * + * This file will hold the prototypes for procedures which allocate, free, + * and evaluate various polynomials. Those polynomial structures are also + * defined here. + * + * @ingroup Stats + * + * @author Someone at IfA + * @author George Gusciora, MHPCC + * + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-14 19:32:42 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + #if !defined(PS_FUNCTIONS_H) #define PS_FUNCTIONS_H -/** \file psFunctions.h - * \brief Standard Mathematical Functions. - * \ingroup Stats - */ /** \addtogroup Stats @@ -23,17 +38,21 @@ -psVector *psGaussianDev(float mean, - float sigma, - int Npts); - -/*****************************************************************************/ +/** Produce a vector of random numbers from a Gaussian distribution with + the specified mean and sigma */ +psVector *psGaussianDev(float mean, ///< The mean of the Gaussian + float sigma, ///< The sigma of the Gaussian + int Npts); ///< The size of the vector + + + + /** One-dimensional polynomial */ typedef struct { - int n; ///< Number of terms - float *coeff; ///< Coefficients - float *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask + int n; ///< Number of terms + float *coeff; ///< Coefficients + float *coeffErr; ///< Error in coefficients + char *mask; ///< Coefficient mask } psPolynomial1D; Index: /trunk/psLib/src/math/psSpline.c =================================================================== --- /trunk/psLib/src/math/psSpline.c (revision 1019) +++ /trunk/psLib/src/math/psSpline.c (revision 1020) @@ -4,8 +4,9 @@ * routines. * - * @author George Gusciora, MHPCC + * This file will hold the functions for allocated, freeing, and evaluating + * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-06-07 00:32:53 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-14 19:32:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -103,4 +104,5 @@ psPolynomial1D *psPolynomial1DAlloc(int n) { + int i = 0; psPolynomial1D *newPoly = NULL; @@ -110,4 +112,9 @@ newPoly->coeffErr = (float *) psAlloc(n * sizeof(float)); newPoly->mask = (char *) psAlloc(n * sizeof(char)); + for (i=0;icoeff[i] = 0.0; + newPoly->coeffErr[i] = 0.0; + newPoly->mask[i] = 0.0; + } return(newPoly); @@ -117,4 +124,5 @@ { int x = 0; + int y = 0; psPolynomial2D *newPoly = NULL; @@ -131,5 +139,11 @@ newPoly->mask[x] = (char *) psAlloc(nY * sizeof(char)); } - + for (x=0;xcoeff[x][y] = 0.0; + newPoly->coeffErr[x][y] = 0.0; + newPoly->mask[x][y] = 0; + } + } return(newPoly); @@ -140,4 +154,5 @@ int x = 0; int y = 0; + int z = 0; psPolynomial3D *newPoly = NULL; @@ -160,4 +175,13 @@ } } + for (x=0;xcoeff[x][y][z] = 0.0; + newPoly->coeffErr[x][y][z] = 0.0; + newPoly->mask[x][y][z] = 0; + } + } + } return(newPoly); @@ -169,4 +193,5 @@ int x = 0; int y = 0; + int z = 0; psPolynomial4D *newPoly = NULL; @@ -195,4 +220,15 @@ } } + for (w=0;wcoeff[w][x][y][z] = 0.0; + newPoly->coeffErr[w][x][y][z] = 0.0; + newPoly->mask[w][x][y][z] = 0; + } + } + } + } return(newPoly); @@ -394,4 +430,5 @@ psDPolynomial1D *psDPolynomial1DAlloc(int n) { + int i = 0; psDPolynomial1D *newPoly = NULL; @@ -401,4 +438,9 @@ newPoly->coeffErr = (double *) psAlloc(n * sizeof(double)); newPoly->mask = (char *) psAlloc(n * sizeof(char)); + for (i=0;icoeff[i] = 0.0; + newPoly->coeffErr[i] = 0.0; + newPoly->mask[i] = 0.0; + } return(newPoly); @@ -408,4 +450,5 @@ { int x = 0; + int y = 0; psDPolynomial2D *newPoly = NULL; @@ -422,4 +465,11 @@ newPoly->mask[x] = (char *) psAlloc(nY * sizeof(char)); } + for (x=0;xcoeff[x][y] = 0.0; + newPoly->coeffErr[x][y] = 0.0; + newPoly->mask[x][y] = 0; + } + } return(newPoly); @@ -430,4 +480,5 @@ int x = 0; int y = 0; + int z = 0; psDPolynomial3D *newPoly = NULL; @@ -450,4 +501,13 @@ } } + for (x=0;xcoeff[x][y][z] = 0.0; + newPoly->coeffErr[x][y][z] = 0.0; + newPoly->mask[x][y][z] = 0; + } + } + } return(newPoly); @@ -459,4 +519,5 @@ int x = 0; int y = 0; + int z = 0; psDPolynomial4D *newPoly = NULL; @@ -485,4 +546,15 @@ } } + for (w=0;wcoeff[w][x][y][z] = 0.0; + newPoly->coeffErr[w][x][y][z] = 0.0; + newPoly->mask[w][x][y][z] = 0; + } + } + } + } return(newPoly); Index: /trunk/psLib/src/math/psSpline.h =================================================================== --- /trunk/psLib/src/math/psSpline.h (revision 1019) +++ /trunk/psLib/src/math/psSpline.h (revision 1020) @@ -1,9 +1,24 @@ +/** @file psFunctions.h + * \brief Standard Mathematical Functions. + * \ingroup Stats + * + * This file will hold the prototypes for procedures which allocate, free, + * and evaluate various polynomials. Those polynomial structures are also + * defined here. + * + * @ingroup Stats + * + * @author Someone at IfA + * @author George Gusciora, MHPCC + * + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-14 19:32:42 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + #if !defined(PS_FUNCTIONS_H) #define PS_FUNCTIONS_H -/** \file psFunctions.h - * \brief Standard Mathematical Functions. - * \ingroup Stats - */ /** \addtogroup Stats @@ -23,17 +38,21 @@ -psVector *psGaussianDev(float mean, - float sigma, - int Npts); - -/*****************************************************************************/ +/** Produce a vector of random numbers from a Gaussian distribution with + the specified mean and sigma */ +psVector *psGaussianDev(float mean, ///< The mean of the Gaussian + float sigma, ///< The sigma of the Gaussian + int Npts); ///< The size of the vector + + + + /** One-dimensional polynomial */ typedef struct { - int n; ///< Number of terms - float *coeff; ///< Coefficients - float *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask + int n; ///< Number of terms + float *coeff; ///< Coefficients + float *coeffErr; ///< Error in coefficients + char *mask; ///< Coefficient mask } psPolynomial1D; Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 1019) +++ /trunk/psLib/src/math/psStats.c (revision 1020) @@ -15,9 +15,7 @@ #include "float.h" -#define ROBUST_SIZE_THRESHOLD 30000 // Vectors that are large than this +#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 PI 3.141592653 +#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. #define CLIPPED_NUM_ITER_LB 1 #define CLIPPED_NUM_ITER_UB 10 @@ -26,5 +24,5 @@ #define true 1 #define false 0 -#define MYMAXFLOAT 1e99 +#define MYMAXFLOAT HUGE #ifndef DOXYGEN @@ -48,5 +46,5 @@ newStruct->sampleUQ = NAN; newStruct->sampleLQ = NAN; - newStruct->sampleLimit = ROBUST_SIZE_THRESHOLD; + newStruct->sampleLimit = DEFAULT_ROBUST_SIZE_THRESHOLD; newStruct->robustMean = NAN; newStruct->robustMedian = NAN; @@ -78,5 +76,14 @@ /****************************************************************************** - HISTOGRAM FUNCTIONS +psHistogramAlloc(lower, upper, n): allocate a uniform histogram structure +with the specifed upper and lower limits, and the specifed number of bins. +This routine will also set the bounds for each of the bins. + +Input: + lower + upper + n +Returns: + The histogram structure *****************************************************************************/ psHistogram *psHistogramAlloc(float lower, @@ -84,26 +91,31 @@ int n) { - int i = 0; - psHistogram *newHist = NULL; - float binSize = 0.0; - + int i = 0; // Loop index variable + psHistogram *newHist = NULL; // The new histogram structure + float binSize = 0.0; // The histogram bin size + + // GUS: Verify that this is the correct action. if (n == 0) { - psAbort(__func__, "psHistogram requested with 0 bins"); + return(NULL); } // GUS: Verify that this is the correct action. if (lower > upper) { - // psAbort(__func__, "psHistogram requested with lower>upper"); return(NULL); } + // Allocate memory for the new histogram structure. If there are N + // bins, then there are N+1 bounds to those bins. newHist = (psHistogram *) psAlloc(sizeof(psHistogram)); newHist->bounds = psVectorAlloc(n+1, PS_TYPE_F32); newHist->bounds->n = newHist->bounds->nalloc; + // Calculate the bounds for each bin. binSize = (upper - lower) / (float) n; for (i=0;ibounds->data.F32[i] = lower + (binSize * (float) i); } + + // Allocate the bins, and initialize them to zero. newHist->nums = psVectorAlloc(n, PS_TYPE_S32); newHist->nums->n = newHist->nums->nalloc; @@ -111,4 +123,6 @@ newHist->nums->data.S32[i] = 0; } + + // Initialize the other members. newHist->minNum = 0; newHist->maxNum = 0; @@ -118,10 +132,17 @@ } -// When this is called, the number of data elements in bounds -// is n. Therefore, the number of bins is n-1. +/****************************************************************************** +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. +Input: + bounds +Returns: + The histogram structure + *****************************************************************************/ psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds) { - psHistogram *newHist = NULL; - int i; + psHistogram *newHist = NULL; // The new histogram structure + int i; // Loop index variable // GUS: Verify that this is the correct action. @@ -131,4 +152,11 @@ } + // GUS: Verify that this is the correct action. + if (bounds->n <= 1) { + // psAbort(__func__, "psHistogram requested with NULL bounds"); + return(NULL); + } + + // Allocate memory for the new histogram structure. newHist = (psHistogram *) psAlloc(sizeof(psHistogram)); newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32); @@ -137,4 +165,7 @@ newHist->bounds->data.F32[i] = bounds->data.F32[i]; } + + // Allocate the bins, and initialize them to zero. If there are N bounds, + // then there are N-1 bins. newHist->nums = psVectorAlloc((bounds->n)-1, PS_TYPE_S32); newHist->nums->n = newHist->nums->nalloc; @@ -143,4 +174,5 @@ } + // Initialize the other members. newHist->minNum = 0; newHist->maxNum = 0; @@ -159,30 +191,30 @@ /***************************************************************************** - NOTE: Can we assume fixed size bins in the myHist data structure? +psHistogramVector(out, in, mask, maskVal): this procedure takes as input a +preallocated and initialized histogram structure. It fills the bins in that +histogram structure in accordance with the input data "in" and the, possibly +NULL, mask vector. - NOTE: This procedure assumes that the following has been already set - myHist->lower - myHist->upper - myHist->minVal - myHist->maxVal - This procedure sets - myHist->nums - myHist->minNum - myHist->maxNum +Inputs: + out + in + mask + maskVal +Returns: + The histogram structure "out". *****************************************************************************/ psHistogram *psHistogramVector(psHistogram *out, - psVector *in, - psVector *mask, + const psVector *restrict in, + const psVector *restrict mask, int maskVal) { - int i = 0; - int j = 0; - float binSize = 0.0; - int binNum = 0; - int numBins = 0; + int i = 0; // Loop index variable + int j = 0; // Loop index variable + float binSize = 0.0; // Histogram bin size + int binNum = 0; // A temporary bin number + int numBins = 0; // The total number of bins // GUS: Verify that this is the correct action. if (out == NULL) { - // psAbort(__func__, "psHistogram requested with NULL struct"); return(NULL); } @@ -190,56 +222,38 @@ // GUS: Verify that this is the correct action. if (in == NULL) { - // psAbort(__func__, "psHistogram requested with NULL struct"); return(out); } + if (in->type.type != PS_TYPE_F32) { + psAbort(__func__, + "Only data type PS_TYPE_F32 is currently supported (0x%x).", + in->type.type); + } + + if (mask != NULL) { + if (in->n != mask->n) { + psAbort(__func__, + "Vector data and vector mask are of different sizes."); + } + if (mask->type.type != PS_TYPE_U8) { + psAbort(__func__, "Vector mask must be type PS_TYPE_U8"); + } + } + // GUS: determine the correct action for a variety of other cases: + // in vector has 0 elements, and histogram structure has zero bins. + numBins = out->nums->n; - if (mask != NULL) { - for (i=0;in;i++) { - // Check if this pixel is masked, and if so, skip it. - if (!(mask->data.U8[i] & maskVal)) { - // Check if this pixel is below the minimum value, and if so - // count it, then skip it. - if (in->data.F32[i] < out->bounds->data.F32[0]) { - out->minNum++; - - // Check if this pixel is above the maximum value, and if so - // count it, then skip it. - } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { - out->maxNum++; - } else { - // If this is a uniform histogram, determining the correct - // number is trivial. - if (out->uniform == true) { - binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; - - binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) / - binSize); - (out->nums->data.S32[binNum])++; - // If this is a non-uniform histogram, determining the correct - // bin number requires a bit more work. - } else { - // GUS: This is slow. Put a smarter algorithm here to - // find the correct bin number (bin search, probably) - for (j=0;j<(out->bounds->n)-1;j++) { - if ((out->bounds->data.S32[j] <= in->data.F32[i]) && - (in->data.F32[i] <= out->bounds->data.S32[j+1])) { - (out->nums->data.S32[j])++; - } - } - } - } - } - } - } else { - for (i=0;in;i++) { + for (i=0;in;i++) { + // Check if this pixel is masked, and if so, skip it. + if ((mask == NULL) || + ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) { + // Check if this pixel is below the minimum value, and if so + // count it, then skip it. if (in->data.F32[i] < out->bounds->data.F32[0]) { - // Check if this pixel is below the minimum value, and if so - // count it, then skip it. out->minNum++; - } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { // Check if this pixel is above the maximum value, and if so // count it, then skip it. + } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { out->maxNum++; } else { @@ -251,6 +265,4 @@ binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) / binSize); - (out->nums->data.S32[binNum])++; - // If this is a non-uniform histogram, determining the correct // bin number requires a bit more work. @@ -267,15 +279,20 @@ } } - } return(out); } -void p_printVector(psVector *myVector, - psVector *maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; +/***************************************************************************** +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++) { @@ -291,5 +308,14 @@ ****************************************************************************** ****************************************************************************** - MISC STATISTICAL FUNCTIONS + 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? ****************************************************************************** ****************************************************************************** @@ -298,6 +324,11 @@ /****************************************************************************** - ASSUMPTION: the mean is always calculated exactly. Robust means are never - calculated in this routine. +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, @@ -306,10 +337,12 @@ psStats *stats) { - int i = 0; - float mean = 0.0; - int count = 0; - float rangeMin = 0.0; - float rangeMax = 0.0; - + 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; @@ -317,4 +350,5 @@ 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]) && @@ -355,4 +389,15 @@ } +/****************************************************************************** +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, @@ -360,9 +405,10 @@ psStats *stats) { - int i = 0; - float max = -MYMAXFLOAT; - float rangeMin = 0.0; - float rangeMax = 0.0; - + 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; @@ -408,4 +454,15 @@ } +/****************************************************************************** +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, @@ -413,8 +470,8 @@ psStats *stats) { - int i = 0; - float min = MYMAXFLOAT; - float rangeMin = 0.0; - float rangeMax = 0.0; + 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) { @@ -462,6 +519,14 @@ /****************************************************************************** - This routine calculates the number of non-masked pixels in the vector. - that fall within the min/max range, if given. +p_psVectorSampleMax(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, @@ -470,8 +535,8 @@ psStats *stats) { - int i = 0; - int numData = 0; - float rangeMin = 0.0; - float rangeMax = 0.0; + 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) { @@ -512,4 +577,15 @@ +/****************************************************************************** +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, @@ -517,12 +593,12 @@ psStats *stats) { - psVector *unsortedVector = NULL; - psVector *sortedVector = NULL; - int count = 0; - int i = 0; - int nValues = 0; - float rangeMin = 0.0; - float rangeMax = 0.0; - psStats *stats2 = NULL; + 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 @@ -531,4 +607,6 @@ 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); @@ -548,5 +626,5 @@ // Determine how many data points fit inside this min/max range - // and are not maxed, IF the maskVector is not NULL> + // and are not masked, IF the maskVector is not NULL> nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); @@ -620,10 +698,8 @@ float sigma) { - int i = 0; - int j = 0; - float tmpf = 0.0; - float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; - // The coefficients used in the histogram - // smoothing calculation. + 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++) { @@ -631,20 +707,19 @@ // If sigma does not equal zero, then we use Gaussian smoothing. #ifdef DARWIN - tmpf = (float) sqrt(2.0f * PI * sigma * sigma); + tmpf = (float) sqrt(2.0f * M_PI * sigma * sigma); #else - tmpf = sqrtf(2.0f * PI * sigma * sigma); + 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; + (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 instead - * use boxcar smoothing. + * the above code will divide by zero. Therefore, we don't need + * to smooth the data. */ - gaussianCoefs[i] = 1.0f / (1.0f + (2.0f * (float) GAUSS_WIDTH)); + return; } } @@ -662,8 +737,13 @@ /****************************************************************************** - p_psVectorSampleQuartiles() - This procedure calculates the upper and/or lower quartiles of the - data set. - +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, @@ -672,35 +752,33 @@ psStats *stats) { - psVector *unsortedVector = NULL; - psVector *sortedVector = NULL; - int count = 0; - int ind = 0; - int i = 0; - int nValues = 0; - float rangeMin = 0.0; - float rangeMax = 0.0; + 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. - /* GUS: When IfA provides the algorithm, insert is here. - if (myVector->n > stats->sampleLimit) { - 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. - psStatsFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } - */ + if (myVector->n > stats->sampleLimit) { + psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); + 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. + psStatsFree(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 @@ -757,8 +835,6 @@ // Calculate the quartile points exactly. - ind = 3 * (nValues / 4); - stats->sampleUQ = sortedVector->data.F32[ind]; - ind = (nValues / 4); - stats->sampleLQ = sortedVector->data.F32[ind]; + stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)]; + stats->sampleLQ = sortedVector->data.F32[nValues / 4]; // Free the temporary data structures. @@ -770,6 +846,6 @@ /****************************************************************************** - p_psVectorRobustStats(): this procedure calculates a variety of robust - stat measures: +p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure +calculates a variety of robust stat measures: PS_STAT_ROBUST_MEAN PS_STAT_ROBUST_MEDIAN @@ -784,4 +860,11 @@ 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, @@ -791,12 +874,13 @@ { psHistogram *robustHistogram = NULL; - float binSize = 0.0; + float binSize = 0.0; // Size of the histogram bins float sigmaE = 0.0; - int LQBinNum = -1; - int UQBinNum = -1; - int i = 0; // Loop index variable. + 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; // NOTE: The SDRS states that the sample quartiles must be used to @@ -827,14 +911,12 @@ } - // Create the histogram structure. + // Create the histogram structure (yes, 2 is necessary, not 1). + numBins = 2 + (int) ((stats->max - stats->min) / binSize); + robustHistogram = psHistogramAlloc(stats->min, stats->max, - 1 + (int) ((stats->max - stats->min) / binSize)); - - + numBins); // Populate the histogram array. - // GUS: fix this - // robustHistogram = psHistogramVector(robustHistogram, myVector); - // + psHistogramVector(robustHistogram, myVector, maskVector, maskVal); // Smooth the histogram. @@ -868,8 +950,6 @@ // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL - // What algorithm should I use to do this? - - // The resulting fit parameters are the robust mean, mean_r, and sigma - // What is the mean_r? + // GUS: This step is dependent on the functions in psMinimize.c being + // implemented. Currently, they are not. if (stats->options & PS_STAT_ROBUST_MEAN) { @@ -895,10 +975,17 @@ } -/***************************************************************************** - NOTE: This function assumes that p_psVectorMean() has already been called - and the correct value is stored in stats->sampleMean. +/****************************************************************************** +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. +NOTE: the mean is always calculated exactly. Robust means are never +calculated in this routine. *****************************************************************************/ void p_psVectorSampleStdev(const psVector *restrict myVector, @@ -907,14 +994,16 @@ psStats *stats) { - int i = 0; - int countInt = 0; - float countFloat = 0.0; - float mean = 0.0; - float diff = 0.0; - float sumSquares = 0.0; - float sumDiffs = 0.0; - float rangeMin = 0.0; - float rangeMax = 0.0; - + 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); @@ -979,4 +1068,16 @@ } +/****************************************************************************** +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, @@ -984,10 +1085,13 @@ psStats *stats) { - int i = 0; - int j = 0; - float clippedMean = 0.0; - float clippedStdev = 0.0; - psVector *tmpMask = NULL; - + 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))) { @@ -996,4 +1100,5 @@ } + // Endure that stats->clipSigma is within the proper range. if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && (stats->clipSigma <= CLIPPED_SIGMA_UB))) { @@ -1002,8 +1107,12 @@ } - + // 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++) { @@ -1026,9 +1135,11 @@ 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++) { - // printf("p_psVectorClippedStats(): Iteration %d (%f %f)\n", i, - // clippedMean, clippedStdev); for (j=0;jn;j++) { // a) Exclude all values x_i for which |x_i - x| > K * stdev @@ -1038,6 +1149,4 @@ } // b) compute new mean and stdev - // GUS: I should probably create a new struct here since the - // following calls will overwrite any old values in sampleMean. p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); @@ -1050,4 +1159,6 @@ } } + stats->sampleMean = oldStanMean; + stats->sampleStdev= oldStanStdev; // 7. The last calcuated value of x is the cliped mean. @@ -1066,9 +1177,19 @@ /****************************************************************************** +psVectorStats(myVector, maskVector, maskVal, stats): this is the public API +function which calls the above private stats functions based on what bits +were set in stats->options. - NOTE: The current strategy is to implement everything assuming that all - input data is of type PS_TYPE_F32. Once the basic code is in place, - we will macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64. +Inputs + myVector + maskVector + maskVal + stats +Returns + The stats structure. +NOTE: The current strategy is to implement everything assuming that all +input data is of type PS_TYPE_F32. Once the basic code is in place, we will +macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64. *****************************************************************************/ psStats *psVectorStats(psStats *stats, @@ -1077,4 +1198,5 @@ unsigned int maskVal) { + // GUS: Verify that this is the correct action. if (in == NULL) { return(stats); @@ -1084,4 +1206,6 @@ } + // 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__, @@ -1090,4 +1214,5 @@ } + // Ensure that the mask vector is of the proper size and type. if (mask != NULL) { if (in->n != mask->n) { @@ -1124,4 +1249,7 @@ } + // Since the various robust stats quantities share much computation, they + // are grouped together in a single private function: + // p_psVectorRobustStats() if ((stats->options & PS_STAT_ROBUST_MEAN) || (stats->options & PS_STAT_ROBUST_MEDIAN) || @@ -1147,8 +1275,4 @@ } - // OLD CODE: Should we check for an unknown option? - // default: - // psAbort(__func__, "Unknown options 0x%x.\n", stats->options); - return(stats); } Index: /trunk/psLib/src/math/psStats.h =================================================================== --- /trunk/psLib/src/math/psStats.h (revision 1019) +++ /trunk/psLib/src/math/psStats.h (revision 1020) @@ -1,11 +1,20 @@ +/** @file psStats.h + * \brief basic statistical operations + * @ingroup Stats + * + * This file will hold the definition of the histogram and stats data + * structures. It also contains prototypes for procedures which operate + * on those data structures. + * + * @author George Gusciora, MHPCC + * + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-14 19:33:09 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ #if !defined(PS_STATS_H) #define PS_STATS_H -/****************************************************************************** - This file will histogram/statistical functions and data structures. - *****************************************************************************/ -/** \file psStats.h - * \brief basic statistical operations - * \ingroup Stats - */ + #include "psVector.h" @@ -36,5 +45,7 @@ -/** This is the generic statistics structure */ +/** This is the generic statistics structure. It contails the data members + for the various statistic values. It contains the options member to + specifiy which statistics must be calculated. */ typedef struct { @@ -74,8 +85,8 @@ ); -/** Constructor */ +/** A constructor for the stats structure.*/ psStats *psStatsAlloc(psStatsOptions options); ///< Statistics to measure -/** Destructor */ +/** A destructor for the stats structure.*/ void psStatsFree(psStats *restrict stats); ///< Stats structure to destroy @@ -83,4 +94,5 @@ Histogram functions and data structures. *****************************************************************************/ +/** The basic histogram structure which contains bounds and bins. */ typedef struct { @@ -111,6 +123,6 @@ /** Calculate a histogram \ingroup MathGroup **/ psHistogram *psHistogramVector (psHistogram *out, ///< Histogram data - psVector *in, ///< Vector to analyse - psVector *mask, ///< Mask dat for input vector + const psVector *restrict in, ///< Vector to analyse + const psVector *restrict mask, ///< Mask dat for input vector int maskVal); ///< Mask value