Changeset 729 for trunk/psLib/src/math/psStats.c

Timestamp:

May 18, 2004, 3:50:57 PM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/math/psStats.c (modified) (7 diffs)

trunk/psLib/src/math/psStats.c

@@ -14 +14 @@
 #include "psSort.h"
 
+#include "float.h"
+#include <math.h>
+#define ROBUST_SIZE_THRESHOLD 10000   // 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 CLIPPED_NUM_ITER_LB 1
+#define CLIPPED_NUM_ITER_UB 10
+#define CLIPPED_SIGMA_LB 1.0
+#define CLIPPED_SIGMA_UB 10.0
 
 /******************************************************************************
@@ -78 +89 @@
     newHist->minNum = 0;
     newHist->maxNum = 0;
-
+    newHist->numBins = numBins;
     return(newHist);
 }
@@ -156 +167 @@
 }
 
+void p_printVector(psVector *myVector,
+                   psVector *maskVector,
+                   unsigned int maskVal,
+                   psStats *newStruct)
+{
+    int i = 0;
+
+    for (i=0;i<myVector->n;i++) {
+        if (maskVector != NULL)
+            printf("Element %d is %f (mask is %d)\n", i, myVector->vec.f[i], maskVector->vec.ui8[i]);
+        else
+            printf("Element %d is %f\n", i, myVector->vec.f[i]);
+    }
+}
+
+
 /******************************************************************************
     MISC STATISTICAL FUNCTIONS
@@ -262 +289 @@
 
 
-#define MEDIAN_SIZE_THRESHOLD 10000
 void p_psArraySampleMedian(const psVector *restrict myVector,
                            const psVector *restrict maskVector,
@@ -278 +304 @@
     }
     unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    unsortedVector->n = unsortedVector->nalloc;
     sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    sortedVector->n = sortedVector->nalloc;
 
     if (maskVector != NULL) {
         for (i=0;i<myVector->n;i++) {
             if (!(maskVal & maskVector->vec.ui8[i])) {
-                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+                unsortedVector->vec.f[count++] = myVector->vec.f[i];
+
             }
         }
@@ -303 +332 @@
 }
 
-void p_psArrayRobustMedian(const psVector *restrict myVector,
-                           const psVector *restrict maskVector,
-                           unsigned int maskVal,
-                           psStats *newStruct)
+void p_psArraysmoothHistGaussian(psHistogram *robustHistogram,
+                                 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.
+
+    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 * PI * sigma * sigma);
+            #else
+
+            tmpf = sqrtf(2.0f * 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 instead
+             * use boxcar smoothing.
+             */
+            gaussianCoefs[i] = 1.0f / (1.0f + (2.0f * (float) GAUSS_WIDTH));
+        }
+    }
+
+    for(i=0;i<robustHistogram->numBins;i++) {
+        for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) {
+            if (((j+i) >= 0) && ((j+i) < robustHistogram->numBins)) {
+                robustHistogram->nums->vec.i32[j+i]+=
+                    (gaussianCoefs[j+GAUSS_WIDTH] *
+                     (float) robustHistogram->nums->vec.i32[j+i]);
+            }
+        }
+    }
+}
+
+/******************************************************************************
+    p_psArraySampleQuartiles()
+ This procedure calculates the upper and/or lower quartiles of the
+ data set.  It is assumed that the data set is small enough that we
+ can sort all the data points and calculate the quartiles exactly.
+ *****************************************************************************/
+void p_psArraySampleQuartiles(const psVector *restrict myVector,
+                              const psVector *restrict maskVector,
+                              unsigned int maskVal,
+                              psStats *newStruct)
+{
+    psVector *unsortedVector = NULL;
+    psVector *sortedVector = NULL;
+    int count = 0;
+    int ind = 0;
+    int i = 0;
+
+    // return is we have already calculated both quartile points.
+    if ((!isnan(newStruct->sampleLQ)) &&
+            (!isnan(newStruct->sampleUQ))) {
+        return;
+    }
+
+    if (-1 == newStruct->nValues) {
+        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
+    }
+
+
+    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+
+    count = 0;
+    if (maskVector != NULL) {
+        for (i=0;i<myVector->n;i++) {
+            if (!(maskVal & maskVector->vec.ui8[i])) {
+                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+            }
+        }
+        psSort(sortedVector, unsortedVector);
+    } else {
+        psSort(sortedVector, myVector);
+    }
+
+    if (newStruct->options & PS_STAT_SAMPLE_LQ) {
+        ind = 3 * (newStruct->nValues / 4);
+        newStruct->sampleUQ = sortedVector->vec.f[ind];
+    }
+
+    if (newStruct->options & PS_STAT_SAMPLE_UQ) {
+        ind = (newStruct->nValues / 4);
+        newStruct->sampleLQ = sortedVector->vec.f[ind];
+    }
+
+    psVectorFree(unsortedVector);
+    psVectorFree(sortedVector);
+}
+
+
+/******************************************************************************
+    p_psArrayRobustStats(): this procedure calcualtes a variety of robust
+    stat measures:
+ PS_STAT_ROBUST_MEAN
+        PS_STAT_ROBUST_MEAN_NVALUES
+        PS_STAT_ROBUST_MEDIAN
+        PS_STAT_ROBUST_MEDIAN_NVALUES
+        PS_STAT_ROBUST_MODE
+        PS_STAT_ROBUST_MODE_NVALUES
+        PS_STAT_ROBUST_STDEV
+    I have included all that computation in a single function, as opposed to
+    breaking it across several functions for one primary reason: the 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.
+ *****************************************************************************/
+void p_psArrayRobustStats(const psVector *restrict myVector,
+                          const psVector *restrict maskVector,
+                          unsigned int maskVal,
+                          psStats *newStruct)
 {
     psHistogram *robustHistogram = NULL;
     float binSize = 0.0;
-
-    //    if (isnan(myVector->robustLQ) ||
-    //        isnan(myVector->robustUQ)) {
-    //        p_psArrayRobustQuartiles(myVector, maskVector, maskVal, newStruct);
-    //    }
-    //    binSize = ((myVector->robustUQ - myVector->robustLQ) / 1.34) / 10.0;
-
+    float sigmaE = 0.0;
+    int LQBinNum = -1;
+    int UQBinNum = -1;
+    int i = 0;    // Loop index variable.
+    int maxBinNum = 0;
+    int maxBinCount = 0;
+    float dL = 0.0;
+
+    // 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 (0.0 == newStruct->sampleUQ) {
+        newStruct->options = newStruct->options | PS_STAT_SAMPLE_UQ;
+        p_psArraySampleQuartiles(myVector,
+                                 maskVector,
+                                 maskVal,
+                                 newStruct);
+    }
+
+    if (isnan(newStruct->sampleLQ)) {
+        newStruct->options = newStruct->options | PS_STAT_SAMPLE_LQ;
+        p_psArraySampleQuartiles(myVector,
+                                 maskVector,
+                                 maskVal,
+                                 newStruct);
+    }
+
+    // Compute the initial bin size of the robust histogram.
+    sigmaE = (newStruct->sampleUQ - newStruct->sampleLQ) / 1.34f;
+    binSize = sigmaE / 10.0f;
+
+    // Detemine minimum and maximum values in the data vector.
+    if (isnan(newStruct->min)) {
+        p_psArrayMin(myVector, maskVector, maskVal, newStruct);
+    }
+    if (isnan(newStruct->max)) {
+        p_psArrayMax(myVector, maskVector, maskVal, newStruct);
+    }
+
+    // Create the histogram structure.
     robustHistogram = psHistogramAlloc(newStruct->min,
                                        newStruct->max,
                                        binSize);
