3 random numbers to describe point on a sphere [duplicate]

Question

I'm currently working on a problem involving computer graphics and got into a discussion about whatever or not constructing a 3d vector out of 3 random points uniformly distributed points between -1 and 1 (and then normalize the vector) to get points uniformly distributed over the surface of the sphere?

Note that I'm a computer science student not a math student and as such might not be able to follow the really complex stuff.

The other student said it would be better to instead pick 2 random points on a 2d plane and then warp those onto a sphere, however this to me seemed needlessly complex which one of us is right?

Answer 1

Methods to do this are discussed here: Wolfram MathWorld: Sphere Point Picking

In particular a programatically elegant approach is suggested by Marsaglia, stating that with $x_1,x_2\in(-1,1)$ chosen uniformly on this interval and throwing away cases where $x_1^2+x_2^2\geq 1$ we can define $$ \begin{align} x&=2x_1\sqrt{1-(x_1^2+x_2^2)}\\ y&=2x_2\sqrt{1-(x_1^2+x_2^2)}\\ z&=1-2(x_1^2+x_2^2) \end{align} $$ to obtain points $(x,y,z)$ uniformly distributed on the unit sphere. Quite clever!

The algorithm takes a point $(x_1,x_2)$ on the unit disc and transforms it to a point $(x,y,z)$ on the unit sphere. Here is a Dynamic Figure Illustrating this Transformation:

The purple circle in my diagram has radius $\sqrt{0.5}$ thus dividing the area of the unit disc in two halves - just like the tranformation of it divides the area of the unit sphere in two halves. This illustrates the main principle well, namely:

A uniform distribution on the unit disc is transformed into a uniform distribution on the unit sphere by Marsaglia's formulas.

---

Update

I just thought of this simple method which needs no advanced proof:

Pick $x,y,z\in[-1,1]$ which is what you already suggested. Now discard them if $x^2+y^2+z^2>1$ meaning that $(x,y,z)$ lie in the part of the unit cube that is outside of the unit sphere so in one of the corners. If not, normalize the vector $(x,y,z)$.

To avoid division by zero and precision errors you could also discard points with $x^2+y^2+z^2<\delta$ for some appropriate $0<\delta<1$ (as others have suggested in various comments).

Answer 2

Choosing the points uniformly from $[-1,1]^3$ and then normalizing does not give the uniform distribution on the sphere. The poles will get less mass than corners.

One way to get a uniform distribution would be to choose the $3$ coordinates with Gaussian distribution in $\mathbb{R}$ and then normalize.

Answer 3

Your vector will be uniformly distributed over a cube. The cube's vertex, however, is farhter from the inscribed sphere than the cube side's center, so a small region near the vertex will project to a smaller part on the sphere than the same volume region near the side's center. So the spherical distribution of the resulting points over the sphere will not be that uniform.

Anyway it may be close enough to your task, depending on what you need...

Or you may use any area-preserving 2-dimensional projection known in cartography - see https://en.wikipedia.org/wiki/Map_projection#Equal-area

Answer 4

If you're stuck to uniformly distributed PRNG, another option would be to use spherical coordinates, but transform uniformly distributed pairs of numbers to compensate for coordinate distortion.

Because of coordinate distortion, if you take just uniformly distributed pairs $(\vartheta,\phi)$ of numbers as spherical coordinates, you'll get too many points near poles. Recall that when integrating over a sphere, we use a Jacobian determinant of $\sin\theta$. So, we should have $P(\theta)\sim\sin\theta$. To get it, we have to take its inverse cumulative distribution function of our $\vartheta$.

Since $$\int_0^\vartheta \sin\vartheta~\text{d}\vartheta\sim\sin^2\frac{\vartheta}2,$$

and its inverse has the form of

$$\text{CDF}^{-1}(P)=\arcsin\sqrt P,$$

we can take for $\vartheta\in[0,\pi]$:

$$\theta=\arcsin\sqrt\frac\vartheta\pi.$$

Since Jacobian determinant doesn't depend on $\phi$, we don't have to transform this variable.

Now the pair $(\theta,\phi)$ is a pair of uniformly distributed on sphere spherical coordinates.

Answer 5

Since Marsaglia's formula for generating random points on a sphere is given in String's answer, I thought it would be useful to justify the formula.

