Index: /branches/rel-1_0/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- /branches/rel-1_0/psLib/src/math/psMinimizePolyFit.c	(revision 12337)
+++ /branches/rel-1_0/psLib/src/math/psMinimizePolyFit.c	(revision 12338)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-01-09 22:38:53 $
+ *  @version $Revision: 1.29.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-03-08 23:23:00 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -879,4 +879,16 @@
         double stdevValue = psStatsGetValue (stats, stdevOption);
 
+	// correct for the number of degrees of freedom absorbed by the fit
+	const int nterm = psPolynomial1DNterm(poly);
+	if (Nkeep <= nterm) {	// we can't estimate the standard deviation
+            psError(PS_ERR_UNKNOWN, true, "Number of data points %d is <= number of degrees of freedom %d",
+		    Nkeep, nterm);
+            psFree(resid);
+            psFree(fit);
+	    
+            return false;
+	}
+	stdevValue *= (Nkeep - 1)/(Nkeep - nterm - 1);
+
         psTrace("psLib.math", 5, "Mean is %f\n", meanValue);
         psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue);
@@ -1045,5 +1057,32 @@
     }
     psFree(xySums);
-
+#if 0
+    /*
+     * psMatrixGJSolve doesn't detect singular matrices very well. I'll do an SVN decomp
+     * to test for this -- but this is quite expensive, and a better solver would be wise. RHL.
+     *
+     * Don't run this test for now as I made the Gauss-Jordan solver require that |pivot| > FLT_EPSILON
+     *
+     * We could use lambda/evectors to solve the system without proceeding to a Gauss-Jordan
+     * solver, but I'm too lazy for now
+     */
+    psVector *lambda = psVectorAlloc(A->numCols, A->type.type); // eigenvalues
+    psImage *evectors = psMatrixSVD(NULL, lambda, A); // eigenvectors
+
+    for (int i = 0; i < lambda->n; i++) {
+	const double l = psVectorGet(lambda, i);
+
+	if (l > -FLT_EPSILON && l < FLT_EPSILON) {
+	    psError(PS_ERR_UNKNOWN, true, "Found almost zero eigenvalue %g.  Returning NULL.\n", l);
+	    psFree(lambda); psFree(evectors);
+	    psFree(A);
+	    psFree(B);
+	    return false;
+	}
+    }
+
+    psFree(lambda);
+    psFree(evectors);
+#endif
     if (!psMatrixGJSolve(A, B)) {
         psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
@@ -1116,5 +1155,5 @@
         result = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64);
         if (!result) {
-            psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
+            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
         }
         break;
@@ -1253,6 +1292,6 @@
         if (!psVectorStats(stats, resid, NULL, mask, maskValue)) {
             psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
-            psFree(resid)
-            psFree(fit)
+            psFree(resid);
+            psFree(fit);
             return false;
         }
@@ -1260,4 +1299,16 @@
         double meanValue = psStatsGetValue (stats, meanOption);
         double stdevValue = psStatsGetValue (stats, stdevOption);
+
+	// correct for the number of degrees of freedom absorbed by the fit
+	const int nterm = psPolynomial2DNterm(poly);
+	if (Nkeep <= nterm) {	// we can't estimate the standard deviation
+            psError(PS_ERR_UNKNOWN, true, "Number of data points %d is <= number of degrees of freedom %d",
+		    Nkeep, nterm);
+            psFree(resid);
+            psFree(fit);
+	    
+            return false;
+	}
+	stdevValue *= (Nkeep - 1)/(Nkeep - nterm - 1);
 
         psTrace("psLib.math", 5, "Mean is %f\n", meanValue);