5

Given a point p exterior to an axially aligned, origin centered ellipse E, find the (upto) four unique normals to E passing through p.

This is not a Mathematica question. Direct computation is too slow; I am willing to sacrifice precision and accuracy for speed.

I have searched the web, but all I found involved overly complex calculations which if implemented directly appear to lack the performance I need. Is there a more "programmatical" way to do this, like using matrices or scaling the ellipse into a circle?

andand
  • 17,134
  • 11
  • 53
  • 79
Timur Nuriyasov
  • 359
  • 2
  • 16
  • 3
    An ellipse has an uncountably infinite number of normals. To which 4 are you referring? Can you draw a picture of what you're looking for? – andand Dec 03 '13 at 20:52
  • 1
    This question appears to be off-topic because it's about a math problem and not about programming. It's better suited for http://math.stackexchange.com/ – datenwolf Dec 03 '13 at 20:55
  • @andand AFAIK, there can only be 4 normals maximum from a certain point not on the ellipse. – Timur Nuriyasov Dec 03 '13 at 21:16
  • @datenwolf yeah, I've seen similar questions on mathexchange, they usually end with an equation of some sort, which has to be solved for roots to get the answer. Now this is not what I need, because it's incredibly slow – Timur Nuriyasov Dec 03 '13 at 21:20
  • @TimurNuriyasov It sounds like you're using a non-standard definition of "normal". The standard definition of a normal is a vector perpendicular to a tangent at a particular point on the curve or surface (in the 3D case). – andand Dec 03 '13 at 21:21
  • @andand I'm not sure if you understand the problem. I am trying to find a normal. Perhaps this will help http://www.mathpages.com/home/kmath505/kmath505.htm – Timur Nuriyasov Dec 03 '13 at 21:25
  • 2
    The picture helps ... your question never described the relationship between the point and the specific normals you're interested in. – andand Dec 03 '13 at 21:34
  • 1
    @TimurNuriyasov: Unfortunately it will boil down to finding roots of the ellipse. Also applying some transformation won't help. However if you can live with an approximate solution, then a simple Newton method solver will quickly converge. (EDIT removed that part about elliptic curves, because it was not really accurate). – datenwolf Dec 03 '13 at 21:49
  • The website tells you what to do: the x-coordinates "of the points on the ellipse where the ellipse is perpendicular to the line through P are the real roots of `f`". "For practical purposes we can use Newton's method to quickly find a real root." Or is there a part of that you're not getting? – Teepeemm Dec 04 '13 at 04:05
  • @Teepemm Yeah, perhaps I have no choice but to do just that. I just thought there is a simpler and faster method that I overlooked, designed specifically for opengl and / or programmer use – Timur Nuriyasov Dec 04 '13 at 06:23
  • How your ellipse is defined? – MBo Dec 04 '13 at 09:46
  • @MBo I have major axis, minor axis and center location. But I can find almost any common thing about the ellipse – Timur Nuriyasov Dec 04 '13 at 11:34
  • Is it known whether the point is outside or inside the ellipse? Or are both cases of interest? The inside cases are slightly easier to code, since we know in advance there will be exactly four normals. The outside case may have two or four normals. – hardmath Dec 04 '13 at 11:50

3 Answers3

5

Let's assume the ellipse E is in "standard position", center at the origin and axes parallel to the coordinate axes:

    (x/a)^2 + (y/b)^2 = 1   where a > b > 0

The boundary cases a=b are circles, where the normal lines are simply ones that pass through the center (origin) and are thus easy to find. So we omit discussion of these cases.

The slope of the tangent to the ellipse at any point (x,y) may be found by implicit differentiation:

    dy/dx = -(b^2 x)/(a^2 y)

For the line passing through (x,y) and a specified point p = (u,v) not on the ellipse, that is normal to ellipse E when its slope is the negative reciprocal of dy/dx:

    (y-v)/(x-u) * (-b^2 x)/(a^2 y) = -1       (N)

which simplifies to:

    (x - (1+g)u) * (y + gv) = -g(1+g)uv  where g = b^2/(a^2 - b^2)

In this form we recognize it is the equation for a right rectangular hyperbola. Depending on how many points of intersection there are between the ellipse and the hyperbola (2,3,4), we have that many normals to E passing through p.

