Index: trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- trunk/psLib/src/math/psMinimizePolyFit.c	(revision 36222)
+++ trunk/psLib/src/math/psMinimizePolyFit.c	(revision 36237)
@@ -770,4 +770,7 @@
 #endif
 	}
+	psFree(Zcoeff);
+	psFree(ZcoeffErr);
+
       }
 	
Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 36222)
+++ trunk/psLib/src/math/psStats.c	(revision 36237)
@@ -141,5 +141,5 @@
 
 // Debug information
-#define CZW 1
+#define CZW 0
 
 /*****************************************************************************/
@@ -425,5 +425,4 @@
     for (long i = 0; i < myVector->n; i++) {
         // Check if the data is with the specified range
-
         if (!isfinite(data[i]))
             continue;
@@ -881,5 +880,5 @@
         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];
+        cumulative->bounds->data.F32[0] = histogram->bounds->data.F32[1];
 
 	// Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins).
@@ -887,8 +886,8 @@
 	// 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++) {
+	for (long i = 1; i < histogram->nums->n - 1; 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];
+	    cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i+1];
 	    Nc ++;
 	}
@@ -898,25 +897,4 @@
 	    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;
-	long delta_x = 0;
-	for (long i = 1; i < histogram->nums->n; i++) {
-	  if (histogram->nums->data.F32[i] > 0.0) { // This bin is bigger than the last one.
-	    cumulative->nums->data.F32[i - delta] = cumulative->nums->data.F32[i - delta - 1] + histogram->nums->data.F32[i];
-	    cumulative->bounds->data.F32[i - delta - 1] = histogram->bounds->data.F32[i - delta_x];
-	    delta_x = 0;
-	  }
-	  else { // This bin is the same as the last one, so we shouldn't count this bound
-	    delta++;
-	    delta_x++;
-	  }
-	}
-	for (long i = histogram->nums->n - delta; i < histogram->nums->n; i++) { // Ensure the unused entries are filled.
-	  cumulative->nums->data.F32[i] = cumulative->nums->data.F32[histogram->nums->n - delta - 1];
-	  cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0;
-	}
-# endif
 	
         if (psTraceGetLevel("psLib.math") >= 8) {
@@ -940,25 +918,4 @@
 	// There's no reason to do a quadratic fit near the 50% bin, as it's approximately linear there.
 	// Instead, do a 5-point linear fit.
-# if (0)
-	{ // Quick 5-point linear fit
-	  double Sx = (cumulative->nums->data.F32[binMedian - 2] + cumulative->nums->data.F32[binMedian - 1] +
-		       cumulative->nums->data.F32[binMedian - 0] +
-		       cumulative->nums->data.F32[binMedian + 1] + cumulative->nums->data.F32[binMedian + 2]);
-	  double Sy = (cumulative->bounds->data.F32[binMedian - 2] + cumulative->bounds->data.F32[binMedian - 1] +
-		       cumulative->bounds->data.F32[binMedian - 0] +
-		       cumulative->bounds->data.F32[binMedian + 1] + cumulative->bounds->data.F32[binMedian + 2]);
-	  double Sxx = (pow(cumulative->nums->data.F32[binMedian - 2],2) + pow(cumulative->nums->data.F32[binMedian - 1],2) +
-			pow(cumulative->nums->data.F32[binMedian - 0],2) +
-			pow(cumulative->nums->data.F32[binMedian + 1],2) + pow(cumulative->nums->data.F32[binMedian + 2],2));
-	  double Sxy = (cumulative->bounds->data.F32[binMedian - 2] * cumulative->nums->data.F32[binMedian - 2] +
-			cumulative->bounds->data.F32[binMedian - 1] * cumulative->nums->data.F32[binMedian - 1] +
-			cumulative->bounds->data.F32[binMedian - 0] * cumulative->nums->data.F32[binMedian - 0] +
-			cumulative->bounds->data.F32[binMedian + 1] * cumulative->nums->data.F32[binMedian + 1] +
-			cumulative->bounds->data.F32[binMedian + 2] * cumulative->nums->data.F32[binMedian + 2]);
-	  double linearMedian = ((Sy * Sxx - Sx * Sxy) + (5 * Sxy - Sx * Sy) * (totalDataPoints/2.0))/(5 * Sxx - Sx * Sx);
-	  // psLogMsg("psLib",5,"Median Comp: %f %f\n",stats->robustMedian,linearMedian);
-	  stats->robustMedian = linearMedian;
-	}
-# endif
         stats->robustMedian = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
 
@@ -1003,4 +960,5 @@
         // (extrapolation should not be needed and will result in errors)
         float binLoF32, binHiF32, binL2F32, binH2F32, binL4F32, binH4F32;
+#if (0)
         PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo,
                             totalDataPoints * 0.158655f);
@@ -1015,5 +973,12 @@
         PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4,
                             totalDataPoints * 0.977250f);
-
+#else
+        binLoF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo, totalDataPoints * 0.158655);
+	binHiF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi, totalDataPoints * 0.841345);	       
+        binL2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL2, totalDataPoints * 0.308538);
+	binH2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH2, totalDataPoints * 0.691462);	       
+        binL4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL4, totalDataPoints * 0.022750);
+	binH4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH4, totalDataPoints * 0.977250);
+#endif	
         // report +/- 1 sigma points
         psTrace(TRACE, 5,
@@ -1027,13 +992,36 @@
                 binL4F32, binH4F32);
 