-    //    p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
-    //    robustHistogram = psGetArrayHistogram(robustHistogram, myVector);
-    //    p_psArraySmooth(robustHistogram, (binSize / 4.0));
-    //    dL = (myVector->robustUQ - myVector->robustLQ) / 8.0;
-
-
-    // BROAD: Calculate the Robust Median
-    // Determine the LQ of the distribution.
-    // Determine the UQ of the distribution.
-    // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
-    // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
-
-    // Find the bin with the peak value between LQ and UQ (the MODE)
-    // dL = (UQ - LQ) / 8
+    // Populate the histogram arrat.
+    robustHistogram = psGetArrayHistogram(robustHistogram, myVector);
+
+    // Smooth the histogram.
+    p_psArraysmoothHistGaussian(robustHistogram, sigmaE/4.0f);
+
+    LQBinNum = -1;
+    UQBinNum = -1;
+    for (i=0;i<robustHistogram->numBins;i++) {
+        if ((robustHistogram->nums->vec.i32[i] <= newStruct->sampleLQ) &&
+                (newStruct->sampleLQ <= robustHistogram->nums->vec.i32[i])) {
+            LQBinNum = i;
+        }
+
+        if ((robustHistogram->nums->vec.i32[i] <= newStruct->sampleUQ) &&
+                (newStruct->sampleUQ <= robustHistogram->nums->vec.i32[i])) {
+            UQBinNum = i;
+        }
+    }
+
+    // Determine the bin with the peak value in the range LQ to UQ.
+    maxBinNum = LQBinNum;
+    maxBinCount = robustHistogram->nums->vec.i32[maxBinNum];
+    for (i=LQBinNum;i<=UQBinNum;i++) {
+        if (robustHistogram->nums->vec.i32[i] > maxBinCount) {
+            maxBinNum = i;
+            maxBinCount = robustHistogram->nums->vec.i32[i];
+        }
+    }
+
+    dL = (newStruct->robustUQ - newStruct->robustLQ) / 8.0;
+
     // 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
