I'm trying to compute the angle and major/minor axes of the perspective projection of a circle. The circle is embedded in $\mathbb R^3$ space with any orientation not perpendicular to the screen, but I use the plane $z=1$ for screen and $(0,0,0)$ as my focal point.
I manage to get the equation of the conic in general form, but in the case it is an ellipse, the computation to get the center axis-aligned form gets out of hand and I'm feeling I'm missing something to compute it.
So here is what I've done. Let be a circle of center $c = (c_x, c_y, c_z)$ and radius $r$, with a normal $\mathbf n = (n_x,n_y,n_z)$. It is described by the following set of equations, one equation of a plane and one of a sphere, \begin{gather} \begin{cases} (x-c_x)^2 + (y-c_y)^2 + (z-c_z)^2 = r^2, \\ n_x (x - c_x) + n_y (x - c_y) + n_z (x - c_z) = 0.\\ \end{cases} \end{gather}
They reduce to
\begin{gather} \newcommand{\nc}{{\langle\mathbf n|c\rangle}} \newcommand{\nxy}{{\langle\mathbf n|(x,y,1)\rangle}} \begin{cases} x^2 - 2c_xx + c_x^2 + y^2 - 2c_yy + c_y^2 + z^2 - 2c_zz + c_z^2 - r^2 = 0, \\ n_xx + n_yy + n_zz = \nc, \\ \end{cases}\\ \end{gather}
where $\langle\_|\_\rangle$ is the scalar product.
My perspective projection is the function
\begin{align*} f: & \mathbb R^3 \to \mathbb R^2 \times \{1\} \\ & (x,y,z) \mapsto \left( \frac xz, \frac yz \right)\\ \end{align*}
And the inverse projection to the plane containing the circle is \begin{align*} f^{-1}: \mathbb R^2 \times \{1\} \to& \mathbb{R}^3 \\ (x,y,1) \mapsto& \left( \frac{x\nc}{\nxy}, \frac{y\nc}{\nxy}, \frac{\nc}{\nxy} \right)\\ \end{align*}
If $(x,y)$ is part of the projected circle, then the inverse projection of that point must also be on the sphere that contained the circle.
\begin{gather*} (f^{-1}(x)-c_x)^2 + (f^{-1}(y)-c_y)^2 + (f^{-1}(z)-c_z)^2 = r^2 \\ \end{gather*}
Expanding that, we have \begin{gather*} f^{-1}(x)^2 - 2c_xf^{-1}(x) + c_x^2 + f^{-1}(y)^2 - 2c_yf^{-1}(y) + c_y^2 + f^{-1}(z)^2 - 2c_zf^{-1}(z) + c_z^2 - r^2 = 0 \\ \end{gather*} \begin{align} \left(\frac{x\nc}{\nxy}\right)^2 - 2c_x\frac{x\nc}{\nxy} + c_x^2 + \left(\frac{y\nc}{\nxy}\right)^2 - 2c_y\frac{y\nc}{\nxy} + c_y^2 &\\ + \left(\frac{\nc}{\nxy}\right)^2 - 2c_z\frac{\nc}{\nxy} + c_z^2 - r^2 &= 0 \\ \end{align}
Still expanding, this already gets quite verbose
\begin{align*} &\left(x\nc\right)^2 &&- 2c_x x\nc \nxy &&+ c_x^2 \nxy^2 \\ +& \left(y\nc\right)^2 &&- 2c_y y\nc \nxy &&+ c_y^2 \nxy^2 \\ +& \left(\nc\right)^2 &&- 2c_z \nc \nxy &&+ c_z^2 \nxy^2 -r^2 \nxy^2 = 0 \end{align*}
\begin{align*} & \nc^2 x^2 &-& 2 c_x \nc n_x x^2 &-& 2 c_x \nc n_y xy &-& 2 c_x \nc n_z x &\\ +& c_x^2 n_x^2 x^2 &+& c_x^2 n_y^2 y^2 &+& c_x^2 n_z^2 &+& 2 c_x^2 n_x n_y xy &+& 2 c_x^2 n_x n_z x &+& 2 c_x^2 n_y n_z y &\\ +& \nc^2 y^2 &-& 2 c_y \nc n_x xy &-& 2 c_y \nc n_y y^2 &-& 2 c_y \nc n_z y &\\ +& c_y^2 n_x^2 x^2 &+& c_y^2 n_y^2 y^2 &+& c_y^2 n_z^2 &+& 2 c_y^2 n_x n_y xy &+& 2 c_y^2 n_x n_z x &+& 2 c_y^2 n_y n_z y &\\ +& \nc^2 &-& 2 c_z \nc n_x x &-& 2 c_z \nc n_y y &-& 2 c_z \nc n_z &\\ +& c_z^2 n_x^2 x^2 &+& c_z^2 n_y^2 y^2 &+& c_z^2 n_z^2 &+& 2 c_z^2 n_x n_y xy &+& 2 c_z^2 n_x n_z x &+& 2 c_z^2 n_y n_z y &\\ -& r^2 n_x^2 x^2 &-& r^2 n_y^2 y^2 &-& r^2 n_z^2 &-& 2 r^2 n_x n_y xy &-& 2 r^2 n_x n_z x &-& 2 r^2 n_y n_z y &=& \quad 0 \\ \end{align*}
Now, writing the quadratic equation by component in the form $A\,xy + B\, x^2 + C\,y^2 + D\,x + E\,y + F$, the factors are as below.
\begin{align*} x y &:A=&& &-& &2& (c_xn_y + c_yn_x) \nc &+& &2& \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y \\ x^2 &:B=&& \nc^2 &-& &2& c_x n_x \nc &+& & & \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) n_x^2 \\ y^2 &:C=&& \nc^2 &-& &2& c_y n_y \nc &+& & & \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) n_y^2 \\ x &:D=&& &-& &2& (c_x n_z + c_z n_x) \nc &+& &2& \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) n_x n_z \\ y &:E=&& &-& &2& (c_y n_z + c_z n_y) \nc &+& &2& \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) n_y n_z \\ 1 &:F=&& \nc^2 &-& &2& c_z n_z \nc &+& & & \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) n_z^2 \\ \end{align*}
According to wikipedia, this quadratic equation is an ellipse if $4BC - A^2 > 0$.
I'll spare you the details, but I check that condition and arrived at the conclusion that
\begin{gather*} r < \frac{ \sqrt{n_x^2 + n_y^2 + n_z^2} }{\sqrt{n_x^2 + n_y^2} } |c_z| \end{gather*}
must be true. It seems reasonable as it means that the circle must not cross the plane $z=0$.
However, now I want to compute the anngle of the projected ellipse relative to the x-axis. Again, according to wikipedia, this can be computed with the following formula : \begin{gather*} \theta = \begin{cases} 0 \text{ if } A = 0, B < C \\ \frac\pi2 \text{ if } A = 0, B > C \\ \arctan\frac{C - B - \sqrt{(C-B)^2 + A^2}}{A} \text{ if } A \ne 0 \end{cases} \end{gather*}
Let's assume for now that $A \ne 0$, and let's compute $\frac{C - B - \sqrt{(C-B)^2 + A^2}}{A}$.
This is where this is getting ugly. So, first I replace the $A$, $B$ and $C$ values. \begin{gather*} \frac{ \left( \nc^2 - 2 c_y n_y \nc + \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) n_y^2 \right) - \left( \nc^2 - 2 c_x n_x \nc + \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) n_x^2 \right) } { - 2 (c_xn_y + c_yn_x) \nc + 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y }\\ +\frac{ -\sqrt{ \begin{gathered} \left( \left( \nc^2 - 2 c_y n_y \nc + \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) n_y^2 \right) - \left( \nc^2 - 2 c_x n_x \nc + \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) n_x^2 \right) \right)^2\\ + \left( - 2 (c_xn_y + c_yn_x) \nc + 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y \right) ^2 \end{gathered} } } { - 2 (c_xn_y + c_yn_x) \nc + 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y } \end{gather*}
The two first lines are one giant numerator as are the three lines that are part of the square root. Now let's compute more (I wont't write the equality symbols).
\begin{gather*} \frac{ \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) (n_x^2 - n_y^2) + 2 ( c_x n_x - c_y n_y ) \nc } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) \nc }\\ +\frac{ -\sqrt{ \begin{gathered} \Big( \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) (n_x^2 - n_y^2) + 2 ( c_x n_x - c_y n_y ) \nc \Big)^2\\ + \Big( 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) \nc \Big) ^2 \end{gathered} } } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) \nc } \end{gather*}
Then, \begin{gather*} \frac{ \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) (n_x^2 - n_y^2) + 2 ( c_x n_x - c_y n_y ) \nc } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) \nc }\\ +\frac{ -\sqrt{ \begin{gathered} \left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 (n_x^2 - n_y^2)^2 + 4 \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) (n_x^2 - n_y^2) ( c_x n_x - c_y n_y ) \nc\\ + 4 ( c_x n_x - c_y n_y )^2 \nc^2 \\ + 4 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right)^2 n_x^2 n_y^2 - 8 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y (c_xn_y + c_yn_x) \nc \\ + 4 (c_xn_y + c_yn_x)^2 \nc^2 \end{gathered} } } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) \nc } \end{gather*}
Then, \begin{gather*} \frac{ \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) (n_x^2 - n_y^2) + 2 ( c_x n_x - c_y n_y ) \nc } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) \nc }\\ +\frac{ -\sqrt{ \begin{gathered} \left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_x^4 -2 \left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_x^2 n_y^2 +\left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_y^4 \\ + 4 \nc \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_x n_x^3 + 4 \nc \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_y n_x^2 n_y \\ + 4 \nc \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_x n_x n_y^2 + 4 \nc \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_y n_y^3 \\ + 4 \nc^2 c_x^2 n_x^2 + 8 \nc^2 c_x c_y n_x n_y + 4 \nc^2 c_y^2 n_y^2 \\ + 4 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right)^2 n_x^2 n_y^2 \\ - 8 \nc \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) c_x n_x n_y^2 - 8 \nc \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) c_y n_x^2 n_y \\ + 4 \nc^2 c_x^2 n_y^2 + 8 \nc^2 c_x c_y n_x n_y + 4 \nc^2 c_y^2 n_x^2 \end{gathered} } } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) \nc } \end{gather*}
Then, \begin{gather*} \frac{ \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) (n_x^2 - n_y^2) + 2 ( c_x n_x - c_y n_y ) \nc } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) \nc }\\ +\frac{ -\sqrt{ \begin{gathered} \left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_x^4 +2 \left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_x^2 n_y^2 +\left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_y^4 \\ + 4 \nc \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_x n_x^3 - 4 \nc \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_y n_x^2 n_y \\ - 4 \nc \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_x n_x n_y^2 + 4 \nc \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_y n_y^3 \\ + 4 \nc^2 c_x^2 n_x^2 + 4 \nc^2 c_x^2 n_y^2 + 4 \nc^2 c_y^2 n_x^2 + 4 \nc^2 c_y^2 n_y^2 + 16\nc^2 c_x c_y n_x n_y \end{gathered} } } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) \nc } \end{gather*}
Then, \begin{gather*} \frac{ \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) (n_x^2 - n_y^2) + 2 ( c_x n_x - c_y n_y ) (n_xc_x+n_yc_y+n_zc_z) } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) (n_xc_x+n_yc_y+n_zc_z) }\\ +\frac{ -\sqrt{ \begin{gathered} \left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_x^4 +2 \left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_x^2 n_y^2 +\left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_y^4 \\ + 4 (n_xc_x+n_yc_y+n_zc_z) \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_x n_x^3 - 4 (n_xc_x+n_yc_y+n_zc_z) \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_y n_x^2 n_y \\ - 4 (n_xc_x+n_yc_y+n_zc_z) \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_x n_x n_y^2 + 4 (n_xc_x+n_yc_y+n_zc_z) \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_y n_y^3 \\ + 4 (n_xc_x+n_yc_y+n_zc_z)^2 c_x^2 n_x^2 + 4 (n_xc_x+n_yc_y+n_zc_z)^2 c_x^2 n_y^2 + 4 (n_xc_x+n_yc_y+n_zc_z)^2 c_y^2 n_x^2 \\ + 4 (n_xc_x+n_yc_y+n_zc_z)^2 c_y^2 n_y^2 + 16(n_xc_x+n_yc_y+n_zc_z)^2 c_x c_y n_x n_y \end{gathered} } } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) (n_xc_x+n_yc_y+n_zc_z) } \end{gather*}
And finally I reach this monster : \begin{gather*} \frac{ \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) (n_x^2 - n_y^2) + 2 ( c_x n_x - c_y n_y ) (n_xc_x+n_yc_y+n_zc_z) } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) (n_xc_x+n_yc_y+n_zc_z) }\\ +\frac{ -\sqrt{ \begin{gathered} c_x^4 n_x^4 + c_y^4 n_x^4 + c_z^4 n_x^4 + n_x^4 r^4 + 2 c_x^2 c_y^2 n_x^4 + 2 c_x^2 c_z^2 n_x^4 + 2 c_y^2 c_z^2 n_x^4 - 2 c_x^2 r^2 n_x^4 - 2 c_y^2 r^2 n_x^4 - 2 c_z^2 r^2 n_x^4 %\left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_x^4 \\ + 2c_x^4 n_x^2 n_y^2 + 2c_y^4 n_x^2 n_y^2 + 2c_z^4 n_x^2 n_y^2 + 2n_x^2 n_y^2 r^4 + 4 c_x^2 c_y^2 n_x^2 n_y^2 + 4 c_x^2 c_z^2 n_x^2 n_y^2 + 4 c_y^2 c_z^2 n_x^2 n_y^2 - 4 c_x^2 r^2 n_x^2 n_y^2 \hfill\\\hfill - 4 c_y^2 r^2 n_x^2 n_y^2 - 4 c_z^2 r^2 n_x^2 n_y^2 %+2 \left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_x^2 n_y^2 \\ + c_x^4 n_y^4 + c_y^4 n_y^4 + c_z^4 n_y^4 + n_y^4 r^4 + 2 c_x^2 c_y^2 n_y^4 + 2 c_x^2 c_z^2 n_y^4 + 2 c_y^2 c_z^2 n_y^4 - 2 c_x^2 r^2 n_y^4 - 2 c_y^2 r^2 n_y^4 - 2 c_z^2 r^2 n_y^4 %+\left(c_x^2 + c_y^2 + c_z^2 - r^2\right)^2 n_y^4 \\ + 4 c_x^4 n_x^4 + 4 c_x^2 c_y^2 n_x^4 + 4 c_x^2 c_z^2 n_x^4 - 4 c_x^2 r^2 n_x^4 + 4 c_x^3 c_y n_x^3 n_y + 4 c_x c_y^3 n_x^3 n_y + 4 c_x c_y c_z^2 n_x^3 n_y \hfill\\\hfill - 4 c_x c_y r^2 n_x^3 n_y + 4 c_x^3 c_z n_x^3 n_z + 4 c_x c_y^2 c_z n_x^3 n_z + 4 c_x c_z^3 n_x^3 n_z - 4 c_x c_z r^2 n_x^3 n_z %- 4 (n_xc_x+n_yc_y+n_zc_z) \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_x n_x^3 \\ - 4 c_x^3 c_y n_x^3 n_y - 4 c_x c_y^3 n_x^3 n_y - 4 c_x c_y c_z^2 n_x^3 n_y + 4 c_x c_y r^2 n_x^3 n_y - 4 c_x^2 c_y^2 n_x^2 n_y^2 - 4 c_y^4 n_x^2 n_y^2 - 4 c_y^2 c_z^2 n_x^2 n_y^2 \hfill\\\hfill + 4 c_y^2 r^2 n_x^2 n_y^2 - 4 c_x^2 c_y c_z n_x^2 n_y n_z - 4 c_y^3 c_z n_x^2 n_y n_z - 4 c_y c_z^3 n_x^2 n_y n_z + 4 c_y c_z r^2 n_x^2 n_y n_z %- 4 (n_xc_x+n_yc_y+n_zc_z) \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_y n_x^2 n_y \\ - 4 c_x^4 n_x^2 n_y^2 - 4 c_x^2 c_y^2 n_x^2 n_y^2 - 4 c_x^2 c_z^2 n_x^2 n_y^2 + 4 c_x^2 r^2 n_x^2 n_y^2 - 4 c_x^3 c_y n_x n_y^3 - 4 c_x c_y^3 n_x n_y^3 - 4 c_x c_y c_z^2 n_x n_y^3 \hfill\\\hfill + 4 c_x c_y r^2 n_x n_y^3 - 4 c_x^3 c_z n_x n_y^2 n_z - 4 c_x c_y^2 c_z n_x n_y^2 n_z - 4 c_x c_z^3 n_x n_y^2 n_z + 4 c_x c_z r^2 n_x n_y^2 n_z %- 4 (n_xc_x+n_yc_y+n_zc_z) \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_x n_x n_y^2 \\ + 4 c_x^3 c_y n_x n_y^3 + 4 c_x c_y^3 n_x n_y^3 + 4 c_x c_y c_z^2 n_x n_y^3 - 4 c_x c_y r^2 n_x n_y^3 + 4 c_x^2 c_y^2 n_y^4 + 4 c_y^4 n_y^4 + 4 c_y^2 c_z^2 n_y^4 \hfill\\\hfill - 4 c_y^2 r^2 n_y^4 + 4 c_x^2 c_y c_z n_y^3 n_z + 4 c_y^3 c_z n_y^3 n_z + 4 c_y c_z^3 n_y^3 n_z - 4 c_y c_z r^2 n_y^3 n_z %+ 4 (n_xc_x+n_yc_y+n_zc_z) \left(c_x^2 + c_y^2 + c_z^2 - r^2\right) c_y n_y^3 \\ + 4 c_x^4 n_x^4 + 4 c_x^2 c_y^2 n_x^2 n_y^2 + 4 c_x^2 c_z^2 n_x^2 n_z^2 + 8 c_x^3 c_y n_x^3 n_y + 8 c_x^3 c_z n_x^3 n_z + 8 c_x^2 c_y c_z n_x^2 n_y n_z %+ 4 (n_xc_x+n_yc_y+n_zc_z)^2 c_x^2 n_x^2 \\ + 4 c_x^4 n_x^2 n_y^2 + 4 c_x^2 c_y^2 n_y^4 + 4 c_x^2 c_z^2 n_y^2 n_z^2 + 8 c_x^3 c_y n_x n_y^3 + 8 c_x^3 c_z n_x n_y^2 n_z + 8 c_x^2 c_y c_z n_y^3 n_z %+ 4 (n_xc_x+n_yc_y+n_zc_z)^2 c_x^2 n_y^2 \\ + 4 c_x^2 c_y^2 n_x^4 + 4 c_y^4 n_x^2 n_y^2 + 4 c_y^2 c_z^2 n_x^2 n_z^2 + 8 c_x c_y^3 n_x^3 n_y + 8 c_x c_y^2 c_z n_x^3 n_z + 8 c_y^3 c_z n_x^2 n_y n_z %+ 4 (n_xc_x+n_yc_y+n_zc_z)^2 c_y^2 n_x^2 \\ + 4 c_x^2 c_y^2 n_x^2 n_y^2 + 4 c_y^4 n_y^4 + 4 c_y^2 c_z^2 n_y^2 n_z^2 + 8 c_x c_y^3 n_x n_y^3 + 8 c_x c_y^2 c_z n_x n_y^2 n_z + 8 c_y^3 c_z n_y^3 n_z %+ 4 (n_xc_x+n_yc_y+n_zc_z)^2 c_y^2 n_y^2 \\ + 16 c_x^3 c_y n_x^3 n_y + 16 c_x c_y^3 n_x n_y^3 + 16 c_x c_y c_z^2 n_x n_y n_z^2 + 32 c_x^2 c_y^2 n_x^2 n_y^2 + 32 c_x^2 n_y c_z n_x^2 n_y n_z + 32 c_x c_y^2 c_z n_x n_y^2 n_z %+ 16(n_xc_x+n_yc_y+n_zc_z)^2 c_x c_y n_x n_y \end{gathered} } } { 2 \left(c_x^2 + c_y^2 + c_z ^2 - r^2\right) n_x n_y - 2 (c_xn_y + c_yn_x) (n_xc_x+n_yc_y+n_zc_z) } \end{gather*}
At which point I give up because I don't know an efficient way to simplify that. So, my question is, do you know an alternate strategy to compute that angle $\theta$ ? Or do you know a way to simplify that ugly square root or a proper way to approximate it ?
Thank you for any help.
Originally I had suggested that starting from the general form you found, you can extract center, orientation and radii of the ellipse following e.g. the steps I wrote in this answer. But I guess that approach misses the point of finding a closed formula for the angle. The comment by Jan-Magnus Økland already suggests an approach for that, namely going from $\tan\theta$ to $\tan2\theta$. To this effect read Wikipedia to find
$$\tan2\theta=\frac{2\tan\theta}{1-\tan^2\theta}$$
Where does the $\theta=\arctan\frac{C - B - \sqrt{(C-B)^2 + A^2}}{A}$ come from? This looks like the solution to some quadratic equation, using the well-known quadratic formula. A bit of experimenting and you find a corresponding quadratic equation
\begin{align*} \frac12A\tan^2\theta+(B-C)\tan\theta-\frac12A&=0\\ 2\tan\theta(B-C)&=\left(1-\tan^2\theta\right)A\\ \tan2\theta=\frac{2\tan\theta}{1-\tan^2\theta}&=\frac{A}{B-C}\\ \theta&=\frac12\arctan\frac{A}{B-C}+k\frac\pi2 \end{align*}
for some $k\in\mathbb Z$ which gives the angle up to some multiple of $\frac\pi2$ as is to be expected. Picking the right one among these algebraically equivalent ones is hard to do with a closed formula.
But I'd also suggest you revisit how you obtain that general form in the first place. Here is a higher level view which avoids much of the tedious coordinate computation.
A point $p$ is on your circle if $\langle p-c|p-c\rangle=r^2$ and $\langle p-c|n\rangle=0$. You already have that much, although not using this notation. Rewrite these conditions to
\begin{align*} \langle p|p\rangle-2\langle p|c\rangle+\langle c|c\rangle-r^2&=0\\ \langle p|n\rangle -\langle c|n\rangle&=0 \end{align*}
The first equation now has a quadratic term, a linear term and a constant term. The second equation has only a linear and a constant term. If your final projected point on the plane is $q=(x,y,1)$, then the corresponding preimage is some multiple thereof, i.e. some $p=tq$ for $t\in\mathbb R$. You can plug this into the equation of the plane and solve for $t$:
\begin{align*} \langle tq|n\rangle-\langle c|n\rangle&=0\\ t\langle q|n\rangle&=\langle c|n\rangle\\ t&=\frac{\langle c|n\rangle}{\langle q|n\rangle} \end{align*}
Then use this in the equation of the sphere:
\begin{align*} t^2\langle q|q\rangle-2t\langle q|c\rangle+\langle c|c\rangle-r^2&=0\\ \langle c|n\rangle^2\langle q|q\rangle -2\langle c|n\rangle\langle q|n\rangle\langle q|c\rangle +\langle q|n\rangle^2\left(\langle c|c\rangle-r^2\right)&=0 \end{align*}
Note how I avoided division by multiplying the whole equation by the square of the denominator of $t$. If you are deep into bra-ket notation, you could also write all of this as
$$ \langle q|\,\left( \langle c|n\rangle^2\cdot\mathbf 1 -\langle c|n\rangle\cdot|c\rangle\langle n| -\langle n|c\rangle\cdot|n\rangle\langle c| +\left(\langle c|c\rangle-r^2\right)\cdot|n\rangle\langle n| \right)\,|q\rangle=0 $$
where $|\_\rangle\langle\_|$ denotes the outer product of two vectors, i.e. a $3\times3$ matrix. So the big parenthesis contains a sum of four terms. Each of these terms is scalar times matrix. The first matrix is the identity matrix $\mathbf 1$ which is symmetric. The second and third term are mutually transpose, so their sum is symmetric. (Writing $\langle c|n\rangle$ as the scalar factor for one term and $\langle n|c\rangle$ for the other is merely to highlight the symmetry, they are of course the same number.) The last matrix uses twice the same vector so it is symmetric, too. Thus the whole big parenthesis describes a symmetric matrix. Using a bit of a more conventional notation and assuming column vectors, you could write this as
$$q^T\,\left( (c^Tn)^2\,\mathbf 1 -(c^Tn)(cn^T) -(n^Tc)(nc^T) +(c^Tc-r^2)(nn^T) \right)\,q=0$$
The matrix computed by the big parenthesis is the homogeneous matrix representation of your conic section. When you wrote “I manage to get the equation of the conic in general form”, this is what I had expected at first. In the notation you use, this matrix corresponds to
$$\begin{pmatrix}B&A/2&D/2\\A/2&C&E/2\\D/2&E/2&F\end{pmatrix}$$
I do hope that the above will give you some idea of how to obtain these parameters in a more high-level view, with less coordinates cluttering the path. The $3\times3$ matrix notation can be useful for other parameters as well. For example the center of the conic has homogeneous coordinates equal to the cross product of the first two columns of that matrix.