Index: /trunk/psLib/src/dataManip/psStats.c
===================================================================
--- /trunk/psLib/src/dataManip/psStats.c	(revision 795)
+++ /trunk/psLib/src/dataManip/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);
+}
Index: /trunk/psLib/src/dataManip/psStats.h
===================================================================
--- /trunk/psLib/src/dataManip/psStats.h	(revision 795)
+++ /trunk/psLib/src/dataManip/psStats.h	(revision 796)
@@ -1,5 +1,7 @@
 #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
@@ -7,28 +9,24 @@
  */
 #include "psVector.h"
-
-/** statistics which may be calculated */
+/******************************************************************************
+    Statistical functions and data structures.
+ *****************************************************************************/
 typedef enum {
     PS_STAT_SAMPLE_MEAN           = 0x000001,
     PS_STAT_SAMPLE_MEDIAN         = 0x000002,
     PS_STAT_SAMPLE_STDEV          = 0x000004,
-    PS_STAT_SAMPLE_UQ             = 0x000008,
-    PS_STAT_SAMPLE_LQ             = 0x000010,
-    PS_STAT_ROBUST_MEAN           = 0x000020,
-    PS_STAT_ROBUST_MEAN_NVALUES   = 0x000040,
-    PS_STAT_ROBUST_MEDIAN         = 0x000080,
-    PS_STAT_ROBUST_MEDIAN_NVALUES = 0x000100,
-    PS_STAT_ROBUST_MODE           = 0x000200,
-    PS_STAT_ROBUST_MODE_NVALUES   = 0x000400,
-    PS_STAT_ROBUST_STDEV          = 0x000800,
-    PS_STAT_ROBUST_UQ             = 0x001000,
-    PS_STAT_ROBUST_LQ             = 0x002000,
-    PS_STAT_CLIPPED_MEAN          = 0x004000,
-    PS_STAT_CLIPPED_MEAN_NVALUES  = 0x008000,
-    PS_STAT_CLIPPED_MEAN_NSIGMA   = 0x010000,
-    PS_STAT_CLIPPED_STDEV         = 0x020000,
-    PS_STAT_MAX                   = 0x040000,
-    PS_STAT_MIN                   = 0x080000,
-    PS_STAT_NVALUES               = 0x100000
+    PS_STAT_SAMPLE_QUARTILE       = 0x000008,
+    PS_STAT_ROBUST_MEAN           = 0x000010,
+    PS_STAT_ROBUST_MEDIAN         = 0x000020,
+    PS_STAT_ROBUST_MODE           = 0x000040,
+    PS_STAT_ROBUST_STDEV          = 0x000080,
+    PS_STAT_ROBUST_QUARTILE       = 0x000100,
+    PS_STAT_CLIPPED_MEAN          = 0x000200,
+    PS_STAT_CLIPPED_STDEV         = 0x000400,
+    PS_STAT_MAX                   = 0x000800,
+    PS_STAT_MIN                   = 0x001000,
+    PS_STAT_USE_RANGE             = 0x002000,
+    PS_STAT_USE_BINSIZE           = 0x004000,
+    PS_STAT_ROBUST_FOR_SAMPLE     = 0x008000
 } psStatsOptions;
 
@@ -42,15 +40,15 @@
     double sampleUQ;            ///< upper quartile of sample
     double sampleLQ;            ///< lower quartile of sample
+    double sampleLimit;         ///<
     double robustMean;          ///< robust mean of array
     int    robustMeanNvalues;   ///< number of measurements used for robust mean
     double robustMedian;        ///< robust median of array
-    int    robustMedianNvalues; ///< number of measurements used for robust median
     double robustMode;          ///< Robust mode of array
-    int    robustModeNvalues;   ///< Number of measurements used for robust mode
     double robustStdev;         ///< robust standard deviation of array
     double robustUQ;            ///< robust upper quartile
     double robustLQ;            ///< robust lower quartile
+    double robustN50;           ///<
+    double robustNfit;          ///<
     double clippedMean;         ///< Nsigma clipped mean
-    int    clippedMeanNvalues;  ///< number of data points used for clipped mean
     double clippedStdev;        ///< standard deviation after clipping
     double clipSigma;           ///< Nsigma used for clipping; user input
@@ -58,5 +56,5 @@
     double min;                 ///< minimum data value in array
     double max;                 ///< maximum data value in array