-}
-
-
-float p_psArrayXXX(const psVector *restrict myVector,
-                   const psVector *restrict maskVector,
-                   unsigned int maskVal)
-{
-    printf("ERROR: Don't call me: p_psArrayXXX()\n");
-    exit(1);
+    // What is the mean_r?
+
+    if (newStruct->options & PS_STAT_ROBUST_MEAN) {
+        newStruct->robustMean = 0.0;
+    }
+    if  (newStruct->options & PS_STAT_ROBUST_MEAN_NVALUES) {
+        newStruct->robustMeanNvalues = 0.0;
+    }
+    if  (newStruct->options & PS_STAT_ROBUST_MEDIAN) {
+        newStruct->robustMedian = 0.0;
+    }
+    if  (newStruct->options & PS_STAT_ROBUST_MEDIAN_NVALUES) {
+        newStruct->robustMedianNvalues = 0.0;
+    }
+    if  (newStruct->options & PS_STAT_ROBUST_MODE) {
+        newStruct->robustMode = maxBinNum;
+    }
+    if  (newStruct->options & PS_STAT_ROBUST_MODE_NVALUES) {
+        newStruct->robustModeNvalues = 0.0;
+    }
+    if  (newStruct->options & PS_STAT_ROBUST_STDEV) {
+        newStruct->robustStdev = 0.0;
+    }
+    if  (newStruct->options & PS_STAT_ROBUST_UQ) {
+        newStruct->robustUQ = 0.0;
+    }
+    if  (newStruct->options & PS_STAT_ROBUST_LQ) {
+        newStruct->robustLQ = 0.0;
+    }
 }
 
@@ -385 +610 @@
     }
     countFloat = (float) countInt;
+
     #ifdef DARWIN
 
