Index: trunk/psLib/src/math/psPolynomialUtils.c
===================================================================
--- trunk/psLib/src/math/psPolynomialUtils.c	(revision 21183)
+++ trunk/psLib/src/math/psPolynomialUtils.c	(revision 21518)
@@ -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,19 @@
 
     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;
+
+    // sigma_xn^2 / xn^2 = sigma_11^2 / C11
+
+    // without a supplied error, we calculate the normalized error
+    double fdetErr2 = 0.5/PS_SQR(poly->coeff[0][2]) + 0.5/PS_SQR(poly->coeff[2][0]) + 4.0/(6.0*PS_SQR(poly->coeff[1][1]));
+    double fxnErr2 = 1.0/(6.0*PS_SQR(poly->coeff[0][1])) + 1.0/(6.0*PS_SQR(poly->coeff[1][0])) + 1.0/(6.0*PS_SQR(poly->coeff[1][1])) + 0.5/PS_SQR(poly->coeff[0][2]);
+    double fynErr2 = 1.0/(6.0*PS_SQR(poly->coeff[0][1])) + 1.0/(6.0*PS_SQR(poly->coeff[1][0])) + 1.0/(6.0*PS_SQR(poly->coeff[1][1])) + 0.5/PS_SQR(poly->coeff[2][0]);
+
+    min.xErr = fabs(min.x) * sqrt(fxnErr2 + fdetErr2);
+    min.yErr = fabs(min.y) * sqrt(fynErr2 + fdetErr2);
 
     return min;