Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 729)
+++ trunk/psLib/src/math/psStats.c	(revision 796)
@@ -15,5 +15,4 @@
 
 #include "float.h"
-#include <math.h>
 #define ROBUST_SIZE_THRESHOLD 10000   // Vectors that are large than this
 // will use robust statistical methods.
@@ -25,5 +24,7 @@
 #define CLIPPED_SIGMA_LB 1.0
 #define CLIPPED_SIGMA_UB 10.0
-
+#define true 1
+#define false 0
+#define MYMAXFLOAT 1e99
 /******************************************************************************
     psStatsAlloc(): This routine must create a new psStats data structure.
@@ -39,15 +40,14 @@
     newStruct->sampleUQ = NAN;
     newStruct->sampleLQ = NAN;
+    newStruct->sampleLimit = NAN;
     newStruct->robustMean = NAN;
-    newStruct->robustMeanNvalues = -1;
     newStruct->robustMedian = NAN;
-    newStruct->robustMedianNvalues = -1;
     newStruct->robustMode = NAN;
-    newStruct->robustModeNvalues = -1;
     newStruct->robustStdev = NAN;
     newStruct->robustUQ = NAN;
     newStruct->robustLQ = NAN;
+    newStruct->robustN50 = NAN;
+    newStruct->robustNfit = NAN;
     newStruct->clippedMean = NAN;
-    newStruct->clippedMeanNvalues = -1;
     newStruct->clippedStdev = NAN;
     newStruct->clipSigma = NAN;
@@ -55,5 +55,5 @@
     newStruct->min = NAN;
     newStruct->max = NAN;
-    newStruct->nValues = -1;
+    newStruct->binsize = NAN;
     newStruct->options = options;
 
@@ -74,48 +74,40 @@
 psHistogram *psHistogramAlloc(float lower,
                               float upper,
-                              float size )
-{
+                              int n)
+{
+    int i = 0;
     psHistogram *newHist = NULL;
-    int numBins = 0;
-
-
-    numBins = 1 + ((upper - lower) / size);
+    float binSize = 0.0;
+
     newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
-    newHist->lower = psVectorAlloc(PS_TYPE_FLOAT, numBins);
-    newHist->upper = psVectorAlloc(PS_TYPE_FLOAT, numBins);
-    newHist->nums =  psVectorAlloc(PS_TYPE_INT32, numBins);
-    newHist->minVal = lower;
-    newHist->maxVal = upper;
+    newHist->bounds = psVectorAlloc(PS_TYPE_FLOAT, n+1);
+    binSize = (upper - lower) / (float) n;
+    for (i=0;i<n+1;i++) {
+        newHist->bounds->vec.f[i] = lower + (binSize * (float) i);
+    }
+    newHist->nums =  psVectorAlloc(PS_TYPE_INT32, n);
     newHist->minNum = 0;
     newHist->maxNum = 0;
-    newHist->numBins = numBins;
+    newHist->uniform = true;
+
     return(newHist);
 }
 
-psHistogram *psHistogramAllocGeneric(const psVector *restrict lower,
-                                     const psVector *restrict upper,
-                                     float minVal,
-                                     float maxVal)
+
+psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds)
 {
     psHistogram *newHist = NULL;
-    int numBins = 0;
     int i;
 
-    if (lower->n != upper->n) {
-        psError(__func__, "There is a different number of upper/lower values");
-    }
-    numBins = lower->n;
     newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
-    newHist->lower = psVectorAlloc(PS_TYPE_FLOAT, numBins);
-    newHist->upper = psVectorAlloc(PS_TYPE_FLOAT, numBins);
-    for (i=0;i<numBins;i++) {
-        newHist->lower->vec.f[i] = lower->vec.f[i];
-        newHist->upper->vec.f[i] = upper->vec.f[i];
-    }
-    newHist->nums = psVectorAlloc(PS_TYPE_INT32, numBins);
-    newHist->minVal = minVal;
-    newHist->maxVal = maxVal;
+    newHist->bounds = psVectorAlloc(PS_TYPE_FLOAT, bounds->n);
+    for (i=0;i<bounds->n;i++) {
+        newHist->bounds->vec.f[i] = bounds->vec.f[i];
+    }
+    newHist->nums = psVectorAlloc(PS_TYPE_INT32, (bounds->n)-1);
+
     newHist->minNum = 0;
     newHist->maxNum = 0;
+    newHist->uniform = false;
 
     return(newHist);
@@ -141,28 +133,53 @@
         myHist->maxNum
  *****************************************************************************/