In this answer, it is shown that the surface area of an annulus $a\le z\le b$ on the sphere is the same as the area of the strip $a\le z\le b$ on the cylinder tangent to the $x$-$y$ equator of the sphere. This means that the cylindrical projection of points uniformly distributed on the sphere, must also be uniformly distributed on the cylinder. That is,

The $z$-coordinate of the random points should be uniformly distributed.

Let $u_1$ and $u_2$ be uniformly distributed in $[-1,1]$. Consider $\left.u_1^2+u_2^2\,\,\middle|\,\,u_1^2+u_2^2\le1\right.$. Given that $u_1^2+u_2^2\le1$, $(u_1,u_2)$ is uniformly distributed in the unit disk, and therefore, $$ P[u_1^2+u_2^2\le r^2]=r^2\tag{1} $$ Equation $(1)$ shows that $u_1^2+u_2^2$ is uniformly distributed in $[0,1]$, and therefore,

$1-2(u_1^2+u_2^2)$ is uniformly distributed in $[-1,1]$

Therefore, we can let $z=1-2(u_1^2+u_2^2)$. Then, $$ \sqrt{1-z^2}=2\sqrt{u_1^2+u_2^2}\sqrt{1-(u_1^2+u_2^2)}\tag{2} $$

Since $(u_1,u_2)$ is uniformly distributed in the unit disk,

$\dfrac{(u_1,u_2)}{\sqrt{u_1^2+u_2^2}}$ is uniformly distributed on the unit circle.

Thus, $(2)$ and the highlighted statements shows that given $u_1^2+u_2^2\le1$, $$ \begin{align} x&=2u_1\sqrt{1-(u_1^2+u_2^2)}\\ y&=2u_2\sqrt{1-(u_1^2+u_2^2)}\\[4pt] z&=1-2(u_1^2+u_2^2) \end{align}\tag{3} $$ is uniformly distributed on the unit sphere.

---

George Marsaglia's formula $(3)$ requires an average of $8/\pi\doteq2.546479$ random numbers per point. If it is faster to evaluate a cosine and a square root than to generate $0.546479$ random numbers, then generate $u_1,u_2$ uniformly in $[-1,1]$.

Compute $\cos(\pi u_2)$, then $\sin(\pi u_2)=\mathrm{sgn}(u_2)\sqrt{1-\cos^2(\pi u_2)}$. $$ \begin{align} x&=\sqrt{1-u_1^2}\cos(\pi u_2)\\ y&=\sqrt{1-u_1^2}\sin(\pi u_2)\\[4pt] z&=u_1 \end{align}\tag{4} $$ is uniformly distributed on the unit sphere, without discarding any random numbers.

Answer 6

Sounds like you want a point on the surface of the sphere, not in its volume. As already pointed out, choosing three numbers in $[-1, 1]$ won't give you an uniform distribution.

A surface has only two degrees of freedom, so you should need only two random numbers. This solution uses just two uniformly distributed numbers and no rejection.

Let's assume you want the sphere of radius $R$. First, pick a number $z$ between $-R$ and $R$ inclusive, and an angle $\theta$ so that $0 \le \theta < 2\pi$. Now calculate $x=\sqrt{R^2-z^2} \cos \theta$ and $y=\sqrt{R^2-z^2} \sin \theta$. The point $(x, y, z)$ is uniformly distributed over the surface of the sphere.

But if you want it in its volume, you can pick a uniform random radius ratio $r$ so that $0 < r \le 1$ and divide the coordinates of the point on the surface already obtained, by the cube root of $r$. That will give you a point uniformly distributed in the volume of the sphere.

Edit: To be clear, that last step is: $q=r^{1/3}, P=(x/q,y/q,z/q)$

Answer 7

Normalizing the vectors amounts to projecting the random 3D points onto the surface.

The probability to reach a given disk on the sphere is proportional to the volume of the portion of the cone from the center through the disk, that is contained in the unit cube. This volume increases as the third power of the length of the ray from the center of the sphere, through the center of the disk, into the first face.

For a disk facing a face, the probability is proportional to $1$; for a disk facing an edge, $\sqrt8\approx2.8$; and for a disk facing a corner, $\sqrt{27}\approx5.2$. [More generally, for a disk centered on $(x, y,z)$ on the sphere, the probability is proportional to $\frac{1}{(\max(|x|,|y|,|z|))^3}$.]

This probability is fairly non-uniform.

Answer 8

