Index: /trunk/psLib/src/dataManip/psStats.c
===================================================================
--- /trunk/psLib/src/dataManip/psStats.c	(revision 698)
+++ /trunk/psLib/src/dataManip/psStats.c	(revision 699)
@@ -270,44 +270,30 @@
     psVector *unsortedVector = NULL;
     psVector *sortedVector = NULL;
-    int dataSize = 0;
     int count = 0;
     int i = 0;
     float median = 0.0;
 
-    p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
-
-    if (dataSize < MEDIAN_SIZE_THRESHOLD) {
-        unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, dataSize);
-        sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, dataSize);
-
-        count = 0;
-        if (maskVector != NULL) {
-            for (i=0;i<myVector->n;i++) {
-                if (!(maskVal & maskVector->vec.ui8[i])) {
-                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
-                }
-            }
-            psSort(sortedVector, unsortedVector);
-        } else {
-            psSort(sortedVector, myVector);
-        }
-
-        if (0 == (dataSize % 2)) {
-            median = 0.5 * (sortedVector->vec.f[(dataSize/2)-1] +
-                            sortedVector->vec.f[dataSize/2]);
-        } else {
-            median = sortedVector->vec.f[dataSize/2];
-        }
-    } else {
-        // BROAD: Calculate the Robust Median
-        // Determine the LQ of the distribution.
-        // Determine the UQ of the distribution.
-        // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
-        // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
-        // Find the bin with the peak value between LQ and UQ (the MODE)
-        // dL = (UQ - LQ) / 8
-        // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
-        // The resulting fit parameters are the robust mean, mean_r, and sigma
-        psError(__func__, "GUS: p_psArraySampleMedian(), large data sample");
+    if (-1 == newStruct->nValues) {
+        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
+    }
+    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+
+    if (maskVector != NULL) {
+        for (i=0;i<myVector->n;i++) {
+            if (!(maskVal & maskVector->vec.ui8[i])) {
+                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+            }
+        }
+        psSort(sortedVector, unsortedVector);
+    } else {
+        psSort(sortedVector, myVector);
+    }
+
+    if (0 == (newStruct->nValues % 2)) {
+        median = 0.5 * (sortedVector->vec.f[(newStruct->nValues/2)-1] +
+                        sortedVector->vec.f[newStruct->nValues/2]);
+    } else {
+        median = sortedVector->vec.f[newStruct->nValues/2];
     }
 
@@ -315,4 +301,39 @@
     psVectorFree(sortedVector);
     newStruct->sampleMedian = median;
+}
+
+void p_psArrayRobustMedian(const psVector *restrict myVector,
+                           const psVector *restrict maskVector,
+                           unsigned int maskVal,
+                           psStats *newStruct)
+{
+    psHistogram *robustHistogram = NULL;
+    float binSize = 0.0;
+
+    //    if (isnan(myVector->robustLQ) ||
+    //        isnan(myVector->robustUQ)) {
+    //        p_psArrayRobustQuartiles(myVector, maskVector, maskVal, newStruct);
+    //    }
+    //    binSize = ((myVector->robustUQ - myVector->robustLQ) / 1.34) / 10.0;
+
+    robustHistogram = psHistogramAlloc(newStruct->min,
+                                       newStruct->max,
+                                       binSize);
+    //    p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
+    //    robustHistogram = psGetArrayHistogram(robustHistogram, myVector);
+    //    p_psArraySmooth(robustHistogram, (binSize / 4.0));
+    //    dL = (myVector->robustUQ - myVector->robustLQ) / 8.0;
+
+
+    // BROAD: Calculate the Robust Median
+    // Determine the LQ of the distribution.
+    // Determine the UQ of the distribution.
+    // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
+    // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
+
+    // Find the bin with the peak value between LQ and UQ (the MODE)
+    // dL = (UQ - LQ) / 8
+    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
+    // The resulting fit parameters are the robust mean, mean_r, and sigma
 }
 
