Index: /branches/eam_branch_20090208/psLib/src/math/psPolynomialUtils.c
===================================================================
--- /branches/eam_branch_20090208/psLib/src/math/psPolynomialUtils.c	(revision 21444)
+++ /branches/eam_branch_20090208/psLib/src/math/psPolynomialUtils.c	(revision 21445)
@@ -139,4 +139,5 @@
 // XXX determine the errors, and propagate to the output of psImageBicubeMin
 // XXX x,y are the image pixel index, but the fit is implied for the image pixel coordinates (0.5,0.5 center)
+// this solution is determined assuming constant weight per pixel
 psPolynomial2D *psImageBicubeFit(const psImage *image, int x, int y)
 {
@@ -179,5 +180,5 @@
 psPlane psImageBicubeMin(const psPolynomial2D *poly)
 {
-    psPlane min = { NAN, NAN, 0, 0 };   // Minimum value to return
+    psPlane min = { NAN, NAN, NAN, NAN };   // Minimum value to return
 
     PS_ASSERT_PTR_NON_NULL(poly, min);
@@ -186,6 +187,17 @@
 
     double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]); // Determinant
-    min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det;
-    min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det;
+    double xn = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]);
+    double yn = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]);
+
+    min.x = xn / det;
+    min.y = yn / det;
+
+    // without a supplied error, we calculate the normalized error
+    double detErr2 = 4.0*PS_SQR(poly->coeff[1][1])/6.0 + 8.0*PS_SQR(poly->coeff[0][2]) + 8.0*PS_SQR(poly->coeff[2][0]);
+    double xnErr2 = PS_SQR(poly->coeff[0][1])/6.0 + PS_SQR(poly->coeff[1][1])/6.0 + 2.0*PS_SQR(poly->coeff[1][0]) + 4.0*PS_SQR(poly->coeff[0][2])/6.0;
+    double ynErr2 = PS_SQR(poly->coeff[1][0])/6.0 + PS_SQR(poly->coeff[1][1])/6.0 + 2.0*PS_SQR(poly->coeff[0][1]) + 4.0*PS_SQR(poly->coeff[2][0])/6.0;
+
+    min.xErr = min.x * sqrt(xnErr2 / PS_SQR(xn) + detErr2 / PS_SQR(det));
+    min.yErr = min.y * sqrt(ynErr2 / PS_SQR(yn) + detErr2 / PS_SQR(det));
 
     return min;