Find the properties of an ellipse from 5 points in 3D space

Problem

I'd like to write code that solves the following problem:

• Take 5 arbitrary points in Cartesian coordinates $(x,y,z)$.
• Check whether there is an ellipse that goes through them all ^{(with some tolerance for floating-point inaccuracies)}.
• If so, find the ellipse's
  • center $\mathbf{c}$,
  • major radius $a$ (length of the semi-major axis),
  • minor radius $b$ (length of the semi-minor axis).

Canonical approach

Similar discussions tend to already start with five points in 2D space.
Extending that to 3D points, I suppose the "canonical" approach would look like this:

1. Check that all five points are coplanar, and determine that plane.
2. Convert the five 3D points to 2D points on that plane.
3. Use the five $(x,y)$ points to solve the conic equation $$ax^2+by^2+cxy+dx+ey+f=0$$ for the coefficients $a, b, c, d, e, f$ using some algorithm for solving linear system of equations ^{(hopefully stable with floating-point numbers)}.
4. Check the coefficients to make sure they represent an ellipse, and not another type of conic.
5. Calculate the ellipse properties from the coefficients ^(formulae).
6. Convert the ellipse center $\mathbf{c}$ back to 3D space.

Is there a shortcut?

The above approach seems cumbersome to implement, and possibly inefficient at run-time.
So, I'm wondering if there's a better way to do this in my case — where input and output is in 3D space, and I'm not actually interested in the full ellipse equation, just the three ellipse properties mentioned above.

I'm holding out hope, because for the simpler but conceptually similar problem of "finding the circle through three 3D points", this Wikipedia section provides a closed-formula solution with just a few dot products and cross products.

Any ideas?

Solution 1:

I don't know if this is simpler or not, but I'll propose a more geometric approach. All you need is a good routine to find the intersection point of lines $AB$ and $CD$ from the coordinates of points $A$, $B$, $C$, $D$.

Let $ABCDE$ be the five given points. I'll suppose in the following they are on the same plane. We can use Pascal's theorem to find the line tangent to the ellipse at $A$: if $F$ is the intersection of $AB$ with $CD$, and $G$ is the intersection of $AC$ with $BE$, then the intersection $T_A$ of $FG$ with $DE$ is a point on the tangent at $A$ to the conic $ABCDE$ (see here for a detailed proof).

Repeat the same construction with points $B$ and $C$, to find $T_B$ and $T_C$ on the respective tangents. Let then $H$ be the intersection of $AT_A$ with $BT_B$ and $L$ the intersection of $BT_B$ with $CT_C$. If $M$ and $N$ are the midpoints of segments $AB$ and $BC$, then the center $O$ of the ellipse is the intersection of lines $HM$ and $LN$ (this is due to a another well-wknon theorem: center, midpoint of a chord and intersection of the tangents at the endpoints of the chord are collinear).

Having found the center you can now compute the lengths of two conjugate semi-diameters: $$ \alpha=\sqrt{\overline{OH}\cdot \overline{OM}},\quad \beta={\alpha\cdot\overline{AM}\over\sqrt{\alpha^2-\overline{OM}^2}}, $$ while the angle $\theta$ between them is simply the angle between lines $HM$ and $AM$. If it is possible to make an ellipse pass through the five points then $\beta$ is a real number, but that is not a sufficient condition. One should also find point $K$, the intersection between line $OM$ and the line through $C$ parallel to $AB$: the ellipse can be constructed only if the number $$ \beta'={\alpha\cdot\overline{CK}\over\sqrt{\alpha^2-\overline{OK}^2}} $$ is the same as $\beta$.

You can finally find semi-axes $a$ and $b$ of the ellipse by solving the system: $$ a^2+b^2=\alpha^2+\beta^2\\ ab=\alpha\beta\sin\theta. $$

Solution 2:

The first question we have to settle is about the "foating point inaccuracies".

Let's work in homogeneous coordinates.

If the $5$ points are coplanar then there is a normal vector to the plane $\bf n$ such that $$ \left( {\matrix{ {x_{\,1} } & {y_{\,1} } & {z_{\,1} } & 1 \cr {x_{\,2} } & {y_{\,2} } & {z_{\,2} } & 1 \cr \vdots & \vdots & \vdots & \vdots \cr \vdots & \vdots & \vdots & \vdots \cr {x_{\,5} } & {y_{\,5} } & {z_{\,5} } & 1 \cr } } \right) \left( {\matrix{ {n_{\,1} } \cr {n_{\,2} } \cr {n_{\,3} } \cr {n_{\,0} } \cr } } \right) = \left( {\matrix{ 0 \cr 0 \cr 0 \cr 0 \cr 0 \cr } } \right) $$ That means that $\bf n$ is in the null-space of $\bf X$, and that this has rank (no greater than) $3$.

The above construction is difficult to check in presence of "inaccuracies".

I would approach the problem by finding $\bf n$ through a least-squares interpolation, i.e. by minimizing the modulus of the vector resulting by the multiplication above.
This task is greatly simplified by taking the Singular Value Decomposition of $\bf X$.
...