shortest distance between a line and an ellipsoid

Question

after a week of fruitless trials I have to ask for help, most probably it is a trivial problem, although I could not find an answer using results of related problems on several forums. The problem I'm struggling with is as follows:

I'm looking for the shortest distance between a line and an ellipsoid. the ellipsoid E is defined by:

$x_p^TAx_p=1$ (1)

In which $A=diag[1/a^2,1/b^2,1/c^2]$ with a,b and c being half the length of the principal axes and $x_p$ is a point on the ellipsoid.

A line L is defined by $x_l=s_a+\lambda r_a$ (2)

with $x_l$ a point in the line.

Now, my general strategy is to find a point on the ellipsoid at which the surface normal is tangent to the directional vector $r_a$ of the line L and then work out the distance between that point (or rather points, as there will be two, logically) and the line as defined.

I followed 2 approaches:

1/ find a line parallel to L for which the dot product of the surface normal on the ellipsoid (which is $2x_p^TA$) and the directional vector of the line is zero. This boils down to finding a plane that is tangent to the ellipsoid and contains the directional vector of the line and a vector defined by the cross product of the surface normal and the directional vector of the line. Once this plane has been found it is trivial to work out the distance between the plane and the original line $x_l=s_a+\lambda r_a$. As far as I can see this should work but for some reason I cannot get correct results out of this (nor analytically nor numerical using either the Nelder-Mead or a Levenberg-Marquart algorithm.

2/ plug in $x_l=\alpha s_a+\beta p+\lambda r_a$ in equation (1), treat $\lambda$ as a unknown and work out the determinant of the resulting quadratic equation for $\lambda$ and find a combination of $\alpha$ and $\beta$ that set the discriminant of this equation to zero (resulting in a value for $\lambda$ that yields the value for a double root for the intersection between the line and the ellipsoid that is automatically fulfilling the demand to have a surface normal tangent to the directional vector of the line). Also following that line of thought did not result in correct answers. As I stated I could not use results from other related problems that are discussed on this and other forums as most of these answers assume that the point sought after and the line lie within a plane that goes through the origin, which is true for an ellipsoid with two equal axis, but not for the general case, as far as I can see).

As I'm new to this forum I apologise in advance for maybe not being 100% clear or posing a trivial question, but I'm really stuck and need someone with experience and/or a new perspective to look at this (due to the virus I'm confined in my house without anyone around having a clue about what I'm doing).

thanks in advance

Francois Clemens

Coming back to the reference (math.stackexchange.com/q/895452 ) using Plucker coordinates: I coded this option in Matlab: %supportvector sa

sa=[8,0,2]; %directional vector ra

ra=[2,2,1]; % axes of the ellipsoid

a=8; b=4; c=1; %define matrix A

A=zeros(3,3); A=diag([1/a^2,1/b^2,1/c^2]);

%derive the Plucker coordinates (e and p) in terms of vectors sa and ra

e=sa-ra; p=cross((sa+ra),sa);

