Changeset 21445

Timestamp:

Feb 11, 2009, 6:40:36 PM (17 years ago)

Author:

eugene

Message:

adding error calculation to psImageBicubeFit

File:

• branches/eam_branch_20090208/psLib/src/math/psPolynomialUtils.c (modified) (3 diffs)

branches/eam_branch_20090208/psLib/src/math/psPolynomialUtils.c

@@ -139 +139 @@
 // 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 +180 @@
 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 +187 @@
 
     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;