+	// If some of the fits failed, attempt to fix this
+	if (!isfinite(binLoF32) && isfinite(binHiF32)) { binLoF32 = -1.0 * binHiF32; }
+	if (!isfinite(binHiF32) && isfinite(binLoF32)) { binHiF32 = -1.0 * binLoF32; }
+	if (!isfinite(binL2F32) && isfinite(binH2F32)) { binL2F32 = -1.0 * binH2F32; }
+	if (!isfinite(binH2F32) && isfinite(binL2F32)) { binH2F32 = -1.0 * binL2F32; }
+	if (!isfinite(binL4F32) && isfinite(binH4F32)) { binL4F32 = -1.0 * binH4F32; }
+	if (!isfinite(binH4F32) && isfinite(binL4F32)) { binH4F32 = -1.0 * binL4F32; }
+	
         // 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;
 
+	// Fix again?
+	if (!isfinite(sigma1) && isfinite(sigma2) && isfinite(sigma4)) { sigma1 = (sigma2 + sigma4) / 2.0; }
+	if (!isfinite(sigma2) && isfinite(sigma1) && isfinite(sigma4)) { sigma2 = (sigma1 + sigma4) / 2.0; }
+	if (!isfinite(sigma4) && isfinite(sigma2) && isfinite(sigma1)) { sigma4 = (sigma2 + sigma1) / 2.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));
+	//        sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
+	// CZW: Instead, take the median.  Taking the MIN forces a bias on unbiased data.
+	//      It seems like occasionally getting the wrong answer on a complex distribution
+	//      is more acceptable than always getting the wrong answer for simple ones.
+
+	
+	sigma = PS_MAX( PS_MIN(sigma1,sigma2),
+			PS_MIN( PS_MAX(sigma1,sigma2),
+				sigma4));
 
         psTrace(TRACE, 6, "The 1x sigma is %f.\n", sigma1);
@@ -1042,8 +1030,27 @@
 
         psTrace(TRACE, 6, "The current sigma is %f.\n", sigma);
-        stats->robustStdev = sigma;
-
-
-#if (CZW && 0)
+	//        stats->robustStdev = sigma;
+	stats->robustStdev = sigma;
+
+#if (CZW && 0) 
+	// Skewness check: Find least biased sample for each pair.
+	sigma1 = 2.0 * PS_MIN(binH2F32 - stats->robustMedian,
+			      stats->robustMedian - binL2F32);
+	sigma2 = 1.0 * PS_MIN(binHiF32 - stats->robustMedian,
+			      stats->robustMedian - binLoF32);
+	sigma4 = 0.5 * PS_MIN(binH4F32 - stats->robustMedian,
+			      stats->robustMedian - binL4F32);
+	// Kurtosis check: Take median sample as the solution.
+	stats->robustStdev = PS_MAX( PS_MIN(sigma1,sigma2),
+				     PS_MIN( PS_MAX(sigma1,sigma2),
+					     sigma4));
+#endif
+
+	
+#if (0)
+	printf("CZW: bad sigma?: %f %f  %f %f  %f %f  %f %f %f  %f\n",
+	       binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32,
+	       sigma1,sigma2,sigma4,sigma);
+	
 	printf("CZW (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
 	       iterate,
@@ -1200,5 +1207,5 @@
  * "vectorFittedStats_v4" all versions of fitted stats now resolve to this function (only v4
  * has really been used) vectorFittedStats requires guess for fittedMean and fittedStdev
- * robustN50 should also be set gaussian fit is performed using 2D polynomial to ln(y) this
+ * robustN50 should also be set gaussian fit is performed using 1D polynomial to ln(y) this
  * version follows the upper portion of the distribution until it passes 0.5*peak
  ********************/
@@ -1235,5 +1242,5 @@
         return true;
     }
-
+    if (myVector->n < 1) { printf("There are no elements in this vector.\n"); abort(); }
     float guessStdev = stats->robustStdev;  // pass the guess sigma
     float guessMean = stats->robustMedian;  // pass the guess mean
@@ -1255,6 +1262,6 @@
             // set roughly so that the lowest bins have about 2 cnts
             // Nsmallest ~ N50 / (4*dN))
-            psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8));
-            binSize = guessStdev / dN;
+            psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8)); 
+           binSize = guessStdev / dN;
         }
 
@@ -1288,5 +1295,5 @@
         psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
 
-	
+
         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
@@ -1551,5 +1558,5 @@
 
 	
-#if (CZW && 0)
+#if (CZW && 1)
 	printf("CZW IN FITTED: iter   %d %f \n"
 	       "               low    %f %f \n"
@@ -2575,5 +2582,6 @@
 	Sxx += PS_SQR(yVec->data.F32[u]);
 
-	deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
+	deltaY = xVec->data.F32[u];
+	//deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
 	Sy += deltaY;
 	Sxy += yVec->data.F32[u] * deltaY;