-psHistogram *psGetArrayHistogram(psHistogram *restrict myHist,
-                                 const psVector *restrict myVector)
-{
-    int i = 0;
+psHistogram *psHistogramVector(psHistogram *out,
+                               psVector *in,
+                               psVector *mask,
+                               int maskVal)
+{
+    int i = 0;
+    int j = 0;
     float binSize = 0.0;
     int binNum = 0;
-
-    binSize = (myHist->maxVal - myHist->minVal) / (float) myHist->nums->n;
-    for (i=0;i<myVector->n;i++) {
-        if (myVector->vec.f[i] < myHist->minVal) {
-            myHist->minNum++;
-        } else if (myVector->vec.f[i] > myHist->maxVal) {
-            myHist->maxNum++;
-        } else {
-            binNum = (int) ((myVector->vec.f[i] - myHist->minVal) / binSize);
-            if ((myVector->vec.f[i] >= myHist->lower->vec.f[i]) &&
-                    (myVector->vec.f[i] <= myHist->upper->vec.f[i])) {
-                myHist->nums->vec.i32[i]++;
+    int numBins = 0;
+
+    numBins = out->nums->n;
+    for (i=0;i<in->n;i++) {
+        // Check if this pixel is masked, and if so, skip it.
+        if (!(mask->vec.i32[i] & maskVal)) {
+            // Check if this pixel is below the minimum value, and if so
+            // count it, then skip it.
+            if (in->vec.f[i] < out->bounds->vec.f[0]) {
+                out->minNum++;
+
+                // Check if this pixel is above the maximum value, and if so
+                // count it, then skip it.
+            } else if (in->vec.f[i] > out->bounds->vec.f[numBins]) {
+                out->maxNum++;
             } else {
-                psError(__func__, "data value was not within the bounds of the bin it should habe been in.");
-            }
-        }
-    }
-    return(myHist);
+                // If this is a uniform histogram, determining the correct
+                // number is trivial.
+                if (out->uniform == true) {
+                    binSize = out->bounds->vec.f[1] - out->bounds->vec.f[0];
+
+                    binNum = (int) ((in->vec.f[i] - out->bounds->vec.f[0]) /
+                                    binSize);
+                    (out->nums->vec.i32[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->vec.i32[j] <= in->vec.f[i]) &&
+                                (in->vec.f[i] <= out->bounds->vec.i32[j+1])) {
+                            (out->nums->vec.i32[j])++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    return(out);
 }
 
@@ -184,43 +201,107 @@
 
 /******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
     MISC STATISTICAL FUNCTIONS
- *****************************************************************************/
-void p_psArraySampleMean(const psVector *restrict myVector,
-                         const psVector *restrict maskVector,
-                         unsigned int maskVal,
-                         psStats *newStruct)
+ ******************************************************************************
+ ******************************************************************************
+ *****************************************************************************/
+void p_psVectorSampleMean(const psVector *restrict myVector,
+                          const psVector *restrict maskVector,
+                          unsigned int maskVal,
+                          psStats *stats)
 {
     int i = 0;
     float mean = 0.0;
     int count = 0;
-
-    if (maskVector != NULL) {
-        for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.ui8[i])) {
+    float rangeMin = 0.0;
+    float rangeMax = 0.0;
+
+    if (stats->options & PS_STAT_USE_RANGE) {
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i]) &&
+                        (rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    mean+= myVector->vec.f[i];
+                    count++;
+                }
+            }
+            mean/= (float) count;
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                if ((rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    mean+= myVector->vec.f[i];
+                    count++;
+                }
+            }
+            mean/= (float) count;
+        }
+    } else {
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i])) {
+                    mean+= myVector->vec.f[i];
+                    count++;
+                }
+            }
+            mean/= (float) count;
+        } else {
+            for (i=0;i<myVector->n;i++) {
                 mean+= myVector->vec.f[i];
-                count++;
-            }
-        }
-        mean/= (float) count;
+            }
+            mean/= (float) myVector->n;
+        }
+    }
+
+    stats->sampleMean = mean;
+}
+
+void p_psVectorMax(const psVector *restrict myVector,
+                   const psVector *restrict maskVector,
+                   unsigned int maskVal,
+                   psStats *stats)
+{
+    int i = 0;
+    float max = -MYMAXFLOAT;
+    float rangeMin = 0.0;
+    float rangeMax = 0.0;
+
+    if (stats->options & PS_STAT_USE_RANGE) {
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i])) {
+                    if ((myVector->vec.f[i] > max) &&
+                            (rangeMin <= myVector->vec.f[i]) &&
+                            (myVector->vec.f[i] <= rangeMax)) {
+                        max = myVector->vec.f[i];
+                    }
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                if ((myVector->vec.f[i] > max) &&
+                        (rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    max = myVector->vec.f[i];
+                }
+            }
+        }
     } else {
-        for (i=0;i<myVector->n;i++) {
-            mean+= myVector->vec.f[i];
-        }
-        mean/= (float) myVector->n;
-    }
-    newStruct->sampleMean = mean;
-}
-
-void p_psArrayMax(const psVector *restrict myVector,
-                  const psVector *restrict maskVector,
-                  unsigned int maskVal,
-                  psStats *newStruct)
-{
-    int i = 0;
-    float max = -1e99;
-
-    if (maskVector != NULL) {
-        for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.ui8[i])) {
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i])) {
+                    if (myVector->vec.f[i] > max) {
+                        max = myVector->vec.f[i];
+                    }
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
                 if (myVector->vec.f[i] > max) {
                     max = myVector->vec.f[i];
@@ -228,25 +309,52 @@
             }
         }
