I would like to simulate ellipsoid fitting.
In the first step I had ellipsoid with centre in 0,0,0 with specific length of axes a, b, c described by eq. $\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1$ and several vectors $v = (v_1, v_2, v_3)$. For these vectors I calculated coordinates x, y, z of intersection of ellipsoid and a straight line with a direction vector $v$ as $x = 0 + v_1t, y = 0 + v_2t, z = 0 + v_3t$ and substituted into the equation of ellipsoid.
By this I created several points on ellipsoid surface in specific directions and I could fit ellipsoid based on these points. I used two algorithms A, B and both worked well I mean eigenvalues were the same as a, b, c and eigenvectors were on x, y, z axes resp.
But I had problem when I wanted to simulate the same ellipsoid but rotated by some angel around some axes, for example beta around y axe. Based on this text I used rotation matrix $R_y (\beta)$ so my new coordinates were $x^{´} = xcos\beta + zsin\beta$, $y^{´}=y$, $z^{´}=zcos\beta - xsin\beta$ so I substituted old coordinates in ellipsoid eq. by new ones and by the same way I created dataset of surface points in specific directions.
But when I fitted the rotated ellipsoid the eigenvalues were not a, b, c and eigenvectors were not in y axe and in rotated x and z axes by $\beta$.
And now I do not know if I made some mistake in creating/calculating my dataset or the fitting algorithms did not work.
Can you help me in this problem? Thank you.
Mark
Yes. You must have made an error somewhere. To simulate identifying a rotated ellipsoid, you must first create points on the unrotated ellipsoid, then rotate and shift them under the rotation matrix $R$ and a shift vector $r_0$, thus generating a new set of points.
Now, the original (basic) ellipsoid had the equation
$$ r^T Q_0 r = 1 $$
The image of coordinate vector $r$ is $ r' = R r + r_0 $. It follows that the equation of the rotated/shifted ellipsoid is
$$ (r' - r_0)^T R Q_0 R^T (r' - r_0) = 1 $$
Let set $Q = R Q_0 R^T $ then the equation is
$$ (r - r_0)^T Q (r - r_0 ) = 1 $$
Note the $3 \times 3$ matrix $Q$ is of the form
$$ Q = \begin{bmatrix} a && d && e \\ d && b && f \\ e && f && c \end{bmatrix} $$
Thus, by multiplying out the above quadratic form, we obtain
$ a (x - x_0)^2 + b (y - y_0)^2 + c (z - z_0)^2 \\+ 2 d (x - x_0)(y - y_0) + 2 e (x - x_0)(z - z_0) + 2 f (y - y_0) (z - z_0) = 1 $
Thus the equation of the ellipsoid is of the form
$$ A x^2 + B y^2 + C z^2 + D x y + E x z + F y z + G x + H y + I z + J = 0 $$
All the points generated by the simulator will satisfy this equation, thus we can write
$$ A x_i^2 + B y_i^2 + C z_i^2 + D x_i y_i + E x_i z_i + F y_i z_i + G x_i + H y_i + I z_i + J = 0 $$
The unknowns here are the constants $A$ through $J$. So we'll define a vector $\theta$ as follows
$ \theta = [ A, B, C, D, E, F, G, H, I, J ]^T $
and we'll define the matrix $\Phi$ to be the $N \times 10$ matrix whose $i$-th row is
$ \Phi_i = \begin{bmatrix}x_i^2 && y_i^2 && z_i^2 && x_i y_i && x_i z_i && y_i z_i && x_i && y_i && z_i && 1 \end{bmatrix} $
Now our equations can be written as
$ \Phi \theta = 0 $
This is an $N$ dimensional vector. To be able to solve for $\theta$, we'll pre-multiply both sides of the equation by $ \Phi^T$, to get
$ \Phi^T \Phi \theta = 0 $
Which is a set of $10$ equation in $10$ unknowns. To get a non-zero solution, we'll truncate the system of equations, and take the first $9$ equations, thus disregarding the last equation. Now, using Gauss-Jordan elimination, we can solve the $ 9 \times 10$ system of equations for the the direction of vector $\theta$.
The solution will be of the form $ \theta = t \varphi $
where $\varphi$ is a fixed vector and $ t $ is an arbitrary scalar.
From here, we now have matrix $Q$ correct up to a scaling factor, as follows
$ Q = t \begin{bmatrix} \varphi_1 && \dfrac{1}{2} \varphi_4 && \dfrac{1}{2} \varphi_5 \\ \dfrac{1}{2} \varphi_4 && \varphi_2 && \dfrac{1}{2} \varphi_6 \\ \dfrac{1}{2} \varphi_5 && \dfrac{1}{2} \varphi_6 && \varphi_3 \end{bmatrix} $
But first, from our equation of the rotated/shifted ellipsoid, we have
$$ (r - r_0)^T Q (r - r_0) = 1 $$
which when expanded becomes
$$ r^T Q r + r^T a + b = 0 $$
where $ a = - 2 Q r_0 $ and $ b = r_0^T Q r_0 - 1 $
It follows that $r_0 = -\dfrac{1}{2} Q^{-1} a $
Now from our identification we have
$ a = t \begin{bmatrix} \varphi_7 && \varphi_8 && \varphi \end{bmatrix} $
Therefore, we can determine the shift vector $r_0$ because the $t$'s cancel out.
Once this is done, we can determine the value of $t$ from
$ b = r_0^T Q r_0 - 1 $
which is linear in $t$. Then our matrix $Q$ is fully specified.
But we're not done yet, to find the semi-axes of the ellipsoid we have to diagonalize $Q$ into $ Q = R D R^T $ and this can done using, for example, MATLAB, or Mathematica. The diagonal matrix $D$ has on its diagonal the reciprocal of the semi-axes lengths squared, i.e.
$ D = \text{diag}(\dfrac{1}{a^2}, \dfrac{1}{b^2}, \dfrac{1}{c^2} ) $
Hence, using this reverse-engineering process (i.e. identification) we were able to retrieve all the parameters of the original ellipsoid, its rotation matrix, and it shift vector.