Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 36121)
+++ trunk/psLib/src/math/psStats.c	(revision 36222)
@@ -141,5 +141,5 @@
 
 // Debug information
-#define CZW 0
+#define CZW 1
 
 /*****************************************************************************/
@@ -232,4 +232,5 @@
         }
         count++;
+
     }
     if (errors) {
@@ -424,4 +425,5 @@
     for (long i = 0; i < myVector->n; i++) {
         // Check if the data is with the specified range
+
         if (!isfinite(data[i]))
             continue;
@@ -1043,6 +1045,6 @@
 
 
-#if (CZW)
-	printf("CZW (%d): median %f sigma %f delta: %f \t %f %f %f %f %f %f %f \t %f %f %f %f %f %f %f\n",
+#if (CZW && 0)
+	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,
 	   stats->robustMedian,stats->robustStdev,
@@ -1060,8 +1062,8 @@
 	   cumulative->nums->data.F32[binMedian+1],
 	   cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]);
-	PS_VECTOR_PRINT_F32(histogram->bounds);
-	PS_VECTOR_PRINT_F32(histogram->nums);
-	PS_VECTOR_PRINT_F32(cumulative->bounds);
-	PS_VECTOR_PRINT_F32(cumulative->nums);
+	//	PS_VECTOR_PRINT_F32(histogram->bounds);
+	//	PS_VECTOR_PRINT_F32(histogram->nums);
+	//	PS_VECTOR_PRINT_F32(cumulative->bounds);
+	//	PS_VECTOR_PRINT_F32(cumulative->nums);
 #endif
 
@@ -1276,5 +1278,5 @@
             return true;
         }
-	//	printf("BINS: %ld %f %f %f %f\n",stats->robustN50,guessStdev,binSize,min,max);
+
         // Calculate the number of bins.
         // XXX can we calculate the binMin, binMax **before** building this histogram?
@@ -1286,5 +1288,4 @@
         psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
 
-	//	printf("BINS2: %ld %f %f %f %f %ld\n",stats->robustN50,guessStdev,binSize,min,max,numBins);
 	
         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
@@ -1385,5 +1386,5 @@
             }
             y->n = x->n = j;
-
+	    
             // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
             // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
@@ -1397,4 +1398,5 @@
                 psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
+
                 psFree(poly);
                 psFree(histogram);
@@ -1478,4 +1480,12 @@
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+#if (CZW && 0)
+	    for (long i = 0; i < x->n; i++) {
+	      printf("CZW: Dcheck: %ld %f %f %f\n",
+		     i,x->data.F32[i],y->data.F32[i],
+		     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
+		     poly->coeff[2] * pow(x->data.F32[i],2));
+	    }
+#endif
             psFree(x);
             psFree(y);
@@ -1493,4 +1503,5 @@
             fullfitStdev = sqrt(-0.5/poly->coeff[2]);
             fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev);
+
 #ifndef PS_NO_TRACE
             psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
@@ -1515,4 +1526,5 @@
             }
 
+	    
             psFree (poly);
         }
@@ -1537,7 +1549,16 @@
 	    done = true;
 	}
-#if (CZW) 
-	printf("CZW IN FITTED? low %f %f full %f %f robust %f %f final %f %f\n",
-	       lowfitMean,lowfitStdev,fullfitMean,fullfitStdev,stats->robustMedian,stats->robustStdev,
+
+	
+#if (CZW && 0)
+	printf("CZW IN FITTED: iter   %d %f \n"
+	       "               low    %f %f \n"
+	       "               full   %f %f \n"
+	       "               robust %f %f \n"
+	       "               final  %f %f\n",
+	       iteration,minValueSym,
+	       lowfitMean,lowfitStdev,
+	       fullfitMean,fullfitStdev,
+	       stats->robustMedian,stats->robustStdev,
 	       guessMean,guessStdev);
 #endif
@@ -2423,12 +2444,4 @@
         psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
 
-#if (CZW)
-	printf("  polyin: %f %f %f == %f %f %f\n",
-	       x->data.F64[0], x->data.F64[1], x->data.F64[2],
-	       y->data.F64[0], y->data.F64[1], y->data.F64[2]);
-	printf("  rawpolyin: %f %f %f == %f %f %f\n",
-	       xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum + 1],
-	       y->data.F64[0], y->data.F64[1], y->data.F64[2]);
-#endif
 
         // Ensure that the y value lies within range of the y values.
@@ -2469,11 +2482,4 @@
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
-#if (CZW)
-	printf("  poly: %f %f %f fit: %f %f %f\n",
-	       myPoly->coeff[0],myPoly->coeff[1],myPoly->coeff[2],
-	       (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
-	       (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
-	       (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
-#endif
 
         psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
@@ -2498,8 +2504,4 @@
 	//        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
 	tmpFloat = binValue;
-		
-#if (CZW)
-	printf("   internal median: %f %f\n",tmpFloat,binValue);
-#endif
 	
     } else {
@@ -2556,36 +2558,39 @@
     )
 {
+
 # if (1)
-    double Sx = (yVec->data.F32[binNum - 2] + yVec->data.F32[binNum - 1] +
-		 yVec->data.F32[binNum - 0] +
-		 yVec->data.F32[binNum + 1] + yVec->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 = (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;
+# define HALF_SIZE 2
+  double Sx = 0.0;
+
+  double Sy = 0.0;
+  double Sxx = 0.0;
+  double Sxy = 0.0;
+  double deltaY = 0.0;
+  int N = 0;
+
+  for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) {
+    if ((u >= 0)&&(u < yVec->n)) {
+      if (u+1 < xVec->n) {
+	Sx += yVec->data.F32[u];
+	Sxx += PS_SQR(yVec->data.F32[u]);
+
+	deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
+	Sy += deltaY;
+	Sxy += yVec->data.F32[u] * deltaY;
+	N += 1;
+      }
+    }
+  }
+  double Det = N * Sxx - Sx * Sx;
+  if (Det == 0.0) return NAN;
+  if (N == 0) return NAN;
+
+  double C0 = (Sy*Sxx - Sx*Sxy) / Det;
+  double C1 = (Sxy*N - Sx*Sy) / Det;
+  
+  double value = C0 + yVal*C1;
+  return value;
+  
+  
 # else
     psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
@@ -2682,7 +2687,4 @@
 	tmpFloat = binValue;
 		
-#if (CZW)
-	printf("   internal median: %f %f\n",tmpFloat,binValue);
-#endif
 	
     } else {