Changeset 10254 for trunk/psLib/src/math/psEllipse.c

Timestamp:

Nov 28, 2006, 4:33:08 PM (19 years ago)

Author:

magnier

Message:

added axis-ratio limits to psEllipse functions, fixed math for psEllipse functions

File:

• trunk/psLib/src/math/psEllipse.c (modified) (4 diffs)

trunk/psLib/src/math/psEllipse.c

@@ -1 +1 @@
 #include <stdio.h>
+#include <assert.h>
 #include "psConstants.h"
 #include "psEllipse.h"
 
-psEllipseAxes psEllipseMomentsToAxes(psEllipseMoments moments)
+// ellipse rotation (major, minor, theta) -> (x2, y2, xy)
+psEllipseMoments psEllipseAxesToMoments(psEllipseAxes axes)
+{
+    psEllipseMoments moments;
+
+    double f1 = PS_SQR(axes.major) + PS_SQR(axes.minor);
+    double f2 = PS_SQR(axes.major) - PS_SQR(axes.minor);
+
+    moments.x2 = +0.5*f1 + 0.5*f2*cos(2*axes.theta);
+    moments.y2 = +0.5*f1 - 0.5*f2*cos(2*axes.theta);
+    moments.xy = -0.5*f2*sin(2*axes.theta);
+
+    assert (isfinite(moments.x2));
+    assert (isfinite(moments.y2));
+    assert (isfinite(moments.xy));
+
+    return moments;
+}
+
+// ellipse rotation (x2, y2, xy) -> (major, minor, theta)
+psEllipseAxes psEllipseMomentsToAxes(psEllipseMoments moments, double maxAR)
 {
     psEllipseAxes axes;
 
-    double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.xy));
-    if (f > (moments.x2 + moments.y2) / 2.0) {
-        f = 0.98*(moments.x2 + moments.y2) / 2.0;
+    double g1 = moments.x2 + moments.y2;
+    double g2 = moments.x2 - moments.y2;
+    double g3 = sqrt(PS_SQR(g2) + 4*PS_SQR(moments.xy));
+
+    axes.major = sqrt (0.5*(g1 + g3));
+    axes.theta = 0.5 * atan2 (-2.0*moments.xy, g2); // theta in radians
+
+    // long, thin objects are likely to have a poorly measured minor axis
+    // the angle and major axis are likely to be ok.
+    // restrict the axis ratio
+    double rAR2 = (g1 - g3) / (g1 + g3);
+    if (rAR2 < 1.0/PS_SQR(maxAR)) {
+        axes.minor = axes.major / maxAR;
+    } else {
+        axes.minor = sqrt (0.5*(g1 - g3));
     }
 
-    axes.major = sqrt (0.5*(moments.x2 + moments.y2) + f);
-    axes.minor = sqrt (0.5*(moments.x2 + moments.y2) - f);
-    axes.theta = atan2 (2*moments.xy, moments.x2 - moments.y2) / 2;
-    // theta in radians
+    assert (isfinite(axes.major));
+    assert (isfinite(axes.minor));
+    assert (isfinite(axes.theta));
 
     return axes;
@@ -25 +57 @@
     psEllipseShape shape;
 
-    double f1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor);
-    double f2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor);
+    double f1 = 1.0 / PS_SQR(axes.minor) + 1.0 / PS_SQR(axes.major);
+    double f2 = 1.0 / PS_SQR(axes.minor) - 1.0 / PS_SQR(axes.major);
 
-    double sxr = 0.5*f1 + 0.5*f2*cos(2*axes.theta);
-    double syr = 0.5*f1 - 0.5*f2*cos(2*axes.theta);
+    double sxr = 0.5*f1 - 0.5*f2*cos(2*axes.theta);
+    double syr = 0.5*f1 + 0.5*f2*cos(2*axes.theta);
 
-    shape.sx = 1.0 / sqrt(sxr);
-    shape.sy = 1.0 / sqrt(syr);
-    shape.sxy = 0.5*f2*sin(2*axes.theta);
+    shape.sx  = +1.0 / sqrt(sxr);
+    shape.sy  = +1.0 / sqrt(syr);
+    shape.sxy = -0.5*f2*sin(2*axes.theta);
+
+    assert (isfinite(shape.sx));
+    assert (isfinite(shape.sy));
+    assert (isfinite(shape.sxy));
 
     return shape;
@@ -39 +75 @@
 
 // ellipse derotation (sx, sy, sxy) -> (major, minor, theta)
-psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape)
+psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape, double maxAR)
 {
     psEllipseAxes axes;
 
-    double f1 = 1.0 / PS_SQR(shape.sx) + 1.0 / PS_SQR(shape.sy);
-    double f2 = 1.0 / PS_SQR(shape.sx) - 1.0 / PS_SQR(shape.sy);
+    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 f3 = sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
 
-    axes.theta = 0.5 * atan2 (2*shape.sxy, f2);
-    axes.major = sqrt (2.0 / (f1 - f3));
     axes.minor = sqrt (2.0 / (f1 + f3));
+    axes.theta = 0.5 * atan2 (-2.0*shape.sxy, f2);
+
+    // 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;
+    } else {
+        axes.major = sqrt (2.0 / (f1 - f3));
+    }
+
+    assert (isfinite(axes.theta));
+    assert (isfinite(axes.major));
+    assert (isfinite(axes.minor));
 
     return axes;
@@ -57 +106 @@
 // force the axis ratio to be less than 10
 // double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2));
+
+//    double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.xy));
+//    if (f > (moments.x2 + moments.y2) / 2.0) {
+//        f = 0.98*(moments.x2 + moments.y2) / 2.0;