Index: /trunk/psLib/src/dataManip/psFunctions.c
===================================================================
--- /trunk/psLib/src/dataManip/psFunctions.c	(revision 2323)
+++ /trunk/psLib/src/dataManip/psFunctions.c	(revision 2324)
@@ -7,12 +7,11 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-04 01:04:59 $
+ *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-10 22:43:48 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  *
- *  XXX: What happens if the polyEVal functions are called with dat of the wrong
+ *  XXX: What happens if the polyEVal functions are called with data of the wrong
  *       type?
- *
  */
 /*****************************************************************************/
@@ -275,5 +274,5 @@
     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  
-    XXX: Should the "coeffErr[]" should be used as well?
+    XXX: Should the "coeffErr[]" be used as well?
  *****************************************************************************/
 static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly)
@@ -329,4 +328,5 @@
     psFree(d);
     return(tmp);
+
     /*
 
@@ -2017,6 +2017,4 @@
 XXX: This stuff does not currently work with a mask.
  
-XXX: nobody asked for us to generate this routine.
- 
 XXX: add another psScalar argument for the result.
  
Index: /trunk/psLib/src/dataManip/psStats.c
===================================================================
--- /trunk/psLib/src/dataManip/psStats.c	(revision 2323)
+++ /trunk/psLib/src/dataManip/psStats.c	(revision 2324)
@@ -9,6 +9,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.85 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-04 01:04:59 $
+ *  @version $Revision: 1.86 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-10 22:43:48 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -55,8 +55,4 @@
 psVector* p_psConvertToF32(psVector* in);
 
-int p_psVectorRobustStats(const psVector* restrict myVector,
-                          const psVector* restrict maskVector,
-                          psU32 maskVal,
-                          psStats* stats);
 /*****************************************************************************/
 /* GLOBAL VARIABLES                                                          */
@@ -129,6 +125,4 @@
 
 /******************************************************************************
- ******************************************************************************
- ******************************************************************************
     MISC PRIVATE STATISTICAL FUNCTIONS
  
@@ -141,8 +135,6 @@
         Is the in data structure NULL?
         Is the in data structure of type PS_TYPE_F32?
- ******************************************************************************
- ******************************************************************************
- *****************************************************************************/
-
+ 
+ *****************************************************************************/
 /******************************************************************************
 p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
@@ -155,6 +147,4 @@
 Returns
     NULL
-ASSUMPTION: the mean is always calculated exactly.  Robust means are never
-calculated in this routine.
  *****************************************************************************/
 
@@ -185,5 +175,10 @@
                 }
             }
-            mean /= (float)count;
+            if (count != 0) {
+                mean /= (float)count;
+            } else {
+                mean = NAN;
+            }
+
         } else {
             for (i = 0; i < myVector->n; i++) {
@@ -193,5 +188,9 @@
                 }
             }
-            mean /= (float)count;
+            if (count != 0) {
+                mean /= (float)count;
+            } else {
+                mean = NAN;
+            }
         }
     } else {
@@ -203,5 +202,9 @@
                 }
             }
-            mean /= (float)count;
+            if (count != 0) {
+                mean /= (float)count;
+            } else {
+                mean = NAN;
+            }
         } else {
             for (i = 0; i < myVector->n; i++) {
@@ -466,5 +469,4 @@
 
     // Calculate the median exactly.
-    // XXX: Is this the correct action?
     if (0 == (nValues % 2)) {
         stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues / 2) - 1] +
@@ -700,13 +702,18 @@
         }
     }
-    countFloat = (float)countInt;
-
-    #ifdef DARWIN
-
-    stats->sampleStdev = (float)sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
-    #else
-
-    stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
-    #endif
+    if (countInt < 2) {
+        stats->sampleStdev = NAN;
+        // XXX PS WARNING
+    } else {
+        countFloat = (float)countInt;
+        #ifdef DARWIN
+
+        stats->sampleStdev = (float)sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
+        #else
+
+        stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
+        #endif
+
+    }
 }
 
@@ -756,6 +763,4 @@
     }
     // 1. Compute the sample median.