@@ -368,10 +389,22 @@
 }
 
-void p_psArraySampleUQ(const psVector *restrict myVector,
-                       const psVector *restrict maskVector,
-                       unsigned int maskVal,
-                       psStats *newStruct)
-{
-    psHistogram    *nonRobustHistogram = NULL;
+/******************************************************************************
+ *****************************************************************************/
+void p_psArraySampleQuartiles(const psVector *restrict myVector,
+                              const psVector *restrict maskVector,
+                              unsigned int maskVal,
+                              psStats *newStruct)
+{
+    psVector *unsortedVector = NULL;
+    psVector *sortedVector = NULL;
+    int count = 0;
+    int ind = 0;
+    int i = 0;
+
+    // return is we have already calculated both quartile points.
+    if ((!isnan(newStruct->sampleLQ)) &&
+            (!isnan(newStruct->sampleUQ))) {
+        return;
+    }
 
     if (-1 == newStruct->nValues) {
@@ -379,24 +412,28 @@
     }
 
-    if (newStruct->nValues < MEDIAN_SIZE_THRESHOLD) {
-
-    }
-    else {
-        if (0 != isnan(newStruct->sampleStdev)) {
-            p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
-        }
-        if (0 != isnan(newStruct->min)) {
-            p_psArrayMin(myVector, maskVector, maskVal, newStruct);
-        }
-        if (0 != isnan(newStruct->max)) {
-            p_psArrayMax(myVector, maskVector, maskVal, newStruct);
-        }
-
-        nonRobustHistogram = psHistogramAlloc(newStruct->min,
-                                              newStruct->max,
-                                              10);
-    }
-}
-
+
+    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+
+    count = 0;
+    if (maskVector != NULL) {
+        for (i=0;i<myVector->n;i++) {
+            if (!(maskVal & maskVector->vec.ui8[i])) {
+                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+            }
+        }
+        psSort(sortedVector, unsortedVector);
+    } else {
+        psSort(sortedVector, myVector);
+    }
+
+    ind = 3 * (newStruct->nValues / 4);
+    newStruct->sampleUQ = sortedVector->vec.f[ind];
+    ind = (newStruct->nValues / 4);
+    newStruct->sampleLQ = sortedVector->vec.f[ind];
+
+    psVectorFree(unsortedVector);
+    psVectorFree(sortedVector);
+}
 
 /******************************************************************************
@@ -436,33 +473,42 @@
     newStruct = psStatsAlloc(stats->options);
 
+    // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_MEAN) {
         p_psArraySampleMean(myVector, maskVector, maskVal, newStruct);
     }
 
+    // ************************************************************************
     if (stats->options & PS_STAT_MAX) {
         p_psArrayMax(myVector, maskVector, maskVal, newStruct);
     }
 
+    // ************************************************************************
     if (stats->options & PS_STAT_MIN) {
         p_psArrayMin(myVector, maskVector, maskVal, newStruct);
     }
 
+    // ************************************************************************
     if (stats->options & PS_STAT_NVALUES) {
         p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
     }
 
-    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {}
-
+    // ************************************************************************
+    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) {
-        p_psArraySampleUQ(myVector, maskVector, maskVal, newStruct);
-    }
-
-    if (stats->options & PS_STAT_SAMPLE_LQ) {
-        newStruct->sampleLQ = p_psArrayXXX(myVector, maskVector, maskVal);
-    }
+    // ************************************************************************
+    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) {
@@ -494,9 +540,6 @@
     }
 
-    if (stats->options & PS_STAT_ROBUST_UQ) {
-        newStruct->robustUQ = p_psArrayXXX(myVector, maskVector, maskVal);
-    }
-
-    if (stats->options & PS_STAT_ROBUST_LQ) {
+    if ((stats->options & PS_STAT_ROBUST_UQ) ||
+            (stats->options & PS_STAT_ROBUST_LQ)) {
         newStruct->robustLQ = p_psArrayXXX(myVector, maskVector, maskVal);
     }
Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 698)
+++ /trunk/psLib/src/math/psStats.c	(revision 699)
@@ -270,44 +270,30 @@
     psVector *unsortedVector = NULL;
     psVector *sortedVector = NULL;
-    int dataSize = 0;
     int count = 0;
     int i = 0;
     float median = 0.0;
 
-    p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
-
-    if (dataSize < MEDIAN_SIZE_THRESHOLD) {
-        unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, dataSize);
-        sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, dataSize);
-
-        count = 0;
-        if (maskVector != NULL) {
-            for (i=0;i<myVector->n;i++) {
-                if (!(maskVal & maskVector->vec.ui8[i])) {
-                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
-                }
-            }
-            psSort(sortedVector, unsortedVector);
-        } else {
-            psSort(sortedVector, myVector);
-        }
-
-        if (0 == (dataSize % 2)) {
-            median = 0.5 * (sortedVector->vec.f[(dataSize/2)-1] +
-                            sortedVector->vec.f[dataSize/2]);
-        } else {
-            median = sortedVector->vec.f[dataSize/2];
-        }
-    } else {
-        // BROAD: Calculate the Robust Median
-        // Determine the LQ of the distribution.
-        // Determine the UQ of the distribution.
-        // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
-        // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
-        // Find the bin with the peak value between LQ and UQ (the MODE)
-        // dL = (UQ - LQ) / 8
-        // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
-        // The resulting fit parameters are the robust mean, mean_r, and sigma
-        psError(__func__, "GUS: p_psArraySampleMedian(), large data sample");
+    if (-1 == newStruct->nValues) {
+        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
+    }
+    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+
+    if (maskVector != NULL) {
+        for (i=0;i<myVector->n;i++) {
+            if (!(maskVal & maskVector->vec.ui8[i])) {
+                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+            }
+        }
+        psSort(sortedVector, unsortedVector);
+    } else {
+        psSort(sortedVector, myVector);
+    }
+
+    if (0 == (newStruct->nValues % 2)) {
+        median = 0.5 * (sortedVector->vec.f[(newStruct->nValues/2)-1] +
+                        sortedVector->vec.f[newStruct->nValues/2]);
+    } else {
+        median = sortedVector->vec.f[newStruct->nValues/2];
     }
 
@@ -315,4 +301,39 @@
     psVectorFree(sortedVector);
     newStruct->sampleMedian = median;
+}
+
+void p_psArrayRobustMedian(const psVector *restrict myVector,
+                           const psVector *restrict maskVector,
+                           unsigned int maskVal,
+                           psStats *newStruct)
+{
+    psHistogram *robustHistogram = NULL;
+    float binSize = 0.0;
+
+    //    if (isnan(myVector->robustLQ) ||
+    //        isnan(myVector->robustUQ)) {
+    //        p_psArrayRobustQuartiles(myVector, maskVector, maskVal, newStruct);
+    //    }
+    //    binSize = ((myVector->robustUQ - myVector->robustLQ) / 1.34) / 10.0;
+
+    robustHistogram = psHistogramAlloc(newStruct->min,
+                                       newStruct->max,
+                                       binSize);
+    //    p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
+    //    robustHistogram = psGetArrayHistogram(robustHistogram, myVector);
+    //    p_psArraySmooth(robustHistogram, (binSize / 4.0));
+    //    dL = (myVector->robustUQ - myVector->robustLQ) / 8.0;
+
+
+    // BROAD: Calculate the Robust Median
+    // Determine the LQ of the distribution.
+    // Determine the UQ of the distribution.
+    // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
+    // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
+
+    // Find the bin with the peak value between LQ and UQ (the MODE)
+    // dL = (UQ - LQ) / 8
+    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
+    // The resulting fit parameters are the robust mean, mean_r, and sigma
 }
 
@@ -368,10 +389,22 @@
 }
 
-void p_psArraySampleUQ(const psVector *restrict myVector,
-                       const psVector *restrict maskVector,
-                       unsigned int maskVal,
-                       psStats *newStruct)
-{
-    psHistogram    *nonRobustHistogram = NULL;
+/******************************************************************************
+ *****************************************************************************/
+void p_psArraySampleQuartiles(const psVector *restrict myVector,
+                              const psVector *restrict maskVector,
+                              unsigned int maskVal,
+                              psStats *newStruct)
+{
+    psVector *unsortedVector = NULL;
+    psVector *sortedVector = NULL;
+    int count = 0;
+    int ind = 0;
+    int i = 0;
+
+    // return is we have already calculated both quartile points.
+    if ((!isnan(newStruct->sampleLQ)) &&
+            (!isnan(newStruct->sampleUQ))) {
+        return;
+    }
 
     if (-1 == newStruct->nValues) {
@@ -379,24 +412,28 @@
     }
 
-    if (newStruct->nValues < MEDIAN_SIZE_THRESHOLD) {
-
-    }
-    else {
-        if (0 != isnan(newStruct->sampleStdev)) {
-            p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
-        }
-        if (0 != isnan(newStruct->min)) {
-            p_psArrayMin(myVector, maskVector, maskVal, newStruct);
-        }
-        if (0 != isnan(newStruct->max)) {
-            p_psArrayMax(myVector, maskVector, maskVal, newStruct);
-        }
-
-        nonRobustHistogram = psHistogramAlloc(newStruct->min,
-                                              newStruct->max,
-                                              10);
-    }
-}
-
+
+    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
+
+    count = 0;
+    if (maskVector != NULL) {
+        for (i=0;i<myVector->n;i++) {
+            if (!(maskVal & maskVector->vec.ui8[i])) {
+                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+            }
+        }
+        psSort(sortedVector, unsortedVector);
+    } else {
+        psSort(sortedVector, myVector);
+    }
+
+    ind = 3 * (newStruct->nValues / 4);
+    newStruct->sampleUQ = sortedVector->vec.f[ind];
+    ind = (newStruct->nValues / 4);
+    newStruct->sampleLQ = sortedVector->vec.f[ind];
+
+    psVectorFree(unsortedVector);
+    psVectorFree(sortedVector);
+}
 
 /******************************************************************************
@@ -436,33 +473,42 @@
     newStruct = psStatsAlloc(stats->options);
 
+    // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_MEAN) {
         p_psArraySampleMean(myVector, maskVector, maskVal, newStruct);
     }
 
+    // ************************************************************************
     if (stats->options & PS_STAT_MAX) {
         p_psArrayMax(myVector, maskVector, maskVal, newStruct);
     }
 
+    // ************************************************************************
     if (stats->options & PS_STAT_MIN) {
         p_psArrayMin(myVector, maskVector, maskVal, newStruct);
     }
 
+    // ************************************************************************
     if (stats->options & PS_STAT_NVALUES) {
         p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
     }
 
-    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {}
-
+    // ************************************************************************
+    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) {
-        p_psArraySampleUQ(myVector, maskVector, maskVal, newStruct);
-    }
-
-    if (stats->options & PS_STAT_SAMPLE_LQ) {
-        newStruct->sampleLQ = p_psArrayXXX(myVector, maskVector, maskVal);
-    }
+    // ************************************************************************
+    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) {
@@ -494,9 +540,6 @@
     }
 
-    if (stats->options & PS_STAT_ROBUST_UQ) {
-        newStruct->robustUQ = p_psArrayXXX(myVector, maskVector, maskVal);
-    }
-
-    if (stats->options & PS_STAT_ROBUST_LQ) {
+    if ((stats->options & PS_STAT_ROBUST_UQ) ||
+            (stats->options & PS_STAT_ROBUST_LQ)) {
         newStruct->robustLQ = p_psArrayXXX(myVector, maskVector, maskVal);
     }