You need $9$ points. The generic equation of an ellipsoid centered at $(x_0, y_0, z_0)$ is of the form \[ \begin{bmatrix} x-x_0 & y-y_0 & z-z_0 \end{bmatrix} \begin{bmatrix} a & d & e \\ d & b & f \\ e & f & c \end{bmatrix} \begin{bmatrix} x-x_0 \\ y-y_0 \\ z-z_0 \end{bmatrix} = 1, \] where the central matrix is positive definite.
Assume for the moment $(x_0, y_0, z_0) = (0,0,0)$. The equation translates to the equation \[ (ax+dy+ez)x + (dx+by+fz)y + (ex+fy+cz)z = 1 \] which we can rewrite as \[ (\bullet) \quad ax^2 + by^2 + cz^2 + 2dxy + 2exz + 2fyz = 1. \] To account for the unknown position of the center we substitute back $x \mapsto x-x_0$ etc, and get \[ a(x-x_0)^2 + b(y-y_0)^2 + c(z-z_0)^2 + 2d(x-x_0)(y-y_0) + 2e(x-x_0)(z-z_0) + 2f(y-y_0)(z-z_0) = 1 \] which expands to \[ ax^2 + by^2 + cz^2 + 2dxy + 2exz + 2fyz + rx + sy + tz + u = 1 \] for some values of $r,s,t,u$ (explicitly computed later - we don't need them for now). Notice that the coefficients of the degree $2$ terms are the same as before. Call $a' = \frac{a}{1-u}$ and the same for the other variables. We can rewrite the equation as \[ (\ast) \quad a' x^2 + b' y^2 + c' z^2 + 2d' xy + 2e' xz + 2f' yz + r'x + s'y + t'z = 1. \]
Now, for any choice of a point $x_i, y_i, z_i$ on the ellipsoid, we can substitute the entries in $(\ast)$ and get a linear equation in $9$ variables. If you choose the points appropriately, so that the corresponding linear system is full rank, you get that the system has a unique solution, which you can compute using any online linear system solver (e.g. https://matrixcalc.org/en/slu.html is an option).
To get the original values, we need to compute $1-u$. We have that \[ u = ax_0^2 + by_0^2 + cz_0^2 + 2dx_0y_0 + 2ex_0z_0 + 2fy_0 z_0 \] and that $r = -2ax_0 -2dy_0 -2ez_0$, $s = -2dx_0 -2by_0 -2fz_0$, $t = -2ex_0 -2fy_0 -2cz_0$. We know the values of $r', s', t'$ from the system above, so we can write the equations $r = r'(1-u)$ (and the same for $s$ and $t$), get a system of three quadratic equations in three variables ($x_0, y_0, z_0$). Again, this can be solved by computer (e.g. https://www.hackmath.net/en/calculator/solving-system-of-linear-equations this should work). This gives you the coordinates of the center.
Once you manage to compute $u$, from the previous solution you can get the values of $a,b,c,d,e,f$. Now, compute the eigenvalues of the matrix given at the beginning. If the eigenvalues are some $u, v, w$, (which are going to be positive because the matrix is positive definite) then the lengths of the semiaxes are $\frac{1}{\sqrt{u}}$, $\frac{1}{\sqrt{v}}$, $\frac{1}{\sqrt{w}}$.
Any positive definite matrix $A$ can be factored as $Q D Q^T$, where $Q$ is an orthonormal matrix (determining a rotation) and $D$ is a diagonal matrix whose entries are $u,v,w$. This is called singular value decomposition and can be computed by computer again (e.g. https://www.emathhelp.net/calculators/linear-algebra/svd-calculator/ should work). The columns of $Q$ represent the (normalized) axes of the ellipsoid, from which you can easily compute the angles (you can e.g. read the entries of the first colum of $Q$ as $(\cos\theta \cos\phi, \cos\theta\sin\phi, \sin\theta)$ and get the angles from there).
Note that doing the computations explicitly is a massive amount of work, the formula unfortunately doesn't get any simpler than this. It would, however, massively simplify if you knew that the ellipsoid is centered at the origin. Then you can get the matrix directly from $(\bullet)$, only needing $6$ points, skipping the quadratic system of three equations, and doing only linear systems.
Poincare
Savionf
