Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 15048)
+++ trunk/psLib/src/math/psStats.c	(revision 15126)
@@ -13,6 +13,6 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.217 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-28 00:36:30 $
+ *  @version $Revision: 1.218 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-29 23:58:55 $
  *
  *  Copyright 2006 IfA, University of Hawaii
@@ -744,4 +744,5 @@
     long binLo, binHi;
     long binL2, binH2;
+    long binL4, binH4;
     long binMedian;
 
@@ -842,4 +843,6 @@
         PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, 0);
         PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, 0);
+        PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022481f, 0);
+        PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977519f, 0);
 
         psTrace(TRACE, 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n",
@@ -864,9 +867,11 @@
         // containing the value of interest.  constrain the answer to be within the current bin.
         // (extrapolation should not be needed and will result in errors)
-        float binLoF32, binHiF32, binL2F32, binH2F32;
-        PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.158655f);
-        PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.841345f);
+        float binLoF32, binHiF32, binL2F32, binH2F32, binL4F32, binH4F32;
+        PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo, totalDataPoints * 0.158655f);
+        PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binHi, totalDataPoints * 0.841345f);
         PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.308538f);
         PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.691462f);
+        PS_BIN_INTERPOLATE (binL4F32, cumulative->nums, cumulative->bounds, binL4, totalDataPoints * 0.022481f);
+        PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4, totalDataPoints * 0.977519f);
 
         // report +/- 1 sigma points
@@ -874,8 +879,24 @@
                 "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
                 binLoF32, binHiF32);
-
-        // ADD step 5: Determine SIGMA as (1/2 of) the distance between these positions.
-        // sigma = (binHiF32 - binLoF32) / 2.0;
-        sigma = (binH2F32 - binL2F32);
+        psTrace(TRACE, 5,
+                "The exact 30.8538 and 69.1462 percent data point positions are: (%f, %f)\n",
+                binL2F32, binH2F32);
+        psTrace(TRACE, 5,
+                "The exact 02.22481 and 97.7519 percent data point positions are: (%f, %f)\n",
+                binL4F32, binH4F32);
+
+        // ADD step 5: Determine SIGMA as the distance between binL2 and binH2 (+/- 0.5 sigma)
+        float sigma1 = (binH2F32 - binL2F32);
+        float sigma2 = (binHiF32 - binLoF32) / 2.0;
+        float sigma4 = (binH4F32 - binL4F32) / 4.0;
+
+	// take the smallest of the three: if we have a clump with wide outliers, sigma2 and
+	// sigma4 will be biased high; if we have a bi-modal distribution, sigma1 and sigma2
+	// will be biased high.
+	sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
+
+        psTrace(TRACE, 6, "The 2x sigma is %f.\n", sigma1);
+        psTrace(TRACE, 6, "The 2x sigma is %f.\n", sigma2);
+        psTrace(TRACE, 6, "The 4x sigma is %f.\n", sigma4);
 
         psTrace(TRACE, 6, "The current sigma is %f.\n", sigma);
@@ -1698,5 +1719,5 @@
         // XXX can we calculate the binMin, binMax **before** building this histogram?
 	// if the range is too absurd, adjust numBins & binSize
-        long numBins = PS_MAX (5, PS_MIN (10000, (max - min) / binSize));
+        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
 	binSize = (max - min) / (float) numBins;
         psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max);
@@ -1846,5 +1867,8 @@
 
 	// if we converge on a solution outside the range binMin - binMax, use a more conservative range
-	if (done && (guessMean > PS_BIN_MIDPOINT(histogram, binMax - 1))) {
+	float minValue = PS_BIN_MIDPOINT(histogram, binMin);
+	float maxValue = PS_BIN_MIDPOINT(histogram, binMax - 1);
+
+	if (done && ((guessMean < minValue) || (guessMean > maxValue))) {
             psTrace(TRACE, 6, "Inconsistent result, re-trying the fit\n");
 
@@ -1904,4 +1928,20 @@
             psTrace(TRACE, 6, "The symmetric stdev is %f.\n", guessStdev);
 #endif
+
+	    // if we converge on a solution outside the range binMin - binMax, use a more conservative range
+	    float minValueSym = PS_BIN_MIDPOINT(histogram, binS);
+	    float maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1);
+
+	    // saturate on min or max value
+	    if (guessMean < minValueSym) {
+		guessMean = minValueSym;
+		psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
+	    }
+		
+	    // saturate on min or max value
+	    if (guessMean > maxValueSym) {
+		guessMean = maxValueSym;
+		psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
+	    }
 
             psFree (poly);