How do I find the tangent line to an ellipsoid nearest to a vector v(x, y, z) emanating from point p(x, y, z)?

Question

I'm trying the find the point on the tangent line horizon of an ellipsoid nearest to a vector emanating from a point p. Assuming the vector does not intersect the ellipse. Basically something like this:

Assume in that picture I know the location of point p and vector v. Also I know the radius components of the ellipse: 1 in the x direction and .99 in the y direction.

Answer 1

1. You can try to compute this algebraically

   but that will most likely lead to system of transcendental equations which we can not solve algebraically anyway unless some clever math tricks can be applied to convert them to more manageable form of equations. For example this is very similar problem (in 2D):

   • Increasing accuracy of solution of transcendental equation

   I was dealing a while ago and as you can see the equation is a bit problem. If you add the 3D it will become even worse not to mention ellipsoidal issues.

   So I strongly recommend not to go this way unless you have no other choice

2. I would use approximation search for the best point possible

   This is easy you just need to compute the distance of iterated solution to your real solution. First some assumptions:

   • p - camera position 3D
   • v - camera look vector 3D
   • r - axis aligned ellipsoid radiuses 3D
   • c - axis aligned ellipsoid center 3D
   • q - point where the tangent line is going through p (your wanted point)
   • q-p - tangent direction

   First we need to know the normal direction n(q) for any point q on the ellipsoid surface. On sphere it would be easy as it will be just q-c but for ellipsoid it is a bit more tricky. I am too lazy to make the math so instead i would use 2 close points to q not lying on the same line and use cross product. As you will use spherical coordinates for this anyway it is easy:

   q(a,b).x=c.x+r.x*cos(a)*cos(b)
   q(a,b).y=c.y+r.y*sin(a)*cos(b)
   q(a,b).z=c.z+r.z*sin(b)
   

   so 2 neighboring points can be q0=q(a+d,b),q1=q(a,b+d) where d is some small angle step constant. Now the normal is easy using cross product

   n(q)=cross(q0-q,q1-q) // perpendicular vectr
   n(q)/=|n(q)| // unit vector
   

   So finally the distance to solution metric e can be computed like this:

   e=|(dot(n(q),(q-p)/|q-p|)|
   

   If (q-p)/|q-p| and n(q) are perpendicular to each other then the dot product is 0.0 the more away from it the closer the result goes to +1.0 or -1.0 so if we use abs value of it then we can use it directly as search metrics.

   So now just try all points on the ellipsoid surface and remember the best solution (where e is minimal). That would be O(n^2) which would be insane on double precision but with use of approximation search or any other approximation technique it will become O(log^2(n)) which is acceptable.

   There is still one problem because this task has more than one solution so you need to throw away unwanted ones (like where the tangent to look vector has too big angle) In such cases just set the e to maximal value (1.0).

   You can further improve performance by not searching the whole surface. If the eccentricity is not too big then you can compute the q for sphere algebraically and then just search near point with same spherical angular coordinates on the ellipsoid.

   Another improvement can be axis separation ... you can convert this to 2x 2D problem solvable in O(log(n))

Answer 2

Two of the coordinate systems used in geodesy are (ECEF) cartesians and geographical coordinates. ECEF stands for earth-centred, earth-fixed. Given a point q outside the spheroid, it's geographical coordinates are the latitude and longitude of a point p on the spheroid, such that q is on the normal to the spheroid from p, and the height coordinate is then the distance along the normal from p to q.

Conversion from geographics to cartesians is straightforward, but going the other way requires some sort of iterative process.

For geographic to cartesian we have

x = r * cos( lambda)
y = r * sin( lambda)
z = ((1 - e2)*nu + h)*sin( phi)
where
r = (nu + h)*cos(phi)
nu = a / sqrt(1 - e2*sin(phi)*sin(phi))
phi the (geodetic) latitude, lambda the longitude, h the height 
(these names are pretty conventional in geodesy)

Going the other way it's easy enough to get lambda ( atan2( x, y)) and r ( hypot( x,y)) but disentangling phi and h from r and z is a bit tricky. Some methods are discussed here

One way to approach your problem would be to find the point on the camera line that has minimal height. The point with the same latitude and longitude and 0 height would be the point you seek. It might be advantageous to use a minimisation routione that usese derivates. It's straightforward -- though tedious -- to compute the derivative matrix of the transformation from geographics to cartesians at a point, and then to invert that to get the derivative matrix of the transformation from cartesians to geographics.

A thing to be beware of is that any calculation involving geographic coordinates is likely to blow up at the poles of the spheroid!