+    }
+
+    stats->max = max;
+}
+
+void p_psVectorMin(const psVector *restrict myVector,
+                   const psVector *restrict maskVector,
+                   unsigned int maskVal,
+                   psStats *stats)
+{
+    int i = 0;
+    float min = MYMAXFLOAT;
+    float rangeMin = 0.0;
+    float rangeMax = 0.0;
+
+    if (stats->options & PS_STAT_USE_RANGE) {
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i])) {
+                    if ((myVector->vec.f[i] < min) &&
+                            (rangeMin <= myVector->vec.f[i]) &&
+                            (myVector->vec.f[i] <= rangeMax)) {
+                        min = myVector->vec.f[i];
+                    }
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                if ((myVector->vec.f[i] < min) &&
+                        (rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    min = myVector->vec.f[i];
+                }
+            }
+        }
     } else {
-        for (i=0;i<myVector->n;i++) {
-            if (myVector->vec.f[i] > max) {
-                max = myVector->vec.f[i];
-            }
-        }
-    }
-    newStruct->max = max;
-}
-
-void p_psArrayMin(const psVector *restrict myVector, /* FLOATS */
-                  const psVector *restrict maskVector, /* INTS */
-                  unsigned int maskVal,
-                  psStats *newStruct)
-{
-    int i = 0;
-    float min = 1e99;
-
-    if (maskVector != NULL) {
-        for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.ui8[i])) {
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i])) {
+                    if (myVector->vec.f[i] < min) {
+                        min = myVector->vec.f[i];
+                    }
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
                 if (myVector->vec.f[i] < min) {
                     min = myVector->vec.f[i];
@@ -254,18 +362,13 @@
             }
         }
-    } else {
-        for (i=0;i<myVector->n;i++) {
-            if (myVector->vec.f[i] < min) {
-                min = myVector->vec.f[i];
-            }
-        }
-    }
-    newStruct->min = min;
+    }
+
+    stats->min = min;
 }
 
 /******************************************************************************
- 
- *****************************************************************************/
-void p_psArrayNValues(const psVector *restrict myVector,
+ This routine calculates the number of non-masked pixels in the vector. 
+ *****************************************************************************/
+int p_psVectorNValues(const psVector *restrict myVector,
                       const psVector *restrict maskVector,
                       unsigned int maskVal,
@@ -284,13 +387,13 @@
         numData = myVector->n;
     }
-    newStruct->nValues = numData;
-}
-
-
-
-void p_psArraySampleMedian(const psVector *restrict myVector,
-                           const psVector *restrict maskVector,
-                           unsigned int maskVal,
-                           psStats *newStruct)
+    return(numData);
+}
+
+
+
+void p_psVectorSampleMedian(const psVector *restrict myVector,
+                            const psVector *restrict maskVector,
+                            unsigned int maskVal,
+                            psStats *newStruct)
 {
     psVector *unsortedVector = NULL;
@@ -299,11 +402,10 @@
     int i = 0;
     float median = 0.0;
-
-    if (-1 == newStruct->nValues) {
-        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
-    }
-    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    int nValues = 0;
+
+    nValues = p_psVectorNValues(myVector, maskVector, maskVal, newStruct);
+    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, nValues);
     unsortedVector->n = unsortedVector->nalloc;
