Deriving the major and minor axis of an ellipse

[frankinstein]

1. Deriving the major and minor axis of an ellipse from conjugate diameters pointsCurrent solution:

double x1 = Diameter1.X - Center.X;
double y1 = Diameter1.Y - Center.Y;
double x2 = Diameter2.X - Center.X;
double y2 = Diameter2.Y - Center.Y;
double xc = (x2 + y1) / 2;
double yc = (y2 - x1) / 2;
double theta = Math.Atan2(yc - y2, x2 - xc);
double phi = Math.Atan2(yc, xc);
double d = xc / Math.Cos(phi);
double yu = yc + d * Math.Sin(theta);
double xv = xc + d * Math.Cos(theta);
double yv = yc - d * Math.Sin(theta);
double alpha = Math.Atan2(yv, xv);
a = (yu - y2) / Math.Sin(theta);
b = 0;
double SemiDiag1 = GraphicTools.GetRadius(Center, Diameter1);
double SemiDiag2 = GraphicTools.GetRadius(Center, Diameter2);
if (a != 0 && !double.IsNaN(a))
{
      b = Math.Sqrt((Math.Pow(SemiDiag1, 2) + Math.Pow(SemiDiag2, 2)) - Math.Pow(a, 2));
}
else
{
       b = d;
       a = Math.Sqrt((Math.Pow(SemiDiag1, 2) + Math.Pow(SemiDiag2, 2)) - Math.Pow(b, 2));
}

The above function where a = (yu - y2) / Math.Sin(theta) can be NaN when the conjugate radii's angles are: 0, 1.57, 3.14, 4.71, 6.28

All other angles the function works perfectly. The alternative process where a = 0 or a = NaN is inaccurate or completely dysfunctional. How does one treat the angles that don't work?

 

[Mark44]

1. Deriving the major and minor axis of an ellipse from conjugate diameters points


Current solution:

double x1 = Diameter1.X - Center.X;
double y1 = Diameter1.Y - Center.Y;
double x2 = Diameter2.X - Center.X;
double y2 = Diameter2.Y - Center.Y;
double xc = (x2 + y1) / 2;
double yc = (y2 - x1) / 2;
double theta = Math.Atan2(yc - y2, x2 - xc);
double phi = Math.Atan2(yc, xc);
double d = xc / Math.Cos(phi);
double yu = yc + d * Math.Sin(theta);
double xv = xc + d * Math.Cos(theta);
double yv = yc - d * Math.Sin(theta);
double alpha = Math.Atan2(yv, xv);
a = (yu - y2) / Math.Sin(theta);
b = 0;
double SemiDiag1 = GraphicTools.GetRadius(Center, Diameter1);
double SemiDiag2 = GraphicTools.GetRadius(Center, Diameter2);
if (a != 0 && !double.IsNaN(a))
{
      b = Math.Sqrt((Math.Pow(SemiDiag1, 2) + Math.Pow(SemiDiag2, 2)) - Math.Pow(a, 2));
}
else
{
       b = d;
       a = Math.Sqrt((Math.Pow(SemiDiag1, 2) + Math.Pow(SemiDiag2, 2)) - Math.Pow(b, 2));
}

The above function where a = (yu - y2) / Math.Sin(theta) can be NaN when the conjugate radii's angles are: 0, 1.57, 3.14, 4.71, 6.28

Whenever theta is $k\pi$, sin(theta) will be zero, so you'll have a division by zero problem at these values. You shouldn't have problems at $\pi/2$ (which you show above as 1.57) or $3\pi/2$ (which you show as 4.71), or at any odd multiple of $\pi/2$.

All other angles the function works perfectly. The alternative process where a = 0 or a = NaN is inaccurate or completely dysfunctional. How does one treat the angles that don't work?

The short answer is - don't use those angles. It would help if you explained what you were doing.

 

[frankinstein]

Whenever theta is $k\pi$, sin(theta) will be zero, so you'll have a division by zero problem at these values. You shouldn't have problems at $\pi/2$ (which you show above as 1.57) or $3\pi/2$ (which you show as 4.71), or at any odd multiple of $\pi/2$.


The short answer is - don't use those angles. It would help if you explained what you were doing.

Actually its the angle between the radii that's the problem and that's when theta is zero and if any of the x or y components zero out, like when the angle between the radii are ∏/2 or some multiple, causes the product of "a" to be zero or division by zero.

The function is used to generate an arc using the parametric equations for an ellipse, why the need for finding the major and minor axis from the conjugate diameter points. The parameters given to the function will have angles that result in products of "a" becoming unusable.

I was hoping that there was a better approach to extracting the major and minor axis from conjugate diameter points that didn't have the problems that my current solution has or a work around for those angles that produce an unusable "a".

 

[Mark44]

I don't know if you can work around the problem, but you can make your code more robust, by guarding against division by zero.

You'll potentially have problems with these lines of code:

double theta = Math.Atan2(yc - y2, x2 - xc);
double phi = Math.Atan2(yc, xc);
double d = xc / Math.Cos(phi);
// <snip>

double alpha = Math.Atan2(yv, xv);
a = (yu - y2) / Math.Sin(theta);

In the first line above, you should make sure that x2 - xc ≠ 0 or close to zero (say 10^-7 or so).
In the second line, you don't want xc to be zero or very close to zero.
In the third line, you don't want phi to be an odd multiple of pi/2.
In the last line, you don't want theta to be any multiple of pi.

For these special cases, I would see if there is some geometry I could exploit, like maybe one point is directly above the other, or on the same horizontal line or what-have-you, and take that into consideration.

BTW, I always avoided using Pow when I was just squaring a number. Rather than use Math.Pow(x, 2), I would just have x * x. I'm not sure the compiler is smart enough to convert a Pow() function call to a simple multiplication, which would be much faster.