-    int    nValues;             ///< number of data values in array
+    double binsize;             ///<
     psStatsOptions options;     ///< bitmask of calculated values
 }
@@ -66,34 +64,29 @@
 /** Do Statistics on an array.  Returns a status value. \ingroup MathGroup */
 psStats *
-psArrayStats(const psVector *restrict myVector, ///< Vector to be analysed: must be FLOAT
-             const psVector *restrict maskVector, ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
-             unsigned int maskVal, ///< Only mask elements with one of these bits set in maskVector
-             psStats *stats  ///< stats structure defines stats to be calculated and how
-            );
+psVectorStats(psStats *stats,       ///< stats structure defines stats to be calculated and how
+              psVector *in,         ///< Vector to be analysed: must be F32
+              psVector *mask,       ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
+              unsigned int maskVal ///< Only mask elements with one of these bits set in maskVector
+             );
 
 /** Constructor */
-psStats *
-psStatsAlloc(psStatsOptions options ///< Statistics to measure
-            );
+psStats *psStatsAlloc(psStatsOptions options); ///< Statistics to measure
 
 /** Destructor */
-void
-psStatsFree(psStats *restrict stats ///< Stats structure to destroy
-           );
+void psStatsFree(psStats *restrict stats); ///< Stats structure to destroy
 
-/** Histograms  */
-// NOTE: I added the numBins member.  Must consult with IfA on this.
+/******************************************************************************
+    Histogram functions and data structures.
+ *****************************************************************************/
 typedef struct
 {
-    const psVector *restrict lower; ///< Lower bounds for the bins (FLOAT)
-    const psVector *restrict upper; ///< Upper bounds for the bins (FLOAT)
-    psVector *nums;                 ///< Number in each of the bins (INT)
-    float minVal;
-    float maxVal;                   ///< Minimum and maximum values
-    int minNum;
-    int maxNum;                     ///< Number below the minimum and above the maximum
-    int numBins;                    ///< Number of bins in this histogram.
+    const psVector *restrict bounds; ///< Bounds for the bins (type F32)
+    psVector *nums;                  ///< Number in each of the bins (INT)
+    int minNum;                      ///< Number below the minimum
+    int maxNum;                      ///< Number above the maximum
+    int uniform;                     ///< Is it a uniform distribution?
 }
 psHistogram;
+
 
 /** Constructor \ingroup MathGroup */
@@ -101,26 +94,21 @@
 psHistogramAlloc(float lower,  ///< Lower limit for the bins
                  float upper,  ///< Upper limit for the bins
-                 float size    ///< Size of the bins
-                );
+                 int n);       ///< Number of bins
+
 
 /** Generic constructor \ingroup MathGroup */
-psHistogram *
-psHistogramAllocGeneric(const psVector *restrict lower, ///< Lower bounds for the bins
-                        const psVector *restrict upper, ///< Upper bounds for the bins
-                        float minVal,                   ///< Minimum value
-                        float maxVal                    ///< Maximum value
-                       );
+psHistogram * psHistogramAllocGeneric(const psVector *restrict bounds); ///< Bounds for the bins
+
 
 /** Destructor \ingroup MathGroup **/
-void
-psHistogramFree(psHistogram *restrict myHist          ///< Histogram to destroy
-               );
+void psHistogramFree(psHistogram *restrict myHist);          ///< Histogram to destroy
 
 
 /** Calculate a histogram \ingroup MathGroup **/
-psHistogram *
-psGetArrayHistogram(psHistogram *restrict myHist,     ///< Histogram data
-                    const psVector *restrict myVector ///< Vector to analyse
-                   );
+psHistogram *psHistogramVector (psHistogram *out,   ///< Histogram data
+                                psVector *in,       ///< Vector to analyse
+                                psVector *mask,     ///< Mask dat for input vector
+                                int maskVal);       ///< Mask value
+
 
 #endif
Index: /trunk/psLib/src/image/psImageStats.c
===================================================================
--- /trunk/psLib/src/image/psImageStats.c	(revision 795)
+++ /trunk/psLib/src/image/psImageStats.c	(revision 796)
@@ -28,5 +28,5 @@
 
     junk = psVectorAlloc(input->numRows * input->numCols, input->type.type);
-    stats2 = psArrayStats(junk, NULL, 0xffffffff, stats);
+    stats2 = psVectorStats(stats, junk, NULL, 0xffffffff);
     psVectorFree(junk);
 
