Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 26481)
+++ trunk/psLib/src/math/psStats.c	(revision 26615)
@@ -749,7 +749,5 @@
     // Iterate to get the best bin size; an iteration limit is enforced at the bottom of the loop.
     for (int iterate = 1; iterate > 0; iterate++) {
-        psTrace(TRACE, 6,
-                "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n",
-                iterate);
+        psTrace(TRACE, 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
 
         // Get the minimum and maximum values
@@ -791,12 +789,15 @@
         psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize);
 
-        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can not specify the bin
-        // size precisely since the argument to psHistogramAlloc() is the number of bins, not the binSize.  If
-        // we get here, we know that binSize != 0.0.
-        long numBins = (max - min) / binSize; // Number of bins
+        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can
+        // not specify the bin size precisely since the argument to psHistogramAlloc() is
+        // the number of bins, not the binSize.  If we get here, we know that binSize !=
+        // 0.0.  We can also have a floating-point round-off error such that the last bin
+        // of the histogram does not correspond exactly with the value of 'max'.  Let's be
+        // a bit generous and extend the histogram by two bins in either direction
+        long numBins = 4 + (max - min) / binSize; // Number of bins
         psTrace(TRACE, 6, "Numbins is %ld\n", numBins);
         psTrace(TRACE, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
         // Generate the histogram
-        histogram = psHistogramAlloc(min, max, numBins);
+        histogram = psHistogramAlloc(min - 2.0*binSize, max + 2.0*binSize, numBins);
         // XXXXX we need to consider this step if errors -> variance
         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
@@ -807,5 +808,4 @@
             psFree(statsMinMax);
             psFree(mask);
-
             return false;
         }
@@ -814,4 +814,37 @@
             PS_VECTOR_PRINT_F32(histogram->nums);
         }
+
+	// perversity check: if most of the values land in a single bin, then we probably
+	// have a perverse case (eg, small number of points at extremely large / small
+	// values; nearly bi-modal distribution).  if so, keep only points within 5? 10?
+	// bins of that excess bin:
+	int nMaxBin = 0;
+	int iMaxBin = 0;
+        for (long i = 1; i < histogram->nums->n; i++) {
+	    if (histogram->nums->data.F32[i] > nMaxBin) {
+		nMaxBin = histogram->nums->data.F32[i];
+		iMaxBin = i;
+	    }
+        }
+	if (nMaxBin > numValid / 2) {
+	    float minKeep = histogram->bounds->data.F32[iMaxBin] - 10*binSize;
+	    float maxKeep = histogram->bounds->data.F32[iMaxBin + 1] + 10*binSize;
+	    int nInvalid = 0;
+	    for (long i = 0; i < myVector->n; i++) {
+		// skip the already-masked values
+		if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskVal) continue;
+		bool invalid = false;
+		invalid |= (myVector->data.F32[i] < minKeep);
+		invalid |= (myVector->data.F32[i] > maxKeep);
+		invalid |= (!isfinite(myVector->data.F32[i]));
+		if (!invalid) continue;
+		mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal;
+		nInvalid ++;
+            }
+	    psTrace(TRACE, 6, "data is concentrated in a single bin, masking %d extreme outliers and retrying\n", nInvalid);
+            psFree(histogram);
+            psFree(cumulative);
+	    continue;
+	}
 
         // ADD step 1: Convert the specific histogram to a cumulative histogram
@@ -1007,4 +1040,5 @@
     stats->robustN50 = N50;
     psTrace(TRACE, 6, "The robustN50 is %ld.\n", N50);
+    psTrace(TRACE, 6, "The robust median and stdev are %f, %f\n", stats->robustMedian, stats->robustStdev);
 
     // Clean up