-    // XXX: This seems odd.  Verify with IfA that we want to calculate the
-    // median here, not the mean.
     p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
 
@@ -828,4 +833,17 @@
     }
 
+    // Ensure that max!=min before we divide by (max-min)
+    if (FLT_EPSILON < fabs(max - min)) {
+        for (i = 0; i < myData->n; i++) {
+            myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * (outHigh - outLow) / (max - min);
+        }
+    } else {
+XXX:
+        PS_WARNING
+        for (i = 0; i < myData->n; i++) {
+            myData->data.F32[i] = outLow;
+        }
+    }
+
     for (i = 0; i < myData->n; i++) {
         myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * (outHigh - outLow) / (max - min);
@@ -849,6 +867,15 @@
     }
 
-    for (i = 0; i < myData->n; i++) {
-        myData->data.F64[i] = outLow + (myData->data.F64[i] - min) * (outHigh - outLow) / (max - min);
+    // Ensure that max!=min before we divide by (max-min)
+    if (FLT_EPSILON < fabs(max - min)) {
+        for (i = 0; i < myData->n; i++) {
+            myData->data.F64[i] = outLow + (myData->data.F64[i] - min) * (outHigh - outLow) / (max - min);
+        }
+    } else {
+XXX:
+        PS_WARNING
+        for (i = 0; i < myData->n; i++) {
+            myData->data.F64[i] = outLow;
+        }
     }
 }
@@ -1157,4 +1184,5 @@
         sumNfit += (float)robustHistogram->nums->data.U32[i];
     }
+    // XXX: divide by zero?
     myMean /= countFloat;
 
Index: /trunk/psLib/src/math/psPolynomial.c
===================================================================
--- /trunk/psLib/src/math/psPolynomial.c	(revision 2323)
+++ /trunk/psLib/src/math/psPolynomial.c	(revision 2324)
@@ -7,12 +7,11 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-04 01:04:59 $
+ *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-10 22:43:48 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  *
- *  XXX: What happens if the polyEVal functions are called with dat of the wrong
+ *  XXX: What happens if the polyEVal functions are called with data of the wrong
  *       type?
- *
  */
 /*****************************************************************************/
@@ -275,5 +274,5 @@
     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  
-    XXX: Should the "coeffErr[]" should be used as well?
+    XXX: Should the "coeffErr[]" be used as well?
  *****************************************************************************/
 static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly)
@@ -329,4 +328,5 @@
     psFree(d);
     return(tmp);
+
     /*
 
@@ -2017,6 +2017,4 @@
 XXX: This stuff does not currently work with a mask.
  
-XXX: nobody asked for us to generate this routine.
- 
 XXX: add another psScalar argument for the result.
  
Index: /trunk/psLib/src/math/psSpline.c
===================================================================
--- /trunk/psLib/src/math/psSpline.c	(revision 2323)
+++ /trunk/psLib/src/math/psSpline.c	(revision 2324)
@@ -7,12 +7,11 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-04 01:04:59 $
+ *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-10 22:43:48 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  *
- *  XXX: What happens if the polyEVal functions are called with dat of the wrong
+ *  XXX: What happens if the polyEVal functions are called with data of the wrong
  *       type?
- *
  */
 /*****************************************************************************/
@@ -275,5 +274,5 @@
     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  
-    XXX: Should the "coeffErr[]" should be used as well?
+    XXX: Should the "coeffErr[]" be used as well?
  *****************************************************************************/
 static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly)
@@ -329,4 +328,5 @@
     psFree(d);
     return(tmp);
+
     /*
 
@@ -2017,6 +2017,4 @@
 XXX: This stuff does not currently work with a mask.
  
-XXX: nobody asked for us to generate this routine.
- 
 XXX: add another psScalar argument for the result.
  
Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 2323)
+++ /trunk/psLib/src/math/psStats.c	(revision 2324)
@@ -9,6 +9,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.85 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-04 01:04:59 $
+ *  @version $Revision: 1.86 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-10 22:43:48 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -55,8 +55,4 @@
 psVector* p_psConvertToF32(psVector* in);
 
-int p_psVectorRobustStats(const psVector* restrict myVector,
-                          const psVector* restrict maskVector,
-                          psU32 maskVal,
-                          psStats* stats);
 /*****************************************************************************/
 /* GLOBAL VARIABLES                                                          */
