Quick contextualisation, I am trying to follow the orientation of a molecule (represented as a rigid object) in respect to the z axis. The molecule is basically shaped as an hexagon, so what i thought was to define 2 vectors: $\vec u$ running from one vertex to the opposite one of the hexagon and a $\vec v $, orthogonal to $\vec u$ . What i do then is simply measuring their angle with the z axis, which i call respectively $\theta$ and $\psi$ and counting the number of occurence of every couple ( $\theta$ , $\psi$ ). With that i fill a 180x180 matrix with the number of occurence of every couple, then divide by the total number of observations. When i plot this matrix i get a typically shaped diamond plot with the probability on the edges showing higher values due to the correlation of the two vectors and the geometrical density of the states (and with this i mean the that the angle of a vector with a fixed axes has a non uniform density distribution).
The specific thing that i'm looking for is to correct for this "geometrical" effect so that when the molecule is freely rotating in respect to the z axis then the 3D plot of the angle probability looks "flat" inside the diamond, or if you want "even".
To better explain, if i was to measure just the probability of occurence of one angle, say $\theta$, i would divide every probability value of $\theta$ by the $sin(\theta)$ to "correct" for the geometrical effect, and than the plot would be evened out (although ofc at 0 you have sin(0)=0 which creates a problem which is still solvable by binning). But if i have two angles i guess a conditional probability is necessary to normalize my distribution since $\theta$ and $\psi$ are not independent, and i can't get my head around what it should be.
Can anyone help by any chance? I attach an image of the angle map that i get and that i would want to flat out:
Edit: After correction i get the map flattened as wanted. The peaks i think are due to the machine calculation of the radicand which is not perfect.


Here's a new answer, based on the clarifications in the comments below.
You clarified that what you're interested in isn't really normalization, but rather a reference density that would be expected for a free molecule so that you can divide this out in order to visualize effects that cause deviations from this free-molecule density.
Let's regard the molecule as fixed and the reference axis for $\theta$ and $\psi$ as freely rotating; that is, let's view the angles $\theta$ and $\psi$ as the angles that two fixed orthogonal vectors form with a uniformly random direction. Let the fixed $x$ and $y$ axes point along the two orthogonal vectors. A uniformly random unit direction vector has independent uniform distributions of $z$ over $[-1,1]$ and of $\phi$ (the angle about the $z$ axis) over $[0,2\pi]$. (Here the $z$ axis is not the random reference axis, but the fixed coordinate axis orthogonal to the $x$ and $y$ axes.) The transformation from $z,\phi$ to $x,y$ is
$$ x=\sqrt{1-z^2}\cos\phi\;,\\ y=\sqrt{1-z^2}\sin\phi\;. $$
Thus
\begin{eqnarray*} \frac{\partial(x,y)}{\partial(z,\phi)} &=&\left|\matrix{\frac{\partial x}{\partial z}&\frac{\partial x}{\partial\phi}\\\frac{\partial y}{\partial z}&\frac{\partial y}{\partial\phi}}\right| \\ &=&\left|\matrix{\frac z{\sqrt{1-z^2}}\cos\phi&-\sqrt{1-z^2}\sin\phi\\\frac z{\sqrt{1-z^2}}\sin\phi&\sqrt{1-z^2}\cos\phi}\right| \\ &=& z \\ &=& \sqrt{1-x^2-y^2}\;. \end{eqnarray*}
The further transform according to $x=\cos\theta$ and $y=\cos\psi$ has Jacobian
$$ \frac{\partial(x,y)}{\partial(\theta,\psi)}=\sin\theta\sin\psi\;, $$
so in all we have
$$ \frac{\partial(z,\phi)}{\partial(\theta,\psi)}=\frac{\sin\theta\sin\psi}{\sqrt{1-\cos^2\theta-\cos^2\psi}}\;. $$
This is the factor you need to divide by, since the density for $z$ and $\phi$ is flat (as these variables are independently uniformly distributed), and the densities are releated by
$$ f_{z,\phi}(z,\phi)\mathrm dz\mathrm d\phi=f_{\theta,\psi}(\theta,\psi)\mathrm d\theta\mathrm d\psi $$
and thus
$$ f_{\theta,\psi}(\theta,\psi)=\frac{\partial(z,\phi)}{\partial(\theta,\psi)}f_{z,\phi}(z,\phi)\propto\frac{\partial(z,\phi)}{\partial(\theta,\psi)}\;. $$
Here $\theta$ and $\psi$ are limited to values for which the radicand in the denominator is non-negative, that is, to the diamond
$$ \left|\,\theta-\frac\pi2\right|+\left|\,\psi-\frac\pi2\right|\le\frac\pi2\;. $$
This result fits with your image. It's restricted to the diamond; it shows roughly a $\sin\theta\sin\psi$ behaviour along the boundary; and the divergence at the boundary (mitigated by your binning procedure) is recognizable.