% define vector r ("A 3D plane tangent to the ellipsoid through the closest point is perpendicular to the plane defined in (2) above. To find the normal of this plane, we use the vector describing the point on the line closest to the origin. This is given by"

r=-cross(e,cross(e,p))/(dot(e,e));

% and finally define the point on the ellipsoids sough for:

l=sqrt(dot(r*A^-1,r));

pos_ell=(r*A^-1)/l;

%the point on the line follows from:

ddv=(sa-pos_ell)-dot((sa-pos_ell),nr)*nr;

dd=sqrt(dot(ddv,ddv));

labda=dot((sa-pos_ell),nr);

pos_line=sa+labda*nr;

The attached figure shows a plot of the results: the blue line is line L, the red arrows represent the surface normal in the two intersection point between L and E and the green line connects the points that are supposed to represent the shortest distance between L and E according to the 'Plucker coordinate based solution'. Where do I go wrong?????

As far as I can see I made no coding errors, yet I obtain incorrect results:

Answer 1

Hint.

Given an ellipsoid and a line

$$ x^T A x = 1\\ x=x_0+\lambda v $$

relaxing the $x_0$ definition as $\bar x$ we have

$$ (\bar x + \lambda v)^TA(\bar x+\lambda v) = 1\Rightarrow\bar x^T A \bar x+2\lambda \bar x A v +\lambda^2v^T A v = 1 $$

solving for $\lambda$

$$ \lambda = \frac{-2\bar x A v\pm\sqrt{4(\bar x A v)^2-4(v^T A v)(\bar x^T A \bar x-1)}}{v^T A v} $$

at tangency we have

$$ (\bar x A v)^2-(v^T A v)(\bar x^T A \bar x-1)=0 $$

now choosing $\bar x$ such that $\bar x^T A \bar x = 1$ we have the condition

$$ \bar x A v = 0 $$

so $\bar x$ and $v$ should be orthogonal regarding $A$. Now the distances are computed as distances between $\bar x$ and the line $x=x_0+\lambda v$

NOTE

The intersection of the line $x = x_0+\lambda v$ with the plane $ x A v = 0$ is at $\lambda^* = -\frac{x_0^T A v}{v^T A v}$ or $x^* = x_0-\frac{x_0^T A v}{v^T A v}v$

This now is a planar problem that can be worked out in a plane as $O x y$. To do that we can rotate the plane $x^T A v$ around a convenient axis as for instance the axis formed by the line $\frac{v_x}{a^2}x + \frac{v_y}{b^2}y=0$ with direction $r = (-\frac{v_y}{b^2},\frac{v_x}{a^2},0)$ by an angle given by $\theta = \arccos\left(\frac{v_x}{||v||}\right)$. This can be done using the Rodrigues rotation formula (see). Now in the $Oxy$ plane we can compute the distance between the projected $x^*$ point and the projected ellipse.

NOTE

Included a MATHEMATICA script showing the adopted solution main lines

parms = {vx -> 2, vy -> 2, vz -> 1, a -> 8, b -> 4, c -> 1, x0 -> 8, y0 -> 0, z0 -> 2};
r = 11;
v = {vx, vy, vz}/Norm[{vx, vy, vz}];
zaxis = {0, 0, 1};
rpt = 0.1;
A = DiagonalMatrix[{1/a^2, 1/b^2, 1/c^2}];
p0 = {x0, y0, z0};
p = {x, y, z};
plane = p.A.v;
ellipsoid = p.A.p - 1;
past = p0 - (p0.A.v)/(v.A.v) v;
line = p0 + lambda v;
line0 = line /. parms;
linerot = plane /. {z -> 0};
axisrot = {-D[linerot, y], D[linerot, x], 0};
vrot = axisrot/Norm[axisrot];
theta = -ArcCos[zaxis.normal/Norm[normal]] /. parms;
R = RotationMatrix[theta, vrot] /. parms;
normal = Grad[plane, p];
rotellipsoid = p.R.A.Transpose[R].p - 1 /. parms // FullSimplify;
ellipsexy = rotellipsoid /. {z -> 0};
rotline = R.line /. parms;
vlrot = R.v;
p0rot = R.p0;
p0rot0 = p0rot /. parms;
pastrot = R.past; 
pastrot0 = pastrot /. parms;

and the graphics.

h = plane /. parms;
g = ellipsoid /. parms;
gr3 = ParametricPlot3D[line0, {lambda, -r, r}];
gr23 = ContourPlot3D[{h == 0, g == 0}, {x, -r, r}, {y, -r, r}, {z, -r, r},
MeshFunctions -> {Function[{x, y, z, f}, h - g]}, 
MeshStyle -> {{Thick, Blue}}, Mesh -> {{0}}, 
ContourStyle -> Directive[Orange, Opacity[0.2], Specularity[White, 30]], 
PlotPoints -> 60, SphericalRegion -> True];
grp0 = Graphics3D[{Red, Sphere[(p0 /. parms), rpt],  PointSize[0.08]}];
grpast = Graphics3D[{Black, Sphere[(past /. parms), rpt],  PointSize[0.08]}];
gr4 = ParametricPlot3D[lambda (vrot /. parms), {lambda, -r, r}];
h = z;
g = rotellipsoid;
gr34 = ContourPlot3D[{h == 0, g == 0}, {x, -r, r}, {y, -r, r}, {z, -r, r},
MeshFunctions -> {Function[{x, y, z, f}, h - g]}, 
MeshStyle -> {{Thick, Blue}}, Mesh -> {{0}}, 
ContourStyle -> Directive[Pink, Opacity[0.2], Specularity[White, 30]], 
PlotPoints -> 60, SphericalRegion -> True];
gr3rot = ParametricPlot3D[rotline, {lambda, -r, r}];
grp0rot = Graphics3D[{Red, Sphere[(p0rot0 /. parms), rpt], PointSize[0.08]}];
grpastrot = Graphics3D[{Black, Sphere[(pastrot0 /. parms), rpt],  PointSize[0.08]}];
Show[grp0, gr23, gr3, grpast, PlotRange -> {{-r, r}, {-r, r}, {-r, r}}]
Show[gr34, gr3rot, gr4, gr3rot, grp0rot, grpastrot, PlotRange -> {{-r, r}, {-r, r}, {-r, r}}]

and the distance calculations

rotline = p0rot + lambda vlrot;
L = (p - rotline).(p - rotline) + mu ellipsexy /. parms;
vars = Join[p, {lambda, mu}];
grad = Grad[L, vars];
sols = Solve[grad == 0, vars] // N;
pts = Table[(rotline /. parms /. sols[[k]]), {k, 1, Length[sols]}];
ptsz0 = Table[({x, y, 0} /. sols[[k]]), {k, 1, Length[sols]}];
grpts = Table[Graphics3D[{Blue, Sphere[pts[[k]], rpt]}], {k, 1, Length[sols]}];
grptsz0 = Table[Graphics3D[{Blue, Sphere[ptsz0[[k]], rpt]}], {k, 1, Length[sols]}];
grdists = Table[ParametricPlot3D[
pts[[k]] mu + ptsz0[[k]] (1 - mu), {mu, 0, 1}, 
PlotStyle -> Dashed], {k, 1, 4}];
Show[gr34, gr3rot, gr4, gr3rot,grpts, grptsz0, grdists, PlotRange -> {{-r, r}, {-r, r}, {-r, r}}]

Answer 2

Hint:

Looking in the direction of the line, you will see a single point and an ellipse, which is the apparent outline of the elliposid (and happens to be the cross-section of the ellipsoid by a plane).

For convenience, rotate space so that the line becomes parallel to the $z$ axis. You will find the apparent outline of the ellipsoid as the locus of the pairs $(x,y)$ such that the ellipsoid equation has a double root in $z$.

Now you are left with a 2D problem, i.e. the distance of a point to an ellipse. This is not very easy. Have a look here: https://stackoverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse