Changeset 37327

Timestamp:

Aug 29, 2014, 8:45:13 AM (12 years ago)

Author:

bills

Message:

Add psEllipsShapeToAxesWithErrors which propaaates the
errors in the shape parameters to the ellipse parameters.

Location:

trunk/psLib/src/math

Files:

• psEllipse.c (modified) (1 diff)
• psEllipse.h (modified) (1 diff)

trunk/psLib/src/math/psEllipse.c

@@ -134 +134 @@
     assert (isfinite(axes.major));
     assert (isfinite(axes.minor));
+
+    return axes;
+}
+// ellipse derotation (sx, sy, sxy) -> (major, minor, theta)
+psEllipseAxes psEllipseShapeToAxesWithErrors(psEllipseShape shape, psEllipseShape shapeErrors, double maxAR, psEllipseAxes *pShapeErrors)
+{
+    psEllipseAxes axes;
+    psEllipseAxes errors;
+
+    double f1 = 1.0 / PS_SQR(shape.sy) + 1.0 / PS_SQR(shape.sx);
+    double f2 = 1.0 / PS_SQR(shape.sy) - 1.0 / PS_SQR(shape.sx);
+
+    double f32 = PS_SQR(f2) + 4*PS_SQR(shape.sxy);
+    double f3 = sqrt(f32);
+
+    double df1 = - 2.0 * shapeErrors.sy / (PS_SQR(shape.sy) * shape.sy)
+                 - 2.0 * shapeErrors.sx / (PS_SQR(shape.sx) * shape.sx);
+    double df2 = - 2.0 * shapeErrors.sy / (PS_SQR(shape.sy) * shape.sy)
+                 + 2.0 * shapeErrors.sx / (PS_SQR(shape.sx) * shape.sx);
+
+    // df3 = 0.5 * df32 / sqrt(f32) = 0.5 * df32 / f3
+    double df3 = 0.5 * (2.0*f2*df2 + 2 * 4 * shape.sxy * shapeErrors.sxy) / f3;
+
+    double f13 = f1 + f3;
+    axes.minor = sqrt (2.0 / f13);
+
+    // dminor = -0.5 * sqrt (2) * df13 * f13**-3/2  = - 0.5 * df13 * minor / f13
+    errors.minor = - 0.5 * (df1 + df3) * axes.minor / f13;
+
+    axes.theta = -0.5 * atan2 (+2.0*shape.sxy, f2);
+
+    // according to wikipedia the derivitive of atan2(y, x) is
+    //
+    //  dAtan2(y, x) = - y * dx / (x**2 + y**2)) + x * dy / (x**2 + y**2) (for x > 0 and y!=0)
+    //  where
+    //  y = 2 * sxy  dy = 2 * dsxy   x = f2 so dx = df2
+
+    // dtheta = -0.5 * dAtan2(y, x)
+
+    errors.theta = -0.5 * ( -2.0 * shape.sxy * df2  +  f2 * 2 *shapeErrors.sxy ) / 
+                            (PS_SQR(f2) + 4 * PS_SQR(shape.sxy)) ;
+
+    // long, thin objects are likely to have a poorly measured major axis
+    // the angle and minor axis are likely to be ok.
+    // restrict the axis ratio
+    double rAR2 = (f1 - f3) / (f1 + f3);
+    if (rAR2 < 1.0/PS_SQR(maxAR)) {
+        axes.major = axes.minor * maxAR;
+        errors.major = errors.minor * maxAR;
+    } else {
+        axes.major = sqrt (2.0 / (f1 - f3));
+
+        // dmajor = -0.5 * (df1 - df3) * sqrt(2) * (f1 - f3)**-3/2
+        //        = -0.5 * (df2 - df3) * major / (f1 - f3)
+        errors.major = -0.5 * axes.major * (df1 - df3) / (f1 - f3);
+    }
+
+    assert (isfinite(axes.theta));
+    assert (isfinite(axes.major));
+    assert (isfinite(axes.minor));
+
+    *pShapeErrors = errors;
 
     return axes;

trunk/psLib/src/math/psEllipse.h

@@ -65 +65 @@
                                    );
 
+/// Convert shape to axes representation and compute errors
+psEllipseAxes psEllipseShapeToAxesWithErrors(psEllipseShape shape, ///< Shape of ellipse
+                                   psEllipseShape errors, ///<errors on shape params
+                                   double maxAR, ///< Maximum allowed axis ratio
+                                   psEllipseAxes *pShapeErrors ///< propagated errors on axes
+                                   );
+
 /// Convert axes to polarization representation
 psEllipsePol psEllipseAxesToPol(psEllipseAxes axes ///< Axes of ellipse