By reflected symmetry, if p is assumed exterior to E, we may take p to be in the first quadrant:

    (u/a)^2 + (v/b)^2 > 1    (exterior to E)
          u,v > 0            (1'st quadrant)

We could have boundary cases where u=0 or v=0, i.e. point p lies on an axis of E, but these cases may be reduced to solving a quadratic, because two normals are the (coinciding) lines through the endpoints of that axis. We defer further discussion of these special cases for the moment.

Here's an illustration with a=u=5,b=v=3 in which only one branch of the hyperbola intersects E, and there will be only two normals:

Ellipse with hyperbola

If the system of two equations in two unknowns (x,y) is reduced to one equation in one unknown, the simplest root-finding method to code is a bisection method, but knowing something about the possible locations of roots/intersections will expedite our search. The intersection in the first quadrant is the nearest point of E to p, and likewise the intersection in the third quadrant is the farthest point of E from p. If the point p were a good bit closer to the upper endpoint of the minor axis, the branches of the hyperbola would shift together enough to create up to two more points of intersection in the fourth quadrant.

One approach would be to parameterize E by points of intersection with the x-axis. The lines from p normal to the ellipse must intersect the major axis which is a finite interval [-a,+a]. We can test both the upper and lower points of intersection q=(x,y) of a line passing through p=(u,v) and (z,0) as z sweeps from -a to +a, looking for places where the ellipse and hyperbola intersect.

In more detail:

1. Find the upper and lower points `q` of intersection of E with the
   line through `p` and `(z,0)` (amounts to solving a quadratic)

3. Check the sign of a^2 y(x-u) - b^2 x(y-v) at `q=(x,y)`, because it
   is zero if and only `q` is a point of normal intersection

Once a subinterval is detected (either for upper or lower portion) where the sign changes, it can be refined to get the desired accuracy. If only modest accuracy is needed, there may be no need to use faster root finding methods, but even if they are needed, having a short subinterval that isolates a root (or root pair in the fourth quadrant) will be useful.

** more to come comparing convergence of various methods **

hardmath
  • 8,753
  • 2
  • 37
  • 65
  • If I understood correctly, this algo requires me to perform several iterations while looking for the derivative sign change. How many iterations should be performed? If it requires many iterations, what makes it better than finding four roots of the quartic c4x^4 + c3x^3 + c2x^2 + c1x + c0 = 0 with Newton's method (as is assumed here http://www.mathpages.com/home/kmath505/kmath505.htm)? – Timur Nuriyasov Dec 04 '13 at 12:37
  • These are good questions. Convergence of bisection is slow-but-sure while Newton's method is fast-but-unsure. So if you require (say) double precision answers, likely you want to start with a rough pigeonholing of roots using bisection (as outlined), then switch to a higher-order method (Newton is especially attractive for polynomials) once the root is "isolated" in a short interval that contains only that "single" root. But your Question led me to believe that simplicity of coding and modest precision were required. – hardmath Dec 04 '13 at 13:56
  • Note that sign change of (N)'s left hand side, rather than "derivative sign change", is sought. If (N) is multiplied out, you'll find it is a fourth degree polynomial in x. Perhaps I should add this to my Answer. – hardmath Dec 04 '13 at 14:09
  • 1
    @TimurNuriyasov: I haven't forgotten this! I'm writing up a detailed comparison of solving the quartic (per Kevin Brow/Mathpages) with Newton's method vs. the bisection method I propose. – hardmath Dec 07 '13 at 13:06
  • looking forward to it :) – Timur Nuriyasov Dec 09 '13 at 11:14
0

I had to solve a problem similar to this, for GPS initialization. The question is: what is the latitude of a point interior to the Earth, especially near the center, and is it single-valued? There are lots of methods for converting ECEF cartesian coordinates to geodetic latitude, longitude and altitude (look up "ECEF to Geodetic"). We use a fast one with only one divide and sqrt per iteration, instead of several trig evaluations like most methods, but since I can't find it in the wild, I can't give it to you here. I would start with Lin and Wang's method, since it only uses divisions in its iterations. Here is a plot of the ellipsoid surface normals to points within 100 km of Earth's center (North is up in the diagram, which is really ECEF Z, not Y): enter image description here

The star-shaped "caustic" in the figure center traces the center of curvature of the WGS-84 ellipsoid as latitude is varied from pole to equator. Note that the center of curvature at the poles is on the opposite side of the equator, due to polar flattening, and that the center of curvature at the equator is nearer to the surface than the axis of rotation.

Wherever lines cross, there is more than one latitude for that cartesian position. The green circle shows where our algorithm was struggling. If you consider that I cut off these normal vectors where they reach the axis, you would have even more normals for a given position for the problem considered in this SO thread. You would have 4 latitudes / normals inside the caustic, and 2 outside.

PaulQ
  • 308
  • 2
  • 9
0

The problem can be expressed as the solution of a cubic equation which gives 1, 2, or 3 real roots. For the derivation and closed form solution see Appendix B of Geodesics on an ellipsoid of revolution. The boundary between 1 and 3 solutions is an astroid.

cffk
  • 1,839
  • 1
  • 12
  • 17