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:

*

*Check that all five points are coplanar, and determine that plane.

*Convert the five 3D points to 2D points on that plane.

*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).

*Check the coefficients to make sure they represent an ellipse, and not another type of conic.

*Calculate the ellipse properties from the coefficients (formulae).

*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?

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$.
...