$\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{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \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]{\vphantom{\large A}\,#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} \newcommand{\wt}[1]{\widetilde{#1}}$ Lets ${\rm P}\pars{\vec{r}}$ the probability distribution of random points uniformly distributed in the surface of a sphere. $\ds{a}$ is the sphere radius. It satisfies ( note that $\ds{{\rm P}\pars{\vec{r}}}$ is, indeed, $\ds{\vec{r}}$-independent ): $$ 1 = \int_{S}{\rm P}\pars{\vec{r}}\,\dd S = {\rm P}\pars{\vec{r}}\int_{S}\dd S ={\rm P}\pars{\vec{r}}\pars{4\pi a^{2}}\quad\imp\quad{\rm P}\pars{\vec{r}} ={1 \over 4\pi a^{2}} $$

Then, $$ 1 = \int_{S}{\dd S \over 4\pi a^{2}} =\int_{0}^{\pi}\half\sin\pars{\theta}\,\dd\theta\int_{0}^{2\pi}{\dd\phi \over 2\pi} $$ which we can visualize as two independent probability distributions $\ds{\half\,\sin\pars{\theta}}$ and $\ds{1 \over 2\pi}$.

Now, we generate random numbers $\ds{\xi_{\theta}}$ and $\ds{\xi_{\phi}}$ uniformly distributed in $\ds{\left[0,1\vphantom{\large A}\right)}$ such that

\begin{align} \int_{0}^{\theta}\half\,\sin\pars{t}\,\dd t=\int_{0}^{\xi_{\theta}}\,\dd t\quad &\imp\quad\sin\pars{\theta \over 2} = \root{\xi_{\theta}}\quad\imp\quad\theta = 2\arcsin\pars{\root{\xi_{\theta}}} \\[3mm] \int_{0}^{\phi}{1 \over 2\pi}\,\dd t=\int_{0}^{\xi_{\phi}}\,\dd t\quad &\imp\quad \phi = 2\pi\xi_{\phi} \end{align} Note that $$ \sin\pars{\theta} = 2\sqrt{\xi_{\theta}\pars{1 - \xi_{\theta}}} \qquad\mbox{and}\qquad\cos\pars{\theta} = 1 - 2\xi_{\theta} $$

Point $\ds{\pars{x,y,z}}$ is given by $$ \begin{array}{rcrcl} x &= &a\sin\pars{\theta}\cos\pars{\phi}&=&2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}} \cos\pars{2\pi\xi_{\phi}} \\ y &=& a\sin\pars{\theta}\sin\pars{\phi}&=&2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}} \sin\pars{2\pi\xi_{\phi}} \\ z &=& a\cos\pars{\theta}&=&a\pars{1 - 2\xi_{\theta}} \end{array} $$ Remember that $\ds{\xi_{\theta}\ \mbox{and}\ \xi_{\phi}}$ are random numbers uniformly distributed in $\ds{\left[0,1\right)}$.

A typical 'run' of the $\tt C++$ script at the bottom yields:

( -0.08187 , -1.981 , -0.2642 )
( 1.562 , 0.7686 , -0.9843 )
( -0.2041 , 0.942 , 1.752 )
( -0.4777 , 1.891 , 0.4411 )
( 1.891 , 0.6153 , 0.2145 )

/* surfaceSphere_0.cc  18-jun-2014 Felix Marin
   http://math.stackexchange.com/a/839189/85343
*/
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;
const double RANDMAX1=double(RAND_MAX) + 1.0,TWOPI=2.0*M_PI;
void surfacePoint(double radius,vector<double> &p);

inline double rand01() // Return uniform values in [0,1)
{
 return rand()/RANDMAX1;
}

int main()
{
 vector<double> point;
 point.resize(3U);
 cout<<setprecision(4);
 srand(size_t(time((time_t *)0)));
 for ( unsigned char n=0 ; n<5 ; ++n ) {
     surfacePoint(2.0,point);
     cout<<"( "<<point[0]<<" , "<<point[1];
     cout<<" , "<<point[2]<<" )\n";
 }
 return 0;
}

void surfacePoint(double radius,vector<double> &p)
{
 static double temp,xiTheta,twoPiXiPhi;

 xiTheta=rand01();
 temp=2.0*radius*sqrt(xiTheta*(1.0 - xiTheta));
 twoPiXiPhi=TWOPI*rand01();
 p[0]=temp*cos(twoPiXiPhi);
 p[1]=temp*sin(twoPiXiPhi);
 p[2]=radius*(1.0 - 2.0*xiTheta);
}