-    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, nValues);
     sortedVector->n = sortedVector->nalloc;
 
@@ -320,9 +422,9 @@
     }
 
-    if (0 == (newStruct->nValues % 2)) {
-        median = 0.5 * (sortedVector->vec.f[(newStruct->nValues/2)-1] +
-                        sortedVector->vec.f[newStruct->nValues/2]);
+    if (0 == (nValues % 2)) {
+        median = 0.5 * (sortedVector->vec.f[(nValues/2)-1] +
+                        sortedVector->vec.f[nValues/2]);
     } else {
-        median = sortedVector->vec.f[newStruct->nValues/2];
+        median = sortedVector->vec.f[nValues/2];
     }
 
@@ -332,6 +434,6 @@
 }
 
-void p_psArraysmoothHistGaussian(psHistogram *robustHistogram,
-                                 float sigma)
+void p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
+                                  float sigma)
 {
     int i = 0;
@@ -365,7 +467,7 @@
     }
 
-    for(i=0;i<robustHistogram->numBins;i++) {
+    for(i=0;i<robustHistogram->nums->n;i++) {
         for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) {
-            if (((j+i) >= 0) && ((j+i) < robustHistogram->numBins)) {
+            if (((j+i) >= 0) && ((j+i) < robustHistogram->nums->n)) {
                 robustHistogram->nums->vec.i32[j+i]+=
                     (gaussianCoefs[j+GAUSS_WIDTH] *
@@ -377,13 +479,13 @@
 
 /******************************************************************************
-    p_psArraySampleQuartiles()
+    p_psVectorSampleQuartiles()
  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)
+void p_psVectorSampleQuartiles(const psVector *restrict myVector,
+                               const psVector *restrict maskVector,
+                               unsigned int maskVal,
+                               psStats *newStruct)
 {
     psVector *unsortedVector = NULL;
@@ -392,4 +494,5 @@
     int ind = 0;
     int i = 0;
+    int nValues = 0;
 
     // return is we have already calculated both quartile points.
@@ -399,11 +502,8 @@
     }
 
-    if (-1 == newStruct->nValues) {
-        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
-    }
-
-
-    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
-    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    nValues = p_psVectorNValues(myVector, maskVector, maskVal, newStruct);
+
+    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, nValues);
+    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, nValues);
 
     count = 0;
@@ -419,13 +519,8 @@
     }
 
-    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];
-    }
+    ind = 3 * (nValues / 4);
+    newStruct->sampleUQ = sortedVector->vec.f[ind];
+    ind = (nValues / 4);
+    newStruct->sampleLQ = sortedVector->vec.f[ind];
 
     psVectorFree(unsortedVector);
@@ -435,12 +530,9 @@
 
 /******************************************************************************
-    p_psArrayRobustStats(): this procedure calcualtes a variety of robust
+    p_psVectorRobustStats(): 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
@@ -452,8 +544,8 @@
     processing would be duplicated.
  *****************************************************************************/
-void p_psArrayRobustStats(const psVector *restrict myVector,
-                          const psVector *restrict maskVector,
-                          unsigned int maskVal,
-                          psStats *newStruct)
+void p_psVectorRobustStats(const psVector *restrict myVector,
+                           const psVector *restrict maskVector,
+                           unsigned int maskVal,
+                           psStats *newStruct)
 {
     psHistogram *robustHistogram = NULL;
@@ -474,17 +566,19 @@
 
     if (0.0 == newStruct->sampleUQ) {
-        newStruct->options = newStruct->options | PS_STAT_SAMPLE_UQ;
-        p_psArraySampleQuartiles(myVector,
-                                 maskVector,
-                                 maskVal,
-                                 newStruct);
+        // GUS: fix this
+        //        newStruct->options = newStruct->options | PS_STAT_SAMPLE_UQ;
+        p_psVectorSampleQuartiles(myVector,
+                                  maskVector,
+                                  maskVal,
+                                  newStruct);
     }
 
     if (isnan(newStruct->sampleLQ)) {
-        newStruct->options = newStruct->options | PS_STAT_SAMPLE_LQ;
-        p_psArraySampleQuartiles(myVector,
-                                 maskVector,
-                                 maskVal,
-                                 newStruct);
+        // GUS: fix this
+        //        newStruct->options = newStruct->options | PS_STAT_SAMPLE_LQ;
+        p_psVectorSampleQuartiles(myVector,
+                                  maskVector,
+                                  maskVal,
+                                  newStruct);
     }
 
@@ -495,8 +589,8 @@
     // Detemine minimum and maximum values in the data vector.
     if (isnan(newStruct->min)) {
-        p_psArrayMin(myVector, maskVector, maskVal, newStruct);
+        p_psVectorMin(myVector, maskVector, maskVal, newStruct);
     }
     if (isnan(newStruct->max)) {
-        p_psArrayMax(myVector, maskVector, maskVal, newStruct);
+        p_psVectorMax(myVector, maskVector, maskVal, newStruct);
     }
 
@@ -506,12 +600,14 @@
                                        binSize);
     // Populate the histogram arrat.
-    robustHistogram = psGetArrayHistogram(robustHistogram, myVector);
+    // GUS: fix this
+    //    robustHistogram = psHistogramVector(robustHistogram, myVector);
+    //
 
     // Smooth the histogram.
-    p_psArraysmoothHistGaussian(robustHistogram, sigmaE/4.0f);
+    p_psVectorsmoothHistGaussian(robustHistogram, sigmaE/4.0f);
 
     LQBinNum = -1;
     UQBinNum = -1;
-    for (i=0;i<robustHistogram->numBins;i++) {
+    for (i=0;i<robustHistogram->nums->n;i++) {
         if ((robustHistogram->nums->vec.i32[i] <= newStruct->sampleLQ) &&
                 (newStruct->sampleLQ <= robustHistogram->nums->vec.i32[i])) {
@@ -546,34 +642,28 @@
         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) {
+    if  (newStruct->options & PS_STAT_ROBUST_QUARTILE) {
         newStruct->robustUQ = 0.0;
-    }
-    if  (newStruct->options & PS_STAT_ROBUST_LQ) {
         newStruct->robustLQ = 0.0;
     }
 }
 
-void p_psArraySampleStdev(const psVector *restrict myVector,
-                          const psVector *restrict maskVector,
-                          unsigned int maskVal,
-                          psStats *newStruct)
+/*****************************************************************************
+    NOTE: This function assumes that p_psVectorMean() has already been called
+    and the correct value is stored in stats->sampleMean.
+ *****************************************************************************/
+void p_psVectorSampleStdev(const psVector *restrict myVector,
+                           const psVector *restrict maskVector,
+                           unsigned int maskVal,
+                           psStats *stats)
 {
     int i            = 0;
@@ -584,28 +674,55 @@
     float sumSquares = 0.0;
     float sumDiffs = 0.0;
-
-    if (0 != isnan(newStruct->sampleMean)) {
-        p_psArraySampleMean(myVector, maskVector, maskVal, newStruct);
-    }
-    mean = newStruct->sampleMean;
-
-    if (maskVector != NULL) {
-        for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.ui8[i])) {
+    float rangeMin = 0.0;
+    float rangeMax = 0.0;
+
+    if (0 != isnan(stats->sampleMean)) {
+        psAbort(__func__, "stats->sampleMean == NAN");
+    }
+    mean = stats->sampleMean;
+
+    if (stats->options & PS_STAT_USE_RANGE) {
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i]) &&
+                        (rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    diff = myVector->vec.f[i] - mean;
+                    sumSquares+= (diff * diff);
+                    sumDiffs+= diff;
+                    countInt++;
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                if ((rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    diff = myVector->vec.f[i] - mean;
+                    sumSquares+= (diff * diff);
+                    sumDiffs+= diff;
+                    countInt++;
+                }
+            }
+            countInt = myVector->n;
+        }
+    } else {
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i])) {
+                    diff = myVector->vec.f[i] - mean;
+                    sumSquares+= (diff * diff);
+                    sumDiffs+= diff;
+                    countInt++;
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
                 diff = myVector->vec.f[i] - mean;
                 sumSquares+= (diff * diff);
                 sumDiffs+= diff;
                 countInt++;
-
-            }
-        }
-    } else {
-        for (i=0;i<myVector->n;i++) {
-            diff = myVector->vec.f[i] - mean;
-            sumSquares+= (diff * diff);
-            sumDiffs+= diff;
-            countInt++;
-        }
-        countInt = myVector->n;
+            }
+            countInt = myVector->n;
+        }
     }
     countFloat = (float) countInt;
