Solving a least squares problem to try and get a best fit ellipsoid (or a general 4D conic section) to a cloud of data points such as the database of nearby stars.
Can solve the linear equation portion, getting a result in the form of:
$a_0x^2 + a_1y^2 + a_2z^2 +a_3xy + a_4xz + a_5yz + a_6x + a_7y + a_8z + 1 = 0$
This question has similarities to: Ellipsoid Equation: converting so I tried to work through that first.
From there, hoping to get a parametric representation similar to General Ellipse Equations
where:
$x_c = $ ellipsoid center in x
$y_c = $ ellipsoid center in y
$z_c = $ ellipsoid center in y
$a = $ ellipsoid major axis length
$b = $ ellipsoid mid axis length
$c = $ ellipsoid minor axis length
$α = $ ellipsoid yaw
$β = $ ellipsoid pitch
$γ = $ ellipsoid roll
$θ = $ parameter, which ranges from 0 to 2π radians
$ϕ = $ parameter, which ranges from π/2 to -π/2 radians
What I've Done
Based on the answers in the previous question, the process looks like convert to:
$ \mathbf x^T M \mathbf x + N \mathbf x + c = 0 $
with:
$\mathbf x = \begin{pmatrix} x \\\\ y \\\\ z \end{pmatrix} $
$ M = \begin{pmatrix} A & D & E \\\\ D & B & F \\\\ E & F & C \end{pmatrix} = \begin{pmatrix} a_0 & \frac{1}{2}a_3 & \frac{1}{2}a_4 \\\\ \frac{1}{2}a_3 & a_1 & \frac{1}{2}a_5 \\\\ \frac{1}{2}a_4 & \frac{1}{2}a_5 & a_2 \end{pmatrix} $
$N = \begin{pmatrix} a_6 \\\\ a_7 \\\\ a_8 \end{pmatrix} = \begin{pmatrix} G \\\\ H \\\\ I \end{pmatrix} $
Then, can find the center:
$\mathbf x_c = \begin{pmatrix} x_c \\\\ y_c \\\\ z_c \end{pmatrix} $
Using: $\mathbf x_c = − \frac{1}{2}M^{-1}N$
Doing so, I get:
$\mathbf x_c = \begin{pmatrix} x_c \\\\ y_c \\\\ z_c \end{pmatrix} = - \frac{1}{2} \begin{pmatrix} \frac{G\left(BC-F^2\right)+H\left(EF-CD\right)+I\left(DF-BE\right)}{2DEF+C\left(BA-D^2\right)-AF^2-BE^2} \\\\ \frac{G\left(EF-CD\right)+H\left(CA-E^2\right)+I\left(DE-AF\right)}{2DEF+C\left(BA-D^2\right)-AF^2-BE^2} \\\\ \frac{G\left(DF-BE\right)+H\left(DE-AF\right)+I\left(BA-D^2\right)}{2DEF+C\left(BA-D^2\right)-AF^2-BE^2} \end{pmatrix}$
From @Jan-Magnus Økland answer in the same thread, and this Wikipedia on Ellipsoids it looks like a,b,c are then the reciprocals of the eigenvalues of M found with $M=R D R^T$ where:
$D = \begin{bmatrix} \lambda _1 & 0 & 0 \\\\ 0 & \lambda _2 & 0 \\\\ 0 & 0 & \lambda _3 \end{bmatrix} $
$a = \frac{1}{\sqrt{- \lambda _1}}$
$b = \frac{1}{\sqrt{- \lambda _2}}$
$c = \frac{1}{\sqrt{- \lambda _3}}$
And then the rotation matrix: $R = R_z(\alpha) \, R_y(\beta) \, R_x(\gamma) = \begin{bmatrix} \cos\alpha\cos\beta & \cos\alpha\sin\beta\sin\gamma - \sin\alpha\cos\gamma & \cos\alpha\sin\beta\cos\gamma + \sin\alpha\sin\gamma \\\\ \sin\alpha\cos\beta & \sin\alpha\sin\beta\sin\gamma + \cos\alpha\cos\gamma & \sin\alpha\sin\beta\cos\gamma - \cos\alpha\sin\gamma \\\\ -\sin\beta & \cos\beta\sin\gamma & \cos\beta\cos\gamma \end{bmatrix} $
is made from the eigenvectors of M, using the same $M=R D R^T$
So then:
$β = - \sin^{-1}(R _{20})$
Then with: $\frac{R _{21}}{R _{22}} = tan(γ)$
$γ = tan^{-1}( \frac{R _{21}}{ R _{22}} ) = atan2(\frac{R _{21}}{cos(β)}, \frac{R _{22}}{cos(β)}) = atan2(\frac{R _{21}}{\sqrt{1-R _{20}^2}}, \frac{R _{22}}{\sqrt{1-R _{20}^2}})$
and then with: $\frac{R _{10}}{R _{00}} = tan(α)$
$α = tan^{-1}( \frac{R _{10}}{ R _{00}} ) = atan2(\frac{R _{10}}{cos(β)}, \frac{R _{00}}{cos(β)}) = atan2(\frac{R _{10}}{\sqrt{1-R _{20}^2}}, \frac{R _{00}}{\sqrt{1-R _{20}^2}})$
All the rotation matrix stuff taken from G. Slabaugh Paper on Computing Euler Angles
Actual Questions
- Are a,b,c really just the reciprocal of the eigenvalues of M?
- Does this work for non-spherical ellipsoids (think that's a hyper-conic section? If it was 2D they would be hyperbola or parabola)
- From the previously quoted answer it looks like @Jan-Magnus Økland and @Amazing Dell Computers are proposing somewhat different methods. Background is engineering, so feel like I'm way out of my depth of math-fu on this.
@Amazing Dell Computers did something with dividing by $(r_0^T A r_0 - c)$, however that seems like it would either result in dividing a 3x3 by an nxn matrix or getting a M' for every single data point?
- If the 3x3 matrix is symmetric, is there a general decomposition for the eigenvalues (as long as its not defective)? (IE: can be written as an equation and does not just involve computing numbers?)
- I'd prefer a result I could write in terms A,B,C,D,E,F
- Something that looks like the ellipse results.
- Just from the center calculation, there's several optimizations obvious. Hoping other terms might simplify.
- Especially since decomposition directly results in rotation matrices made of cos() sin(), seems like there ought to be some relation for $M=R D R^T$
- I'd prefer a result I could write in terms A,B,C,D,E,F
- Are there any really obvious errors I've missed, or bad calculation assumptions?
Edit 12/29/2022
After a bit of trial and tribulations, appear to have a functioning result.
- Notably, dividing through by $(r_0^T A r_0 - c)$ did end up necessary, as otherwise the center offset is not taken into account correctly.
- Also, a, b, c appears to require Abs( eigenvalue )
$a = \frac{1}{\sqrt{\lvert \lambda _1 \rvert}}$
$b = \frac{1}{\sqrt{\lvert \lambda _2 \rvert}}$
$c = \frac{1}{\sqrt{\lvert \lambda _3 \rvert}}$
Otherwise, results like shown below, with obvious best fit spheres result in mixed sign eigenvalues.
Also, the original result I was trying to get, best fit sphere to stars within 100 LY of Earth
The equation for an ellipsoid in canonical form is given by:
$$\frac{(x-xc)^2}{a^2}+\frac{(y-yc)^2}{b^2}+\frac{(z-zc)^2}{c^2}=1$$
where $(xc,yc,zc)$ is the center of the ellipsoid, $a$, $b$, and $c$ are the lengths of the semi-axes, and the angles $\alpha$, $\beta$, and $\gamma$ determine the orientation of the ellipsoid. The angles are defined by:
$$\alpha=\arctan\left(\frac{D}{E}\right)$$ $$\beta=\arctan\left(\frac{F}{\sqrt{E^2+D^2}}\right)$$ $$\gamma=\arctan\left(\frac{E}{D}\right)$$
where the constants $A$, $B$, $C$, $D$, $E$, $F$, $G$, $H$, $I$, and $J$ are defined as follows:
$$A=\frac{a_0}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$ $$B=\frac{a_1}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$ $$C=\frac{a_2}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$ $$D=\frac{a_3}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$
$$E=\frac{a_4}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$ $$F=\frac{a_5}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$ $$G=\frac{a_6}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$ $$H=\frac{a_7}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$ $$I=\frac{a_8}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$ $$J=\frac{1}{(a_3/2)^2+(a_4/2)^2+(a_5/2)^2}$$
The center of the ellipsoid $(xc,yc,zc)$ is given by:
$$xc=-\frac{G}{2A}$$ $$yc=-\frac{H}{2B}$$ $$zc=-\frac{I}{2C}$$
The lengths of the semi-axes are given by:
$$a=\sqrt{-\frac{AH^2+BG^2-2GHD}{4AB-D^2}}$$
$$b=\sqrt{-\frac{AH^2+BG^2-2GHD}{4AB-D^2}}$$
$$c=\sqrt{-\frac{AI^2+CG^2-2GIH}{4AC-E^2}}$$
Note that these formulas assume that $A$, $B$, and $C$ are all positive. If any of these coefficients are negative, the corresponding axis length and angle will need to be negated.