IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 28, 2006, 10:12:11 AM (20 years ago)
Author:
magnier
Message:

fixed up the math on these; ADD text to follow

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/math/psEllipse.c

    r7766 r9768  
    2020}
    2121
     22// ellipse rotation (major, minor, theta) -> (sx, sy, sxy)
    2223psEllipseShape psEllipseAxesToShape(psEllipseAxes axes)
    2324{
    2425    psEllipseShape shape;
    2526
    26     double r1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor);
    27     double r2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor);
     27    double f1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor);
     28    double f2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor);
    2829
    29     double sxr = r1 + r2*cos(2*axes.theta);
    30     double syr = r1 - r2*cos(2*axes.theta);
     30    double sxr = 0.5*f1 + 0.5*f2*cos(2*axes.theta);
     31    double syr = 0.5*f1 - 0.5*f2*cos(2*axes.theta);
    3132
    3233    shape.sx = 1.0 / sqrt(sxr);
    3334    shape.sy = 1.0 / sqrt(syr);
    34     shape.sxy = r2*sin(2*axes.theta);
     35    shape.sxy = 0.5*r2*sin(2*axes.theta);
    3536
    3637    return shape;
    3738}
    3839
     40// ellipse derotation (sx, sy, sxy) -> (major, minor, theta)
    3941psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape)
    4042{
     
    4345    double f1 = 1.0 / PS_SQR(shape.sx) + 1.0 / PS_SQR(shape.sy);
    4446    double f2 = 1.0 / PS_SQR(shape.sx) - 1.0 / PS_SQR(shape.sy);
     47    double f3 = PS_SQR(f2) + 4*PS_SQR(shape.sxy);
    4548
    46     // force the axis ratio to be less than 10
    47     double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2));
    48 
    49     shape.sxy = PS_MIN(PS_MAX(shape.sxy, -r1), r1);
    50 
    51     axes.theta = atan2 (-2.0*shape.sxy, f2) / 2.0;
    52 
    53     double Ar = 0.25*f1 + 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
    54     double Br = 0.25*f1 - 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
    55 
    56     axes.minor = 1.0 / sqrt (Ar);
    57     axes.major = 1.0 / sqrt (Br);
     49    axes.theta = 0.5 * atan2 (-2.0*shape.sxy, f2) / 2.0;
     50    axes.major = sqrt (2.0 / (f1 - f3));
     51    axes.minor = sqrt (2.0 / (f1 + f3));
    5852
    5953    return axes;
    6054}
     55
     56// XXX keep this construction?
     57// force the axis ratio to be less than 10
     58// double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2));
Note: See TracChangeset for help on using the changeset viewer.