@@ -613,17 +730,17 @@
     #ifdef DARWIN
 
-    newStruct->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs *
-                                            sumDiffs/countFloat))/ (countFloat-1));
+    stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs *
+                                        sumDiffs/countFloat))/ (countFloat-1));
     #else
 
-    newStruct->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
-                                     sumDiffs/countFloat))/ (countFloat-1));
+    stats->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)
+void p_psVectorClippedStats(const psVector *restrict myVector,
+                            const psVector *restrict maskVector,
+                            unsigned int maskVal,
+                            psStats *newStruct)
 {
     int i = 0;
@@ -653,8 +770,8 @@
 
     // 1. Compute the sample median.
-    p_psArraySampleMedian(myVector, maskVector, maskVal, newStruct);
+    p_psVectorSampleMedian(myVector, maskVector, maskVal, newStruct);
 
     // 2. Compute the sample standard deviation.
-    p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
+    p_psVectorSampleStdev(myVector, maskVector, maskVal, newStruct);
 
     // 3. Use the sample median as the first estimator of the mean X.
@@ -676,6 +793,6 @@
             // 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);
+            p_psVectorSampleMedian(myVector, tmpMask, maskVal, newStruct);
+            p_psVectorSampleStdev(myVector, tmpMask, maskVal, newStruct);
 
             // c) Use the new mean for x
@@ -695,12 +812,4 @@
     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?
     }
 
