Pinhole camera projection of 3D multivariate Gaussian

Consider a 3D Gaussian with $$3\times 1$$ mean $$\boldsymbol \mu$$ and $$3\times 3$$ covariance $$\boldsymbol \Sigma$$ (which is symmetric positive semidefinite):

$$p(\mathbf x) = \frac{1}{\sqrt{\operatorname{det}(2\pi\boldsymbol \Sigma)}}\exp\left(-\frac{1}{2}(\mathbf x - \boldsymbol \mu)^T \boldsymbol \Sigma^{-1} (\mathbf x - \boldsymbol \mu)\right)$$

We want to apply perspective projection to project it to 2D space with some homogeneous projection matrix $$\mathbf A$$ which is $$3 \times 4$$ which projects the 3D vector $$\begin{bmatrix}x & y & z\end{bmatrix}^T$$ to 2D vector $$\begin{bmatrix}u & v\end{bmatrix}^T$$. In homogeneous coordinates,

$$\begin{bmatrix}u\\ v\\ 1\end{bmatrix} \equiv \mathbf A \begin{bmatrix}x \\ y\\ z\\ 1\end{bmatrix}$$

where $$\equiv$$ means "equal up to multiplication by some scalar".

The question is, of course, assuming some 3D random variable $$\mathbf x$$ is distributed with some 3D Gaussian, what is the distribution of its projection via $$\mathbf A$$?

If it makes it easier, we may simplify the problem by looking at the special case of the "pinhole camera model" --- that is, $$\mathbf A = \mathbf K \begin{bmatrix}\mathbf R & \mathbf t\end{bmatrix}$$ where

• $$\mathbf K$$ is a $$3\times 3$$ matrix of the form $$\mathbf K = \begin{bmatrix}f_x & 0 & m_x \\ 0 & f_y & m_y \\ 0 & 0 & 1\end{bmatrix}$$ where $$f_x$$, $$f_y$$, $$m_x$$, $$m_y$$ are arbitrary scalars,
• $$\mathbf R$$ is a $$3\times 3$$ matrix in the special orthogonal group $$\operatorname{SO}(3)$$ (i.e. it is orthonormal and its determinant is $$+1$$), and
• $$\mathbf t$$ is some arbitrary $$3\times 1$$ vector.

For this special case, here is what I have so far:

The mean is just typical projection.

$$\boldsymbol \mu_{2D} = \begin{bmatrix}a/c \\ b/c\end{bmatrix} \text{ where } \begin{bmatrix}a \\ b \\ c \end{bmatrix} = \mathbf A \begin{bmatrix}\boldsymbol \mu \\ 1\end{bmatrix}$$

The covariance would seem to be

$$\boldsymbol \Sigma_{2D} = \frac{1}{c^2} \mathbf K_{2\times 2} \boldsymbol \Sigma_{\text{rotated }2\times 2} \mathbf K_{2\times 2}^T$$

where $$c$$ is as above and $$\boldsymbol \Sigma_{\text{rotated }2\times 2}$$ is the top-left $$2\times 2$$ submatrix of:

$$\boldsymbol \Sigma_{\text{rotated}} = \mathbf R \boldsymbol \Sigma \mathbf R^T$$ and $$\mathbf K_{2\times 2}$$ is the top-left submatrix of $$\mathbf K$$: $$\mathbf K_{2\times 2} = \begin{bmatrix}f_x & 0\\0 & f_y\end{bmatrix}$$

However, I'm not sure if my calculation for the covariance is correct. In particular, intuitively, variance in the "depth" direction should still exist in the projected "image-space" but this is not captured in my calculation.

In general, the distribution after perspective projection will no longer be a Gaussian.

However, if you want the output to be a Gaussian, we can approximate it by linearization.

Let the Jacobian of $$\mathbf F(\mathbf x) = \begin{bmatrix}a / c \\ b / c\end{bmatrix}$$ where $$\begin{bmatrix}a & b & c\end{bmatrix}^T = \mathbf A \mathbf x$$ be

$$\mathbf J = \begin{bmatrix} \frac{\partial \mathbf F(\mathbf x)}{\partial x_1} & \frac{\partial \mathbf F(\mathbf x)}{\partial x_2} & \frac{\partial \mathbf F(\mathbf x)}{\partial x_3} \end{bmatrix}$$

Then the covariance would simply be $$\mathbf J \boldsymbol \Sigma \mathbf J^T$$.

In the case where $$\mathbf A$$ is the pinhole camera model as above, let us compute the Jacobian of projection by the intrinsic matrix, i.e. of the function $$\mathbf F(\mathbf x) = \begin{bmatrix}a / c \\ b / c\end{bmatrix}$$ where $$\begin{bmatrix}a & b & c\end{bmatrix}^T = \mathbf K \mathbf x$$ and $$\mathbf x = \begin{bmatrix}x & y & z\end{bmatrix}^T$$. Recall from the original question that

$$\mathbf K = \begin{bmatrix}f_x & 0 & m_x\\0 & f_y & m_y\\0 & 0 & 1\end{bmatrix}$$

such that

$$\mathbf F\left(\begin{bmatrix}x\\y\\z\end{bmatrix}\right) = \begin{bmatrix}\frac{f_x x}{z} + m_x \\\frac{f_y y}{z} + m_y\end{bmatrix}$$

In this case, the Jacobian is straightforwardly computed:

$$\mathbf J_{\mathbf K} = \begin{bmatrix} \frac{f_x}{z} & 0 & -\frac{f_x x}{z^2}\\ 0 & \frac{f_y}{z} & -\frac{f_y y}{z^2} \end{bmatrix}$$

If you want to take into account the extrinsic matrix $$\begin{bmatrix}\mathbf R & \mathbf t\end{bmatrix}$$, you can simply use $$\boldsymbol \Sigma_{\text{rotated}}$$ as above.