Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 36114)
+++ trunk/psLib/src/math/psStats.c	(revision 36121)
@@ -879,5 +879,23 @@
         cumulative = psHistogramAlloc(min, max, numBins);
         cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
-
+        cumulative->bounds->data.F32[0] = histogram->bounds->data.F32[0];
+
+	// Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins).
+	long Nc = 1;  // track the current bin of cumulative
+	// the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the 
+	// upper end of the current histogram bin
+	for (long i = 1; i < histogram->nums->n; i++) {
+	    if (histogram->nums->data.F32[i] == 0.0) continue;
+	    cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i];
+	    cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i];
+	    Nc ++;
+	}
+	long Nlast = Nc - 1;  // last valid cumulative bin 
+	for (long i = Nc; i < histogram->nums->n; i++) { // Ensure the unused entries are filled.
+	    cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast];
+	    cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0;
+	}
+
+# if (0)
 	// Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins).
 	long delta = 0;
@@ -898,4 +916,5 @@
 	  cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0;
 	}
+# endif
 	
         if (psTraceGetLevel("psLib.math") >= 8) {
@@ -1051,20 +1070,31 @@
         if (sigma < (3.0 * binSize)) {
             psTrace(TRACE, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
-            long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
-            long maskHi = PS_MIN(histogram->bounds->n - 1, (binMedian + 25)); // High index for masking
-            psF32 medianLo = histogram->bounds->data.F32[maskLo]; // Value at low index
-            psF32 medianHi = histogram->bounds->data.F32[maskHi]; // Value at high index
+
+	    // these limits are supposed to be 25 x the raw bin size, NOT 25 of the cumulative histogram bins
+	    psF32 medianLo = stats->robustMedian - 25*binSize;
+	    psF32 medianHi = stats->robustMedian + 25*binSize;
+
+            // long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
+            // long maskHi = PS_MIN(cumulative->bounds->n - 1, (binMedian + 25)); // High index for masking
+            // psF32 medianLo = cumulative->bounds->data.F32[maskLo]; // Value at low index
+            // psF32 medianHi = cumulative->bounds->data.F32[maskHi]; // Value at high index
             psTrace(TRACE, 6, "Masking data more than 25 bins from the median\n");
-            psTrace(TRACE, 6,
-                    "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n",
-                    binMedian, maskLo, maskHi);
+            // psTrace(TRACE, 6, "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n", binMedian, maskLo, maskHi);
             psTrace(TRACE, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
+	    int Nmasked = 0;
             for (long i = 0 ; i < myVector->n ; i++) {
                 if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
-                    // XXXX is this correct?  is MASK_MARK safe?
+                    if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK) continue;
                     mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK;
                     psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]);
+		    Nmasked ++;
                 }
             }
+
+	    if (Nmasked == 0) {
+		// no significant change to the sigma & binsize -- we are done here 
+		iterate = -1;
+		continue;
+	    }
 
             // Free the histograms; they will be recreated on the next iteration, with new bounds
@@ -1108,5 +1138,5 @@
     PS_BIN_FOR_VALUE (binLo25, cumulative->nums, totalDataPoints * 0.25f, 0);
     PS_BIN_FOR_VALUE (binHi25, cumulative->nums, totalDataPoints * 0.75f, 0);
-    psTrace(TRACE, 6, "The 25-percent and 75-precent data point bins are (%ld, %ld).\n", binLo25, binHi25);
+    psTrace(TRACE, 6, "The 25-percent and 75-percent data point bins are (%ld, %ld).\n", binLo25, binHi25);
 
     // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
@@ -2530,16 +2560,30 @@
 		 yVec->data.F32[binNum - 0] +
 		 yVec->data.F32[binNum + 1] + yVec->data.F32[binNum + 2]);
-    double Sy = (xVec->data.F32[binNum - 2] + xVec->data.F32[binNum - 1] +
-		 xVec->data.F32[binNum - 0] +
-		 xVec->data.F32[binNum + 1] + xVec->data.F32[binNum + 2]);
+
+    double YM = 0.5*(xVec->data.F32[binNum - 2] + xVec->data.F32[binNum - 1]);
+    double Ym = 0.5*(xVec->data.F32[binNum - 1] + xVec->data.F32[binNum - 0]);
+    double Yo = 0.5*(xVec->data.F32[binNum - 0] + xVec->data.F32[binNum + 1]);
+    double Yp = 0.5*(xVec->data.F32[binNum + 1] + xVec->data.F32[binNum + 2]);
+    double YP = 0.5*(xVec->data.F32[binNum + 2] + xVec->data.F32[binNum + 3]);
+
+    double Sy = YM + Ym + Yo + Yp + YP;
+
     double Sxx = (PS_SQR(yVec->data.F32[binNum - 2]) + PS_SQR(yVec->data.F32[binNum - 1]) +
 		  PS_SQR(yVec->data.F32[binNum - 0]) +
 		  PS_SQR(yVec->data.F32[binNum + 1]) + PS_SQR(yVec->data.F32[binNum + 2]));
-    double Sxy = (xVec->data.F32[binNum - 2] * yVec->data.F32[binNum - 2] +
-		  xVec->data.F32[binNum - 1] * yVec->data.F32[binNum - 1] +
-		  xVec->data.F32[binNum - 0] * yVec->data.F32[binNum - 0] +
-		  xVec->data.F32[binNum + 1] * yVec->data.F32[binNum + 1] +
-		  xVec->data.F32[binNum + 2] * yVec->data.F32[binNum + 2]);
-    double value = ((Sy * Sxx - Sx * Sxy) + (5 * Sxy - Sx * Sy) * yVal)/(5 * Sxx - Sx * Sx);
+
+    double Sxy = (yVec->data.F32[binNum - 2] * YM +
+		  yVec->data.F32[binNum - 1] * Ym +
+		  yVec->data.F32[binNum - 0] * Yo +
+		  yVec->data.F32[binNum + 1] * Yp +
+		  yVec->data.F32[binNum + 2] * YP);
+
+    double Det = 5*Sxx - Sx*Sx;
+    if (Det == 0.0) return NAN;
+
+    double C0 = (Sy*Sxx - Sx*Sxy) / Det;
+    double C1 = (Sxy*5 - Sx*Sy) / Det;
+
+    double value = C0 + yVal*C1;
 
     return value;