Find Matrix Describing The 2D Ellipse Made by Intersecting a 3D Ellipsoid w/ a Plane

Question

I have a covariance matrix (P) describing an ellipsoid as well as an observer. The ellipsoid is centered at the origin and the observer is far away. I am fairly certain that from the observers perspective the ellipsoid is an ellipse and it is the outer perimeter made by projecting the ellipsoid onto a plane perpendicular to the line connecting the observer to the center of the ellipse.

In equations I have a semi-positive definite matrix P. Let L be a 3x3 matrix whose diagonal elements are the eigen values of P and v be a 3x3 matrix whose columns are the eigen vectors of P. The ellipsoid is then described by x^T A x = 1 with A = v^T L^-1 v. If I project this onto the plane using T (a simple 2x3 transform matrix) it produces a 2d ellipse of the form m^T B m = 1 where m are 2d vectors in the orthonormal basis of the plane. I want the matrix B.

An equivalent way I attempted to solve the problem geometrically is to find the 1D line characterized by the intersection of the ellipsoid and the plane. If I can parameterize it I can take the derivative and find the max/min of that in terms of distance from the line. Those points then define the major and minor axis of the projected ellipse and I'm set. Once I find the major point I can find the minor analytically since there's only one line perpendicular to the line of sight and the major axis and that line is easy to parameterize. But I can't find a way to parameterize the 2D ellipse in that plane itself.

A third equivalent question, which may be the simplest. I have a 3D ellipsoid and a line through the center. What are the two points (really 4 because symmetry) that are the farthest/closest to that line.

Any help on this would be appreciated as I'm really stuck finding the farthest point on the ellipsoid from the line analytically. I wrote a quick gradient descent solver but that's it.

For context, if it helps, I have a covariance describing the position uncertainty of an object. I want the covariance describing that uncertainty to an observer looking at the object. So a 2D covariance in az/el of the observer or in distance units in the focal plane of the observer. I can move between the two easily.

Answer 1

Since the observer is far away, we can assume the observer is at infinity.

Therefore what we're looking for is the orthogonal projection of the ellipsoid onto the plane $n^T p = d$ (we can take $d = 0$, this does not change the projection)

The ellipsoid itself is given by

$r^T Q r = 1$

where $r = [x,y,z]^T$ and $Q$ is a $3 \times 3$ symmetric and positive definite matrix.

We'll attach a reference frame $Ox'y'z'$ to the projection plane, with its $Ox'y'$ plane being the projection plane itself.

Using this reference frame, a point $r$ (in world coordinates) in space is described by

$r = R \ r'$

where $r' = \begin{bmatrix} x' \\ y' \\ z' \end{bmatrix} = \begin{bmatrix} p \\ t \end{bmatrix} $

where $p =\begin{bmatrix} x' \\ y' \end{bmatrix}$ is two-dimensional, and $z' = t $ with $t \in \mathbb{R}$.

The projection lines coming out of the projection plane are described by

$r = R \begin{bmatrix} p \\ t \end{bmatrix}$

Substitute this into the equation of the ellipsoid:

$\begin{bmatrix} p^T && t \end{bmatrix} R^T Q R \begin{bmatrix} p \\ t \end{bmatrix} = 1$

we now have to partition $ R^T Q R$ as follows

$ R^T Q R = \begin{bmatrix} A && b \\ b^T && c \end{bmatrix} $

where $A $ is $2 \times 2$, $b $ is $2 \times 1$, and $c$ is a scalar.

From the above equation, we now have

$p^T A p + 2 t \ p^T b + t^2 \ c = 1$

We want the discriminant of this quadratic in $t$ to be zero because we want a single root.

so

$ ( p^T b )^2 - (p^T A p - 1) c = 0 $

From which we finally get,

$ p^T ( {b b}^T - c A ) p + c = 0 $

Dividing through by $-c$ and re-arranging,

$ p^T B p = 1 $

where $ B = A - \dfrac{ { b b}^T}{c} $

And this is the equation of projection of the ellipsoid onto the plane.

There is a name for matrix $B$. It is known as the Schur's complement of $c$.

If, on the other hand, the observer is not at infinity, but at some point P, looking straight towards the center of the ellipsoid, then the analysis departs from the above. Attach a reference frame $Ox'y'z'$ with its origin at point $P$, and its $z'$ axis pointing away from world origin (which is the center of the ellipsoid) along the vector $P$. The $x'$ and $y'$ can be chosen freely. The world coordinate of a ray emanating from $P$ towards the ellipsoid intersecting the plane of projection at a point whose coordinates relative to the frame $O'x'y'z'$ is $ r' = \begin{bmatrix} p \\ -\alpha \end{bmatrix} $ where $\alpha$ is the fixed distance between the plane and the observer. The minus sign means the projection plane lies between the observation point and world origin (the center of the ellipsoid), because the $z'$ axis of the frame $O'x'y'z'$ is pointing away from the world origin.

