Question
I have an ellipsoid described by a covariance matrix $\Sigma \in \mathbb{R}^{3\times3}$ and its centroid location $\mu \in \mathbb{R}^{3}$ (yes; notation from statistics as I want to interpret it as a Gaussian later). I want to obtain the projection of it in a geometry given by the projective mapping $P \in \mathbb{R}^{3\times4}$. Within the constraint geometrical setup detailed later, i believe the conic projection should be an ellipse again and that should be expressable with a 2D ($\Sigma, \mu$).
Available as Notebook.
Approach so far
I followed this answer but had no luck transferring the approach to a non-parrallel projection using a quadric where the translation is non-zero.
homogenous represention of quadric
I start by expressing the ellipsoid as a quadric following the answer in [4]. A quadric is described by a symmetrical $4\times4$ matrix $Q$:
$$\hat{\mathbf x}^T\mathtt Q\hat{\mathbf x} = \begin{bmatrix}\mathbf x&1\end{bmatrix} \begin{bmatrix}\mathtt A & \mathbf b \\ \mathbf b^T & c\end{bmatrix} \begin{bmatrix}\mathbf x\\1\end{bmatrix} = 0.$$
I plug in the covariance and mean vector for the respective variables in $Q$ to yield:
$$\mathtt Q = \begin{bmatrix}\mathtt \Sigma & \mathbf \mu \\ \mathbf \mu^T & -1\end{bmatrix}.$$
construction of projection matrix
I want to create projection matrices that describe a projective geometry given by the focal length $f$, the pierce point $(p_x, p_y)$ the source position $C$ and a camera rotation $R_c$. I obtain the projection matrix $P$ following [1]:
$$
P=
\begin{bmatrix}
f & 0 & p_{x} \\
0 & f & p_{y} \\
0 & 0 & 1
\end{bmatrix}
\times
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0\\
0 & 0 & 1 & 0
\end{bmatrix}
\times
\begin{bmatrix}
\mathtt R_{c}^T & -R_{c}^TC \\
\mathtt 0 & 1
\end{bmatrix}
$$
projection of quadric
I followed the answer in [4] which quotes the book [2] Hartley Zisserman chapter 8.3 to get this simple formula for the projection of the quadric $Q$ onto the conic $C$:
$$\mathtt C^*=\mathtt P\mathtt Q^*\mathtt P^T$$
The dual-forms (asterisk) are obtained by the use of the inverse ($C^{*}=C^{-1}$):
$$C = (PQ^{-1}P^T)^{-1}$$
Example
The projective geometry is given by $f=3718$, $C=[0, 0, -1987]^T$, $R=I$. This describes a camera looking into positive z-direction sitting at some point on the z-axis.
I define an ellipsoid aligned with the coordinate axis with dimensions $(300, 100, 100)$ and position $(255.5, 255.5, 255.5)$:
$$\mathtt Q =
\begin{bmatrix}
300 & 0 & 0 & 255.5 \\
0 & 100 & 0 & 255.5 \\
0 & 0 & 100 & 255.5 \\
255.5 & 255.5 & 255.5 & -1
\end{bmatrix},$$
To verify the formulation of the quadric, I generate a set of points on the unit-sphere and multiply them with the quadric $\hat{x}'=Q\hat{x}$. The set of points traces out the desired ellipsoid nicely.
To verify the projection matrices, I project these points onto the detector. The resulting scatter plot looks convincing.
To test if the computed Conic is valid, I create points on the unit-circle and multiply them with the conic $\hat{u}'=C\hat{u}$. This however looks false. I observe, that the Covariance of the projected conic is almost zer0.
Does someone have a clue what could be going wrong? Is the theory consistent, or am I missing something important?
Sources
[1]: Construct Projection Matrix, https://ksimek.github.io/2012/08/22/extrinsic/
[2]: Hartley & Zisserman, Multiple View Geometry In Computer Vision, 2004
[3]: Project Quadric onto Conic, https://math.stackexchange.com/a/2757064
[4]: Construct Quadric from Cov and mean, https://math.stackexchange.com/a/3201266/1018514
Best Answer
In general the approach $C = (P^{\top}Q^{-1}P)^{-1}$ is valid, however there are some intricacies that need to be dealt with carefully:
1. Quadric from Gaussian
I cite my other post that deals with this. To describe the surface-ellipsoid of a 3D Gaussian thresholded at $\sigma$, use this parametrization of $Q$:
$$(x-\mu)^\top\Sigma^{-1}(x-\mu)-\sigma=0$$
$$Q = \begin{bmatrix}A & b \\ b^⊤ & c\end{bmatrix}, \quad \hat{x}^{\top}Q\hat{x}=0, \quad \hat{x}=[x^{\top},1]^{\top}$$
$$A=\Sigma^{-1}, \quad b=-A\mu , \quad c=\mu^⊤\Sigma^{-1}\mu-\sigma$$
2. Gaussian from Conic
The conic $C$ is derived from $Q$ by $C = (P^{\top}Q^{-1}P)^{-1}$. To represent is as a 2D Gaussian distribution $(\Sigma', \mu')$ thresholded at $\sigma$, we derive the parameters analogously to above.
$$C = \begin{bmatrix}a&b/2&d/2\\b/2&c&e/2\\d/2&e/2&f\end{bmatrix} = \begin{bmatrix}A'&b'\\b'^{\top} & f\end{bmatrix}$$
$$\Sigma=A^{-1}, \quad \mu =-\Sigma b', \quad \sigma = \mu^{\top} \Sigma^{-1} \mu - f$$
What I have overlooked so far is to scale the Covariance by $\sigma$ (at least thats my interpretation of it, but please correct me if I am wrong):
$$\Sigma'=\Sigma\sigma$$
3. Sampling the $\sigma$-surface/line of a 3D/2D Gaussian
Finding points on the surface of a 3D Gaussian or analogously the circumference of a 2D Gaussian can be done by using the Cholesky Factor of the Covariance. I like to think of it as the square root for matrices:
$$\Sigma=LL^{\top}$$
Points on the Ellipsoid/Ellipse at $\sigma=1$ can now be obtained by transforming points on the unit-sphere/circle $\hat{x}$ like so:
$$x=L\hat{x}+\mu'$$
Sampling directly using $\Sigma$ will also produce an ellipsoid/ellipse, but the radii of the resulting object will be squared and thus not correspond with the formulation of the quadric/conic $Q$/$C$.