-    newStruct->sampleStdev = sqrt(
-                                 #else
-                                 newStruct->sampleStdev = sqrtf(
-                                                              #endif
-                                                              (sumSquares-(sumDiffs * sumDiffs/countFloat))/ (countFloat-1));
-                             }
-
-                             /******************************************************************************
-                              *****************************************************************************/
-                             void p_psArraySampleQuartiles(const psVector *restrict myVector,
-                                                           const psVector *restrict maskVector,
-                                                           unsigned int maskVal,
-                                                           psStats *newStruct)
-                             {
-                                 psVector *unsortedVector = NULL;
-                                 psVector *sortedVector = NULL;
-                                 int count = 0;
-                                 int ind = 0;
-                                 int i = 0;
-
-                                 // return is we have already calculated both quartile points.
-                                 if ((!isnan(newStruct->sampleLQ)) &&
-                                         (!isnan(newStruct->sampleUQ))) {
-                                     return;
-                                 }
-
-                                 if (-1 == newStruct->nValues) {
-                                     p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
-                                 }
-
-
-                                 unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
-                                 sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
-
-                                 count = 0;
-                                 if (maskVector != NULL) {
-                                     for (i=0;i<myVector->n;i++) {
-                                         if (!(maskVal & maskVector->vec.ui8[i])) {
-                                             unsortedVector->vec.f[count++] = maskVector->vec.f[i];
-                                         }
-                                     }
-                                     psSort(sortedVector, unsortedVector);
-                                 } else {
-                                     psSort(sortedVector, myVector);
-                                 }
-
-                                 ind = 3 * (newStruct->nValues / 4);
-                                 newStruct->sampleUQ = sortedVector->vec.f[ind];
-                                 ind = (newStruct->nValues / 4);
-                                 newStruct->sampleLQ = sortedVector->vec.f[ind];
-
-                                 psVectorFree(unsortedVector);
-                                 psVectorFree(sortedVector);
-                             }
-
-                             /******************************************************************************
-                              
-                                 NOTE: The current strategy is to implement everything assuming that all
-                                 input data is of type PS_TYPE_FLOAT.  Once the basic code is in place,
-                                 we will macro-ize everything and add PS_TYPE_UINT16 and PS_TYPE_DOUBLE.
-                              
-                              *****************************************************************************/
-                             psStats *psArrayStats(const psVector *restrict myVector,
-                                                   const psVector *restrict maskVector,
-                                                   unsigned int maskVal,
-                                                   psStats *stats)
-                             {
-                                 psStats *newStruct = NULL;
-
-                                 if (myVector == NULL) {
-                                     psAbort(__func__,
-                                             "Input data array (myVector) was NULL.");
-                                 }
-
-                                 if (myVector->type.type != PS_TYPE_FLOAT) {
-                                     psAbort(__func__,
-                                             "Only data type PS_TYPE_FLOAT is currently supported.");
-                                 }
-
-                                 if (maskVector != NULL) {
-                                     if (myVector->n != maskVector->n) {
-                                         psAbort(__func__,
-                                                 "Vector data and vector mask are of different sizes.");
-                                     }
-                                     if (maskVector->type.type != PS_TYPE_UINT8) {
-                                         psAbort(__func__,
-                                                 "Vector mask must be type PS_TYPE_UINT8");
-                                     }
-                                 }
-                                 newStruct = psStatsAlloc(stats->options);
-
-                                 // ************************************************************************
-                                 if (stats->options & PS_STAT_SAMPLE_MEAN) {
-                                     p_psArraySampleMean(myVector, maskVector, maskVal, newStruct);
-                                 }
-
-                                 // ************************************************************************
-                                 if (stats->options & PS_STAT_MAX) {
-                                     p_psArrayMax(myVector, maskVector, maskVal, newStruct);
-                                 }
-
-                                 // ************************************************************************
-                                 if (stats->options & PS_STAT_MIN) {
-                                     p_psArrayMin(myVector, maskVector, maskVal, newStruct);
-                                 }
-
-                                 // ************************************************************************
-                                 if (stats->options & PS_STAT_NVALUES) {
-                                     p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
-                                 }
-
-                                 // ************************************************************************
-                                 if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
-                                     p_psArraySampleMedian(myVector, maskVector, maskVal, newStruct);
-                                 }
-
-                                 // ************************************************************************
-                                 if (stats->options & PS_STAT_SAMPLE_STDEV) {
-                                     p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
-                                 }
-
-                                 // ************************************************************************
-                                 if ((stats->options & PS_STAT_SAMPLE_UQ) ||
-                                         (stats->options & PS_STAT_SAMPLE_LQ)) {
-                                     p_psArraySampleQuartiles(myVector, maskVector, maskVal, newStruct);
-                                 }
-
-
-
-
-                                 if (stats->options & PS_STAT_ROBUST_MEAN) {
-                                     newStruct->robustMean = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_ROBUST_MEAN_NVALUES) {
-                                     newStruct->robustMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_ROBUST_MEDIAN) {
-                                     newStruct->robustMedian = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_ROBUST_MEDIAN_NVALUES) {
-                                     newStruct->robustMedianNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_ROBUST_MODE) {
-                                     newStruct->robustMode = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_ROBUST_MODE_NVALUES) {
-                                     newStruct->robustModeNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_ROBUST_STDEV) {
-                                     newStruct->robustStdev = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if ((stats->options & PS_STAT_ROBUST_UQ) ||
-                                         (stats->options & PS_STAT_ROBUST_LQ)) {
-                                     newStruct->robustLQ = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_CLIPPED_MEAN) {
-                                     newStruct->clippedMean = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_CLIPPED_MEAN_NVALUES) {
-                                     newStruct->clippedMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_CLIPPED_MEAN_NSIGMA) {
-                                     newStruct->clipSigma = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-                                 if (stats->options & PS_STAT_CLIPPED_STDEV) {
-                                     newStruct->clippedStdev = p_psArrayXXX(myVector, maskVector, maskVal);
-                                 }
-
-
-                                 //    OLD CODE: Should we check for an unknown option?
-                                 //    default:
-                                 //        psAbort(__func__, "Unknown options 0x%x.\n", stats->options);
-
-                                 return(newStruct);
-                             }
+    newStruct->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs *
+                                            sumDiffs/countFloat))/ (countFloat-1));
+    #else
+
+    newStruct->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
+                                     sumDiffs/countFloat))/ (countFloat-1));
+    #endif
+}
+
+void p_psArrayClippedStats(const psVector *restrict myVector,
+                           const psVector *restrict maskVector,
+                           unsigned int maskVal,
+                           psStats *newStruct)
+{
+    int i = 0;
+    int j = 0;
+    float clippedMean = 0.0;
+    float clippedStdev = 0.0;
+    psVector *tmpMask = NULL;
+
+    if (!((CLIPPED_NUM_ITER_LB <= newStruct->clipIter ) &&
+            (newStruct->clipIter <= CLIPPED_NUM_ITER_UB))) {
+        psAbort(__func__, "Unallowed value for clipIter (%d).\n",
+                newStruct->clipIter);
+    }
+
+    if (!((CLIPPED_SIGMA_LB <= newStruct->clipSigma ) &&
+            (newStruct->clipSigma <= CLIPPED_SIGMA_UB))) {
+        psAbort(__func__, "Unallowed value for clipSigma (%f).\n",
+                newStruct->clipSigma);
+    }
+
+    tmpMask = psVectorAlloc(maskVector->type.type, myVector->nalloc);
+
+    tmpMask->n = maskVector->n;
+    for (i=0;i<tmpMask->n;i++) {
+        tmpMask->vec.ui8[i] = maskVector->vec.ui8[i];
+    }
+
+    // 1. Compute the sample median.
+    p_psArraySampleMedian(myVector, maskVector, maskVal, newStruct);
+
+    // 2. Compute the sample standard deviation.
+    p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
+
+    // 3. Use the sample median as the first estimator of the mean X.
+    clippedMean = newStruct->sampleMean;
+
+    // 4. Use the sample stdev as the first estimator of the mean stdev.
+    clippedStdev = newStruct->sampleStdev;
+
+    // 5. Repeat N times:
+
+    for (i=0;i<newStruct->clipIter;i++) {
+        for (j=0;j<myVector->n;j++) {
+            // a) Exclude all values x_i for which |x_i - x| > K * stdev
+            if ( fabs(myVector->vec.f[j] - clippedMean) >
+                    (newStruct->clipSigma * clippedStdev)) {
+                tmpMask->vec.ui8[i] = 0xff;
+            }
+            // 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_psArraySampleMedian(myVector, tmpMask, maskVal, newStruct);
+            p_psArraySampleStdev(myVector, tmpMask, maskVal, newStruct);
+
+            // c) Use the new mean for x
+            clippedMean = newStruct->sampleMean;
+
+            // d) Use the new stdev for stdev
+            clippedStdev = newStruct->sampleStdev;
+        }
+    }
+
+    // 7. The last calcuated value of x is the cliped mean.
+    if (newStruct->options & PS_STAT_CLIPPED_MEAN) {
+        newStruct->clippedMean = clippedMean;
+    }
+
+    // 8. The last calcuated value of stdev is the cliped stdev.
+    if (newStruct->options & PS_STAT_CLIPPED_STDEV) {
+        newStruct->clippedStdev = clippedStdev;
+    }
+
+    if (newStruct->options & PS_STAT_CLIPPED_MEAN_NVALUES) {
+        p_psArrayNValues(myVector, tmpMask, maskVal, newStruct);
+    }
+
+    if (newStruct->options & PS_STAT_CLIPPED_MEAN_NSIGMA) {
+        // GUS: What to do here?
+    }
+
+    psVectorFree(tmpMask);
+}
+
+
+/******************************************************************************
+ 
+    NOTE: The current strategy is to implement everything assuming that all
+    input data is of type PS_TYPE_FLOAT.  Once the basic code is in place,
+    we will macro-ize everything and add PS_TYPE_UINT16 and PS_TYPE_DOUBLE.
+ 
+ *****************************************************************************/
+psStats *psArrayStats(const psVector *restrict myVector,
+                      const psVector *restrict maskVector,
+                      unsigned int maskVal,
+                      psStats *stats)
+{
+    psStats *newStruct = NULL;
+
+    if (myVector == NULL) {
+        psAbort(__func__,
+                "Input data array (myVector) was NULL.");
+    }
+
+    if (myVector->type.type != PS_TYPE_FLOAT) {
+        psAbort(__func__,
+                "Only data type PS_TYPE_FLOAT is currently supported.");
+    }
+    if (maskVector != NULL) {
+        if (myVector->n != maskVector->n) {
+            psAbort(__func__,
+                    "Vector data and vector mask are of different sizes.");
+        }
+        if (maskVector->type.type != PS_TYPE_UINT8) {
+            psAbort(__func__,
+                    "Vector mask must be type PS_TYPE_UINT8");
+        }
+    }
+    newStruct = psStatsAlloc(stats->options);
+
+    // ************************************************************************
+    if (stats->options & PS_STAT_SAMPLE_MEAN) {
+        p_psArraySampleMean(myVector, maskVector, maskVal, newStruct);
+    }
+
+    // ************************************************************************
+    if (stats->options & PS_STAT_MAX) {
+        p_psArrayMax(myVector, maskVector, maskVal, newStruct);
+    }
+
+    // ************************************************************************
+    if (stats->options & PS_STAT_MIN) {
+        p_psArrayMin(myVector, maskVector, maskVal, newStruct);
+    }
+
+    // ************************************************************************
+    if (stats->options & PS_STAT_NVALUES) {
+        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
+    }
+
+    // ************************************************************************
+    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
+        p_psArraySampleMedian(myVector, maskVector, maskVal, newStruct);
+    }
+
+    // ************************************************************************
+    if (stats->options & PS_STAT_SAMPLE_STDEV) {
+        p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
+    }
+
+    // ************************************************************************
+    if ((stats->options & PS_STAT_SAMPLE_UQ) ||
+            (stats->options & PS_STAT_SAMPLE_LQ)) {
+        p_psArraySampleQuartiles(myVector, maskVector, maskVal, newStruct);
+    }
+
+    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
+            (stats->options & PS_STAT_ROBUST_MEAN_NVALUES) ||
+            (stats->options & PS_STAT_ROBUST_MEDIAN) ||
+            (stats->options & PS_STAT_ROBUST_MEDIAN_NVALUES) ||
+            (stats->options & PS_STAT_ROBUST_MODE) ||
+            (stats->options & PS_STAT_ROBUST_MODE_NVALUES) ||
+            (stats->options & PS_STAT_ROBUST_STDEV) ||
+            (stats->options & PS_STAT_ROBUST_UQ) ||
+            (stats->options & PS_STAT_ROBUST_LQ)) {
+        p_psArrayClippedStats(myVector, maskVector, maskVal, newStruct);
+    }
+
+    if ((stats->options & PS_STAT_CLIPPED_MEAN) ||
+            (stats->options & PS_STAT_CLIPPED_MEAN_NVALUES) ||
+            (stats->options & PS_STAT_CLIPPED_MEAN_NSIGMA) ||
+            (stats->options & PS_STAT_CLIPPED_STDEV)) {
+        p_psArrayClippedStats(myVector, maskVector, maskVal, newStruct);
+    }
+
+    //    OLD CODE: Should we check for an unknown option?
+    //    default:
+    //        psAbort(__func__, "Unknown options 0x%x.\n", stats->options);
+
+    return(newStruct);
+}