$\newcommand{\+}{^{\dagger}}% \newcommand{\angles}[1]{\left\langle #1 \right\rangle}% \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}% \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}% \newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,}% \newcommand{\dd}{{\rm d}}% \newcommand{\ds}[1]{\displaystyle{#1}}% \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}% \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}% \newcommand{\fermi}{\,{\rm f}}% \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}% \newcommand{\half}{{1 \over 2}}% \newcommand{\ic}{{\rm i}}% \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow}% \newcommand{\isdiv}{\,\left.\right\vert\,}% \newcommand{\ket}[1]{\left\vert #1\right\rangle}% \newcommand{\ol}[1]{\overline{#1}}% \newcommand{\pars}[1]{\left( #1 \right)}% \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}}% \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}% \newcommand{\sech}{\,{\rm sech}}% \newcommand{\sgn}{\,{\rm sgn}}% \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}}% \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$ Let's assume we have a mechanism to generate random numbers in $\left[0,1\right)$ and ${\rm P}\pars{\Omega_{\vec{r}}}$ is a distribution function for random points in a sphere of radius $a > 0$. $\Omega_{\vec{r}}$ is the solid angle. In this case, ${\rm P}\pars{\Omega_{\vec{r}}}$ is, indeed, $\Omega_{\vec{r}}\,\,$-independent: $ 1 = \int{\rm P}\pars{\Omega_{\vec{r}}}\,\dd\Omega_{\vec{r}} = {\rm P}\pars{\Omega_{\vec{r}}}\int\dd\Omega_{\vec{r}} = {\rm P}\pars{\vec{r}}\pars{4\pi} \quad\imp\quad{\rm P}\pars{\vec{r}} = {1 \over 4\pi} $ Then, $ 1 = \int_{0}^{\pi}\half\,\sin\pars{\theta}\,\dd\theta\int_{0}^{2\pi} \,{\dd\phi \over 2\pi} $ We can generate random numbers $\xi_{\theta}$ and $\xi_{\phi}$ such that: $ \bracks{~\half\,\sin\pars{\theta}\,\dd\theta = \dd\xi_{\theta}\,, \quad\xi_{0} = 0 \imp \theta = 0~}\ \mbox{and}\ \bracks{~{\dd\phi \over 2\pi} = \dd\xi_{\phi}\,,\quad\xi_{0} = 0 \imp \phi = 0~} $ Those relations yield: $\ds{\half\bracks{-\cos\pars{\theta} + 1} = \xi_{\theta}}$ $\ds{\pars{~\mbox{or/and}\ \sin\pars{\theta/2} = \root{\xi_{\theta}}~}}$ and $\ds{\phi = 2\pi\,\xi_{\phi}}$: $\color{#0000ff}{\large% \theta = 2\arcsin\pars{\root{\xi_{\theta}}}\,,\qquad \phi = 2\pi\xi_{\phi}} $ As we mentioned above, $\xi_{\theta}$ and $\xi_{\phi}$ are uniformly distributed in $\left[0, 1\right)$.
For a sphere of radio $a$ the random points are given by: $ \left\lbrace% \begin{array}{rcl} \color{#0000ff}{\large x} & = & a\sin\pars{\theta}\cos\pars{\phi} = \color{#0000ff}{\large 2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}\cos\pars{2\pi\xi_{\phi}}} \\ \color{#0000ff}{\large y} & = & a\sin\pars{\theta}\sin\pars{\phi} = \color{#0000ff}{\large 2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}\sin\pars{2\pi\xi_{\phi}}} \\ \color{#0000ff}{\large z} & = & a\cos\pars{\theta} = \color{#0000ff}{\large a\pars{1 - 2\xi_{\theta}}} \end{array}\right. $