Index: /trunk/psLib/src/imageops/psImageStats.c
===================================================================
--- /trunk/psLib/src/imageops/psImageStats.c	(revision 795)
+++ /trunk/psLib/src/imageops/psImageStats.c	(revision 796)
@@ -28,5 +28,5 @@
 
     junk = psVectorAlloc(input->numRows * input->numCols, input->type.type);
-    stats2 = psArrayStats(junk, NULL, 0xffffffff, stats);
+    stats2 = psVectorStats(stats, junk, NULL, 0xffffffff);
     psVectorFree(junk);
 
Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 795)
+++ /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);
+}
Index: /trunk/psLib/src/math/psStats.h
===================================================================
--- /trunk/psLib/src/math/psStats.h	(revision 795)
+++ /trunk/psLib/src/math/psStats.h	(revision 796)
@@ -1,5 +1,7 @@
 #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
@@ -7,28 +9,24 @@
  */
 #include "psVector.h"
-
-/** statistics which may be calculated */
+/******************************************************************************
+    Statistical functions and data structures.
+ *****************************************************************************/
 typedef enum {
     PS_STAT_SAMPLE_MEAN           = 0x000001,
     PS_STAT_SAMPLE_MEDIAN         = 0x000002,
     PS_STAT_SAMPLE_STDEV          = 0x000004,
-    PS_STAT_SAMPLE_UQ             = 0x000008,
-    PS_STAT_SAMPLE_LQ             = 0x000010,
-    PS_STAT_ROBUST_MEAN           = 0x000020,
-    PS_STAT_ROBUST_MEAN_NVALUES   = 0x000040,
-    PS_STAT_ROBUST_MEDIAN         = 0x000080,
-    PS_STAT_ROBUST_MEDIAN_NVALUES = 0x000100,
-    PS_STAT_ROBUST_MODE           = 0x000200,
-    PS_STAT_ROBUST_MODE_NVALUES   = 0x000400,
-    PS_STAT_ROBUST_STDEV          = 0x000800,
-    PS_STAT_ROBUST_UQ             = 0x001000,
-    PS_STAT_ROBUST_LQ             = 0x002000,
-    PS_STAT_CLIPPED_MEAN          = 0x004000,
-    PS_STAT_CLIPPED_MEAN_NVALUES  = 0x008000,
-    PS_STAT_CLIPPED_MEAN_NSIGMA   = 0x010000,
-    PS_STAT_CLIPPED_STDEV         = 0x020000,
-    PS_STAT_MAX                   = 0x040000,
-    PS_STAT_MIN                   = 0x080000,
-    PS_STAT_NVALUES               = 0x100000
+    PS_STAT_SAMPLE_QUARTILE       = 0x000008,
+    PS_STAT_ROBUST_MEAN           = 0x000010,
+    PS_STAT_ROBUST_MEDIAN         = 0x000020,
+    PS_STAT_ROBUST_MODE           = 0x000040,
+    PS_STAT_ROBUST_STDEV          = 0x000080,
+    PS_STAT_ROBUST_QUARTILE       = 0x000100,
+    PS_STAT_CLIPPED_MEAN          = 0x000200,
+    PS_STAT_CLIPPED_STDEV         = 0x000400,
+    PS_STAT_MAX                   = 0x000800,
+    PS_STAT_MIN                   = 0x001000,
+    PS_STAT_USE_RANGE             = 0x002000,
+    PS_STAT_USE_BINSIZE           = 0x004000,
+    PS_STAT_ROBUST_FOR_SAMPLE     = 0x008000
 } psStatsOptions;
 
@@ -42,15 +40,15 @@
     double sampleUQ;            ///< upper quartile of sample
     double sampleLQ;            ///< lower quartile of sample
+    double sampleLimit;         ///<
     double robustMean;          ///< robust mean of array
     int    robustMeanNvalues;   ///< number of measurements used for robust mean
     double robustMedian;        ///< robust median of array
-    int    robustMedianNvalues; ///< number of measurements used for robust median
     double robustMode;          ///< Robust mode of array
-    int    robustModeNvalues;   ///< Number of measurements used for robust mode
     double robustStdev;         ///< robust standard deviation of array
     double robustUQ;            ///< robust upper quartile
     double robustLQ;            ///< robust lower quartile
