Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 6322)
+++ trunk/psLib/src/math/psStats.c	(revision 6437)
@@ -16,6 +16,6 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.167 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-03 22:05:22 $
+ *  @version $Revision: 1.168 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 00:56:48 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -440,11 +440,10 @@
 unmasked element within the specified min/max range).  Otherwise, return
 "false".
- 
-XXX: Can you use psVectorCountPixelMask here?
  *****************************************************************************/
-bool p_psVectorCheckNonEmpty(const psVector* myVector,
-                             const psVector* maskVector,
-                             psU32 maskVal,
-                             psStats* stats)
+bool p_psVectorCheckNonEmpty(
+    const psVector* myVector,
+    const psVector* maskVector,
+    psU32 maskVal,
+    psStats* stats)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -496,6 +495,4 @@
 number of non-masked pixels in the vector that fall within the min/max
 range, if specified.
- 
-XXX: Can you use psVectorCountPixelMask here?
  *****************************************************************************/
 psS32 p_psVectorNValues(const psVector* myVector,
@@ -644,9 +641,4 @@
 robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
 smoothed data.
- 
-XXX: Only PS_TYPE_F32 is supported.
- 
-XXX: Write a general routine which smoothes a psVector.  This routine should
-call that.  Is that possible?
  *****************************************************************************/
 psVector *p_psVectorSmoothHistGaussian(
@@ -985,6 +977,4 @@
     -1: error
     -2: warning
- 
-XXX: Use static vectors for tmpMask.
  *****************************************************************************/
 psS32 p_psVectorClippedStats(const psVector* myVector,
@@ -1224,8 +1214,5 @@
 
         //
-        // Ensure that the y values are monotonic.
-        //
-        // XXX: This routine should probably be rewritten in a more general fashion
-        // so that the following checks are not necessary.
+        // Ensure that the y value lies within range of the y values.
         //
         if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
@@ -1236,4 +1223,7 @@
         }
 
+        //
+        // Ensure that the y values are monotonic.
+        //
         if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
                 ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
@@ -1289,6 +1279,5 @@
         } else if (binNum == (xVec->n - 2)) {
             // XXX: Is this right?
-            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
-                              xVec->data.F32[binNum + 1]);
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]);
         }
     }
@@ -1677,5 +1666,5 @@
     tmpScalar.data.F32 = totalDataPoints * 0.25f;
     if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
-        // XXX: Special case where its in first bin.  Must code last bin.
+        // Special case where its in first bin.  The last bin should be handled automatically.
         binLo25 = 0;
     } else {
@@ -1684,5 +1673,5 @@
     tmpScalar.data.F32 = totalDataPoints * 0.75f;
     if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
-        // XXX: Special case where its in first bin.  Must code last bin.
+        // Special case where its in first bin.  The last bin should be handled automatically.
         binHi25 = 0;
     } else {
@@ -1704,5 +1693,4 @@
     // Interpolate to find these two positions exactly: these are the upper
     // and lower quartile positions.
-    // XXX: Check for errors.
     //
     psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(
@@ -1711,4 +1699,5 @@
                            binLo25,
                            totalDataPoints * 0.25f);
+
     psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(
                            *(psVector* *)&cumulativeRobustHistogram->bounds,
@@ -1716,8 +1705,18 @@
                            binHi25,
                            totalDataPoints * 0.75f);
+
+    if (isnan(binLo25F32) || isnan(binHi25F32)) {
+        psError(PS_ERR_UNKNOWN, false, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n");
+        psFree(tmpStatsMinMax);
+        psFree(robustHistogram);
+        psFree(cumulativeRobustHistogram);
+        psFree(tmpMaskVec);
+        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
+        return(1);
+    }
+
     stats->robustLQ = binLo25F32;
     stats->robustUQ = binHi25F32;
     psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32);
-
     psS32 N50 = 0;
     for (psS32 i = 0 ; i < myVector->n ; i++) {
@@ -1863,15 +1862,16 @@
         }
 
-        // XXX: Use the min/max routines for this
-        psF32 minY = FLT_MAX;
-        psF32 maxY = -FLT_MAX;
-        for (psS32 i = 0 ; i < y->n ; i++) {
-            if (y->data.F32[i] > maxY) {
-                maxY = y->data.F32[i];
-            }
-            if (y->data.F32[i] < minY) {
-                minY = y->data.F32[i];
-            }
-        }
+        rc = p_psVectorMin(y, NULL, 0, tmpStatsMinMax);
+        rc|= p_psVectorMax(y, NULL, 0, tmpStatsMinMax);
+        if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+            psFree(tmpMaskVec);
+            psFree(robustHistogram);
+            psFree(cumulativeRobustHistogram);
+            psFree(tmpStatsMinMax);
+            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
+            return(1);
+        }
+
         //
         // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
@@ -1879,5 +1879,5 @@
         //
         for (psS32 i = 0 ; i < y->n ; i++) {
-            y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY);
+            y->data.F32[i]= (y->data.F32[i] - tmpStatsMinMax->min) / (tmpStatsMinMax->max - tmpStatsMinMax->min);
         }
 
@@ -2295,5 +2295,5 @@
                         }
                     } else {
-                        // XXX: This if-statement really shouldn't be necessary.
+                        // This if-statement really shouldn't be necessary.
                         // However, due to numerical lack of precision, we
                         // occasionally produce a binNum outside the range.
@@ -2435,9 +2435,10 @@
 XXX: Should we free stats if the asserts fail?
  *****************************************************************************/
-psStats* psVectorStats(psStats* stats,
-                       const psVector* in,
-                       const psVector* errors,
-                       const psVector* mask,
-                       psMaskType maskVal)
+psStats* psVectorStats(
+    psStats* stats,
+    const psVector* in,
+    const psVector* errors,
+    const psVector* mask,
+    psMaskType maskVal)
 {
     psTrace(__func__, 3,"---- %s() begin  ----\n", __func__);