Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 2273)
+++ 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;