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.