@@ -716,85 +825,69 @@
  
  *****************************************************************************/
-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) {
+psStats *psVectorStats(psStats *stats,
+                       psVector *in,
+                       psVector *mask,
+                       unsigned int maskVal)
+{
+    if (in == NULL) {
+        psAbort(__func__, "Input data vector (in) was NULL.");
+    }
+
+    if (in->type.type != PS_TYPE_FLOAT) {
         psAbort(__func__,
                 "Only data type PS_TYPE_FLOAT is currently supported.");
     }
-    if (maskVector != NULL) {
-        if (myVector->n != maskVector->n) {
+    if (mask != NULL) {
+        if (in->n != mask->n) {
             psAbort(__func__,
                     "Vector data and vector mask are of different sizes.");
         }
-        if (maskVector->type.type != PS_TYPE_UINT8) {
+        if (mask->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);
+        p_psVectorSampleMean(in, mask, maskVal, stats);
+    }
+
+    // ************************************************************************
+    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
+        p_psVectorSampleMedian(in, mask, maskVal, stats);
+    }
+
+    // ************************************************************************
+    if (stats->options & PS_STAT_SAMPLE_STDEV) {
+        p_psVectorSampleStdev(in, mask, maskVal, stats);
+    }
+
+    // ************************************************************************
+    if (stats->options & PS_STAT_SAMPLE_QUARTILE) {
+        p_psVectorSampleQuartiles(in, mask, maskVal, stats);
+    }
+
+    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
+            (stats->options & PS_STAT_ROBUST_MEDIAN) ||
+            (stats->options & PS_STAT_ROBUST_MODE) ||
+            (stats->options & PS_STAT_ROBUST_STDEV) ||
+            (stats->options & PS_STAT_ROBUST_QUARTILE)) {
+        p_psVectorRobustStats(in, mask, maskVal, stats);
+    }
+
+    if ((stats->options & PS_STAT_CLIPPED_MEAN) ||
+            (stats->options & PS_STAT_CLIPPED_STDEV)) {
+        p_psVectorClippedStats(in, mask, maskVal, stats);
     }
 
     // ************************************************************************
     if (stats->options & PS_STAT_MAX) {
-        p_psArrayMax(myVector, maskVector, maskVal, newStruct);
+        p_psVectorMax(in, mask, maskVal, stats);
     }
 
     // ************************************************************************
     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);
+        p_psVectorMin(in, mask, maskVal, stats);
     }
 
@@ -803,4 +896,4 @@
     //        psAbort(__func__, "Unknown options 0x%x.\n", stats->options);
 
-    return(newStruct);
-}
+    return(stats);
+}