@@ -129,6 +125,4 @@
 
 /******************************************************************************
- ******************************************************************************
- ******************************************************************************
     MISC PRIVATE STATISTICAL FUNCTIONS
  
@@ -141,8 +135,6 @@
         Is the in data structure NULL?
         Is the in data structure of type PS_TYPE_F32?
- ******************************************************************************
- ******************************************************************************
- *****************************************************************************/
-
+ 
+ *****************************************************************************/
 /******************************************************************************
 p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
@@ -155,6 +147,4 @@
 Returns
     NULL
-ASSUMPTION: the mean is always calculated exactly.  Robust means are never
-calculated in this routine.
  *****************************************************************************/
 
@@ -185,5 +175,10 @@
                 }
             }
-            mean /= (float)count;
+            if (count != 0) {
+                mean /= (float)count;
+            } else {
+                mean = NAN;
+            }
+
         } else {
             for (i = 0; i < myVector->n; i++) {
@@ -193,5 +188,9 @@
                 }
             }
-            mean /= (float)count;
+            if (count != 0) {
+                mean /= (float)count;
+            } else {
+                mean = NAN;
+            }
         }
     } else {
@@ -203,5 +202,9 @@
                 }
             }
-            mean /= (float)count;
+            if (count != 0) {
+                mean /= (float)count;
+            } else {
+                mean = NAN;
+            }
         } else {
             for (i = 0; i < myVector->n; i++) {
@@ -466,5 +469,4 @@
 
     // Calculate the median exactly.
-    // XXX: Is this the correct action?
     if (0 == (nValues % 2)) {
         stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues / 2) - 1] +
@@ -700,13 +702,18 @@
         }
     }
-    countFloat = (float)countInt;
-
-    #ifdef DARWIN
-
-    stats->sampleStdev = (float)sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
-    #else
-
-    stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
-    #endif
+    if (countInt < 2) {
+        stats->sampleStdev = NAN;
+        // XXX PS WARNING
+    } else {
+        countFloat = (float)countInt;
+        #ifdef DARWIN
+
+        stats->sampleStdev = (float)sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
+        #else
+
+        stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
+        #endif
+
+    }
 }
 
@@ -756,6 +763,4 @@
     }
     // 1. Compute the sample median.
-    // XXX: This seems odd.  Verify with IfA that we want to calculate the
-    // median here, not the mean.
     p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
 
@@ -828,4 +833,17 @@
     }
 
+    // Ensure that max!=min before we divide by (max-min)
+    if (FLT_EPSILON < fabs(max - min)) {
+        for (i = 0; i < myData->n; i++) {
+            myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * (outHigh - outLow) / (max - min);
+        }
+    } else {
+XXX:
+        PS_WARNING
+        for (i = 0; i < myData->n; i++) {
+            myData->data.F32[i] = outLow;
+        }
+    }
+
     for (i = 0; i < myData->n; i++) {
         myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * (outHigh - outLow) / (max - min);
@@ -849,6 +867,15 @@
     }
 
-    for (i = 0; i < myData->n; i++) {
-        myData->data.F64[i] = outLow + (myData->data.F64[i] - min) * (outHigh - outLow) / (max - min);
+    // Ensure that max!=min before we divide by (max-min)
+    if (FLT_EPSILON < fabs(max - min)) {
+        for (i = 0; i < myData->n; i++) {
+            myData->data.F64[i] = outLow + (myData->data.F64[i] - min) * (outHigh - outLow) / (max - min);
+        }
+    } else {
+XXX:
+        PS_WARNING
+        for (i = 0; i < myData->n; i++) {
+            myData->data.F64[i] = outLow;
+        }
     }
 }
@@ -1157,4 +1184,5 @@
         sumNfit += (float)robustHistogram->nums->data.U32[i];
     }
+    // XXX: divide by zero?
     myMean /= countFloat;