In case anyone is once interested in this topic, I think I was now able to answer my question. I had especially problems with finding the right description for a "random orientation" when asking this question. The clue idea was that we can normalize, transform and turn two vectors in space, such that both vectors start from the origin and one of the two vectors is aligned with the z-axis, without affecting the probability distribution of the angle between the two vectors. The "random-orientation" of the second vector is then a point density on the unit sphere and the angle between the two vectors coincides with the polar angle of the spherical coordinate system.
If we now assume the two independent random angles $\Theta'$ and $\Phi'$, where $\Theta'$ is distributed in the yz-plane according to $f_{\Theta'}(\theta')$ and $\Phi'$ is distributed in the xy-plane according to $f_{\Phi'}(\phi')$, we find the probability of a point on the unit sphere contained in a differential area according to
$$f(\theta', \phi') sin(\theta') d\theta' d\phi' = f_{\Theta'}(\theta') f_{\Phi'}(\phi') sin(\theta') d\theta' d\phi',$$
where $sin(\theta') d\theta' d\phi'$ is the surface area element of the spherical coordinates. Therefore we can get the probability distribution of $\theta$ from
$$f_{\Theta}(\theta) d\theta = \int_{\phi} f_{\Theta'}(\theta) f_{\Phi'}(\phi') sin(\theta) d\theta' d\phi' = f_{\Theta'}(\theta) sin(\theta) d\theta$$ as $$f_{\Theta}(\theta) = f_{\Theta'}(\theta) sin(\theta).$$
In particular we can asume $f_{\Theta'}(\theta')$ to approximately be Gaussian in the interval $[0,\pi]$.
Please correct me if this is wrong..
You can’t really measure the “vertical” angle between the vectors via orthogonal projection because that distorts the angle. The underlying reason that projecting onto the $x$-$y$ plane (in the second coordinate system in the question) works for the horizontal angles but in general there’s no orthogonal projection that will get the vertical angles right is that all of the yaw rotations use the same rotation axis, but pitch rotations don’t.
You can, however, project the vectors/points onto the unit sphere, i.e., convert to spherical coordinates: $$\theta=\arccos{z\over\sqrt{x^2+y^2+z^2}} \\ \varphi = \arctan\frac yx.$$ (We don’t care about $\rho$, the distance from the origin, here.) You need to select the correct quadrant for the arc tangent, but math packages have a built-in variant of this function that takes care of this for you. Note, too, that $\theta$ increases downwards from the $z$-axis. Since there’s no camera roll, the horizontal visual angular separation is just the difference between the azimuths $\varphi$ and the vertical visual angular separation the difference in inclination $\theta$.
Using your example of $a=(4,3,1)$ and $b=(2,7,9)$, we have $\varphi_a=\arctan\frac34\approx36.870°$ and $\varphi_b=\arctan\frac72\approx74.055°$, for a difference of approximately $37.185°$, well within the field of vision. This matches the angle obtained from the projections onto the $x$-$y$ plane. It’s probably less computationally expensive to use the latter method since you can compare directly against $\cos 80°$ instead of computing an arc cosine, but this really only works because the visual field is horizontally symmetric. If you had different angles to the left and right, using the arc tangent is better since cosine doesn’t distinguish between different directions.
For the vertical angles, we get $\theta_a=\arctan 5\approx78.690°$ and $\theta_b=\arctan{\sqrt{53}\over9}\approx38.969°$ for a difference of $39.721°$, also within the field of view.
It’s also possible to compare the vertical angles by rotating $b$ so that it lies in the same vertical plane as $a$, but that seems like more work than just converting to spherical coordinates.
If you have to check a lot of points for visibility, computing all of those inverse trigonometric functions could get expensive. The field of view is bounded by four planes through the origin, so an alternative approach is to compute the equations $ax+by+cz=0$ ($\mathbf n\cdot\mathbf x=0$ in vector form) of these planes that make up the sides of this pyramid. The sign of dot product of the normal to the plane $\mathbf n$ with an arbitrary point tells you on which side of the plane the point lies: if positive, then the point lies in the direction of $\mathbf n$ from the plane; if negative, it’s on the opposite side; if zero, the point, of course, lies on the plane. You can arrange for all of the normals to be inward-pointing, so that a point is visible iff all four of these test values are non-negative.
Best Answer
One physicist to another, let's prove a generalization, $\Bbb E\cos^2\theta=\frac1n$ in $n$ dimensions, with an argument any physicist would love. Without loss of generality, one vector is $\hat{e}_1$ and the other is a unit vector, $\sum_{i=1}^nx_i\hat{e}_i$ with $\sum_{i=1}^nx_i^2=1$. Then $\cos\theta$ is the dot product $x_1$, so we want $\mu:=\Bbb Ex_1^2$. By rotational symmetries, all the $x_i^2$ have the same mean $\mu$. Since means are additive, $n\mu=\Bbb E\underbrace{\sum_ix_i^2}_{1}=1$, so $\mu=\frac1n$ as claimed.