The world coordinate corresponding to $r'$ is

$ r = P + R r' $

Now, the ray originating at $P$ and passing through $r$ has the parametric equation

$ q(t) = P + t (r - P) = P + t R r' $

Substitute this into the equation of the ellipsoid, and set the discriminat to zero as before.

$ ( P + t R r')^T Q (P + t R r' ) = 1 $

From which,

$ P^T Q P + 2 t P^T Q R r' + t^2 r'^T R^T Q R r' = 1 $

Therefore, by setting the discriminant to zero, we obtain the following expression

$ ( P^T Q R r' )^2 - ( P^T Q P - 1) ( r'^T R^T Q R r' ) = 0 $

Now note that

$P^T Q R r' = P^T {R R}^T Q R r' $

$P^T Q P = P^T {R R}^T Q {R R}^T P $

As before, partition $R^T Q R$ into

$ R^T Q R = \begin{bmatrix} A && b \\ b^T && c \end{bmatrix} $

And note that from the definition of $R$ as described above, we have,

$ P^T R = [ 0^T , f ] $ where $ f= \| P \| $.

Therefore, $P^T Q P - 1 = c f^2 - 1 $

Now, substitute $r'$ and $P^T Q P - 1$ and $P^T R$ you get

$ \left( [0^T , f] \begin{bmatrix} A && b \\ b^T && c \end{bmatrix} \begin{bmatrix} p \\ -\alpha \end{bmatrix} \right)^2 = (c f^2 - 1) [ p^T , -\alpha ]\begin{bmatrix} A && b \\ b^T && c \end{bmatrix} \begin{bmatrix} p \\ -\alpha \end{bmatrix} $

This simplifies to,

$ ( f b^T p - f c \alpha)^2 = (cf^2 -1) ( p^T A p - 2 \alpha b^T p + c \alpha^2 ) $

Hence,

$ p^T ( (c f^2 - 1) A - f^2 {bb}^T) p + 2 p^T b ( - \alpha (c f^2 - 1) + \alpha f^2 c) - f^2 c^2 \alpha^2 + (c f^2 - 1) ( c \alpha^2 ) = 0 $

This simplifies to

$ p^T ( (c f^2 - 1) A - f^2 {bb}^T ) p + 2 \ \alpha \ p^T b - c \alpha^2 = 0$

And this is the equation of the ellipse that results from intersecting the cone of sight of the ellipsoid with the projection plane.

I have numerically verified the correctness of this equation for a randomly selected $Q$ and $P$.

And for your third question, if you have the line $ \ell(t) = t \ V $ where $V$ is a direction vector. The minimum distance points are the points of intersection of this line with the ellipsoid, and there are two of them that are symmetric about the origin (which is the center of the ellipsoid). For the maximum distance points $r_1$, there is a corresponding value for the parameter $t$ (let's call it $t_1$) such that from perpendicularity we have

$ (r_1 - t_1 V) \cdot V = 0 $

and

$ (r_1 - t_1 V ) = K Q r_1 $

We can assume that $V$ is a unit vector, then

$ t_1 = V^T r_1 $

And it follows that

$ (r_1 - t_1 V) = (I - {V V}^T ) r_1 $

Hence we want

$(I - {VV}^T) r_1 = K Q r_1 $

The constant $K$ can be eliminated using the cross product operation and we get

$ \bigg( Q r_1 \bigg) \times \bigg( (I - {VV}^T ) r_1 \bigg) = \mathbf{0} $

And we need to solve this equation for $r_1$. For this, define the coordinate vectors $e_1 = [1, 0, 0]^T , e_2 = [0,1,0]^T , e_3 = [0,0,1]^T $. Further define the following matrices

$ E_1 = e_2 e_3^T - e_3 e_2^T $

$ E_2 = e_1 e_3^T - e_3 e_1^T $

$ E_3 = e_1 e_2^T - e_2 e_1^T $

Then the above vector equation becomes the following three scalar equations

$ r_1^T Q E_1 (I - {VV}^T) r_1 = 0 $

$ r_1^T Q E_2 (I - {VV}^T) r_1 = 0 $

$ r_1^T Q E_3 (I - {VV}^T) r_1 = 0 $

From these three equations, pick only two of them. In addition, we have the original equation

$ r_1^T Q r_1 = 1 $

Making a total of three quadratic equations in the coordinate vector $r_1$. Having formulated these equations, what remains is to solve them, and this can be done using any of several mathematical apps that are available online, such as Mathematica, or Sage.