Parameters of the conic curve resulting from the perspective projection of a circle

Question

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.

Answer 1

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.

Answer 2

I don't know if this answer will be of help to you, nonetheless I would like to show a more geometric approach.

In diagram below, $A=(0,0,0)$ is the projection point, the yellow plane is the screen ($z=1$) while the circle (of center $G$) lies on the blue plane. The two planes meet at a certain line $r$ (dashed red) and let $BC$ be that diameter of the circle such that line $BC$ is perpendicular to $r$. Moreover, let $F$ be the intersection of line $BC$ with plane $z=0$.

If projection lines $AB$ and $AC$ meet the screen at $M$ and $N$, $MN$ is a diameter of the projected ellipse and its midpoint $O$ is thus the center of the ellipse. Its conjugate diameter $ED$ lies on the line through $O$ parallel to $r$, and $$ OD=OE=OM{\sqrt{BF\cdot CF}\over AF}. $$

Once you have two conjugate diameters, you can find the principal axes as follows (see diagram below). In the plane of the ellipse, draw line $QQ'$ through $N$ perpendicular to $DE$. Points $Q$ and $Q'$ must be chosen such that $NQ=NQ'=OD$. Principal axes $IR$ and $ST$ lie then on the bisectors of the angles formed by lines $OQ$ and $OQ'$, and: $$ IR=OQ'+OQ,\quad TS=OQ'-OQ. $$

Answer 3

$\require{begingroup}\begingroup $To recapitulate what you have already found, $$ \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 \neq 0 \end{cases} $$ where \begin{align} \newcommand{\n}{{\mathbf n}} \newcommand{\braket}[2]{\langle #1 | #2 \rangle} A &= 2 n_x n_y \left(\braket cc - r^2\right) - 2 (c_xn_y + c_yn_x) \braket\n c, \\ B &= \braket\n c^2 - 2 c_x n_x \braket\n c + \left(\braket cc - r^2\right) n_x^2, \\ C &= \braket\n c^2 - 2 c_y n_y \braket\n c + \left(\braket cc - r^2\right) n_y^2. \\ \end{align}

It seems to me that this is a simplified formula for the angle you are trying to compute. (By "this" I mean everything after the word found and before the word it in the text above. I have nearly copied this directly from the question, except that I have used the bra-ket notation to express $c_x^2+c_y^2+c_z^2$ a little more succinctly and I have rearranged some terms.) Considering the number of inputs into the problem ($r$ and all the components of $c$ and $\n$), the formula is not even especially complicated. It only becomes unwieldy when you expand it.

Moreover, I would argue that this is a closed-form expression in the strictest sense, because it gives you a bounded number of steps in which to compute the value, every step consisting of an operation such as addition, subtraction, multiplication, division, the square root, and the arc tangent function (and we could get rid of the arc tangent if you are willing to accept coefficients of an equation of the major axis instead of requiring an angle in radians).

If you are willing to consider a separate formula for the $A\neq 0$ case, you can write $$ \tan\theta = \frac{H - \sqrt{H^2 + A^2}}{A} $$ where \begin{align} A &= 2 n_x n_y \left(\braket cc - r^2\right) - 2 (c_x n_y + c_y n_x) \braket\n c, \\ H &= \left(n_y^2 - n_x^2\right)\left(\braket cc - r^2\right) - 2 (c_y n_y - c_x n_x) \braket\n c. \end{align}

If you don't like the use of the word where and the way I go in and out of "display" mode within the blocks of text that I am calling closed-form expressions, you can adopt a more formal notation that will allow you to have the effect of the word where within one continuous symbolic expression. An example of such a notation is the LET construct in Leslie Lamport's TLA+ specification language.

---

These formulas are simple enough to start plugging in numbers to explore the behavior of the conic. For example, let's take $\n = (0.6,0,0.8)$ and $c = (0,3,4).$ Since $n_y = 0,$ these conditions give us $A = -11.52$ independent of $r.$ But $\left(n_y^2 - n_x^2\right) = -0.36,$ so $H$ varies as $r$ varies. A little calculation (I find a software spreadsheet is useful for this) shows that $\tan\theta \approx 1.13278$ when $r = 1$ but $\tan\theta \approx 0.71368$ when $r = 4.$ That is, the angle $\theta$ does depend on $r.$

$\endgroup$

Answer 4

In this solution, I used the parametric equation of the circle, to generate the equation of the cone whose base is the circle. Then I intersect this cone, with the projection plane $z = 1$.

The center of the circle is

$C = (C_x, C_y, C_z) $

and the unit normal vector to its plane is

$\hat{n} = (n_x, n_y, n_z) $

It is assumed that the circle does not lie in a plane not perpendicular to the projection plane, i.e. $ n_z \ne 0 $

Define the following two unit vectors

$u_1 = \dfrac{(n_z, 0, -n_x)}{ \sqrt{n_x^2 + n_z^2 } } $

$u_2 = \hat{n} \times u_1 $

Then, the circle can be written in parametric form as

$ P(t) = C + r \cos(t) u_1 + r \sin(t) u_2 $

Since the focal point is at $(0,0,0)$ then, by connecting $P(t)$ to the origin we get a cone, where each line on its lateral side is given by $s P(t) $ with $s \in [0, 1] $

Let $ \omega(t) =[ 1 , \cos(t), \sin(t) ] $, then

$P(t) = [ C , r u_1, r u_2 ] \omega(t) $

Note that $ \omega^T(t) Q_0 \omega(t) $ is identically zero, where

$Q_0 = \text{diag}( -1, 1, 1 )$

Now, define the $3 \times 3 $ matrix $A = [C, r u_1, r u_2] $, then

$\omega(t) = A^{-1} P(t)$

And since $ \omega^T(t) Q_0 \omega(t) = 0$, then the equation of the cone is

$ P^T A^{-T} Q_0 A^{-1} P = 0 $

Now, let $Q = A^{-T} Q_0 A^{-1} $, then the equation of our cone is

$ P^T Q P = 0 $

Now intersect this cone with the plane $z = 1$ by setting the $z$ coordinate of $P$ to $1$, you get the equation of a conic in $x$ and $y$, namely, if you segment $Q$ as follows

$ Q = \begin{bmatrix} Q_{22} && b \\ b^T && c \end{bmatrix} $

and define $r = [x, y]^T $, then the equation of the conic section is

$ r^T Q_{22} r + 2 r^T b + c = 0 $

And this conic equation results in the following general-form equation

$ A x^2 + B xy + C y^2 + D x + E y + F = 0 $

which can be analysed using standard methods.