+    double robustN50;           ///<
+    double robustNfit;          ///<
     double clippedMean;         ///< Nsigma clipped mean
-    int    clippedMeanNvalues;  ///< number of data points used for clipped mean
     double clippedStdev;        ///< standard deviation after clipping
     double clipSigma;           ///< Nsigma used for clipping; user input
@@ -58,5 +56,5 @@
     double min;                 ///< minimum data value in array
     double max;                 ///< maximum data value in array
-    int    nValues;             ///< number of data values in array
+    double binsize;             ///<
     psStatsOptions options;     ///< bitmask of calculated values
 }
@@ -66,34 +64,29 @@
 /** Do Statistics on an array.  Returns a status value. \ingroup MathGroup */
 psStats *
-psArrayStats(const psVector *restrict myVector, ///< Vector to be analysed: must be FLOAT
-             const psVector *restrict maskVector, ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
-             unsigned int maskVal, ///< Only mask elements with one of these bits set in maskVector
-             psStats *stats  ///< stats structure defines stats to be calculated and how
-            );
+psVectorStats(psStats *stats,       ///< stats structure defines stats to be calculated and how
+              psVector *in,         ///< Vector to be analysed: must be F32
+              psVector *mask,       ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
+              unsigned int maskVal ///< Only mask elements with one of these bits set in maskVector
+             );
 
 /** Constructor */
-psStats *
-psStatsAlloc(psStatsOptions options ///< Statistics to measure
-            );
+psStats *psStatsAlloc(psStatsOptions options); ///< Statistics to measure
 
 /** Destructor */
-void
-psStatsFree(psStats *restrict stats ///< Stats structure to destroy
-           );
+void psStatsFree(psStats *restrict stats); ///< Stats structure to destroy
 
-/** Histograms  */
-// NOTE: I added the numBins member.  Must consult with IfA on this.
+/******************************************************************************
+    Histogram functions and data structures.
+ *****************************************************************************/
 typedef struct
 {
-    const psVector *restrict lower; ///< Lower bounds for the bins (FLOAT)
-    const psVector *restrict upper; ///< Upper bounds for the bins (FLOAT)
-    psVector *nums;                 ///< Number in each of the bins (INT)
-    float minVal;
-    float maxVal;                   ///< Minimum and maximum values
-    int minNum;
-    int maxNum;                     ///< Number below the minimum and above the maximum
-    int numBins;                    ///< Number of bins in this histogram.
+    const psVector *restrict bounds; ///< Bounds for the bins (type F32)
+    psVector *nums;                  ///< Number in each of the bins (INT)
+    int minNum;                      ///< Number below the minimum
+    int maxNum;                      ///< Number above the maximum
+    int uniform;                     ///< Is it a uniform distribution?
 }
 psHistogram;
+
 
 /** Constructor \ingroup MathGroup */
@@ -101,26 +94,21 @@
 psHistogramAlloc(float lower,  ///< Lower limit for the bins
                  float upper,  ///< Upper limit for the bins
-                 float size    ///< Size of the bins
-                );
+                 int n);       ///< Number of bins
+
 
 /** Generic constructor \ingroup MathGroup */
-psHistogram *
-psHistogramAllocGeneric(const psVector *restrict lower, ///< Lower bounds for the bins
-                        const psVector *restrict upper, ///< Upper bounds for the bins
-                        float minVal,                   ///< Minimum value
-                        float maxVal                    ///< Maximum value
-                       );
+psHistogram * psHistogramAllocGeneric(const psVector *restrict bounds); ///< Bounds for the bins
+
 
 /** Destructor \ingroup MathGroup **/
-void
-psHistogramFree(psHistogram *restrict myHist          ///< Histogram to destroy
-               );
+void psHistogramFree(psHistogram *restrict myHist);          ///< Histogram to destroy
 
 
 /** Calculate a histogram \ingroup MathGroup **/
-psHistogram *
-psGetArrayHistogram(psHistogram *restrict myHist,     ///< Histogram data
-                    const psVector *restrict myVector ///< Vector to analyse
-                   );
+psHistogram *psHistogramVector (psHistogram *out,   ///< Histogram data
+                                psVector *in,       ///< Vector to analyse
+                                psVector *mask,     ///< Mask dat for input vector
+                                int maskVal);       ///< Mask value
+
 
 #endif