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

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?

-

## marked as duplicate by leonbloy, vonbrand, Hakim, user88595, user61527 Jun 19 '14 at 22:55

I remember being bombarded with this question at a job interview years ago...they were interested in my on-the-spot analytical skills. I decided they were a bunch of train-spotters. (They were actually just a 2D graphics company...which was even more ironic). I started with x,y,z and then moved to polar coords. They weren't happy with these answers. I just wasn't happy and left ;) – cmroanirgo Jun 18 '14 at 11:44
If you're a computer science student who wants to do graphics, you better start thinking about becoming a math student or you're gonna have a bad time. – corsiKa Jun 18 '14 at 15:36
The link here: math.stackexchange.com/questions/44689/… , may be informative for your application. – jrand Jun 18 '14 at 17:00
I once answered a related question on Mathematica: mathematica.stackexchange.com/questions/13038/… – David Z Jun 18 '14 at 22:33
The problem is you wanting to use three uniformly distributed random numbers when you, in fact, only need two. By a quirk, the $z$-coordinate (or any fixed one) IS uniformly distributed. So think in terms of: A) selecting $z$, B) selecting a longitude. Both uniformly and, Voila! – Jyrki Lahtonen Jun 19 '14 at 11:29

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.

-
Why does it work to choose three normally-distributed points? – MJD Jun 18 '14 at 13:02
@MJD It works because $x^2+y^2+...+\xi^2=r^2$ is an equation of hypersphere with hyperradius $r$, and the LHS of this equation is an argument of the multivariate normal distribution. – Ruslan Jun 18 '14 at 13:29
@Ruslan - it took me a couple of minutes to grasp this, but I once I did it seems easy to understand - a simplification of the problem would be to choose a random point on a circle by picking (x,y) pairs. If x and y are normally distributed, then their root-sum-of-squares (r) is a normal distribution in R2. Thus, the probability of picking a point on any direction is equal! – ysap Jun 18 '14 at 13:44
In the old days, we would pick a point in the cube, throw away the ones outside the sphere (and the point (0,0,0) if we were unlucky enough to get it) and normalize it. This is Knuth's method (Art of Computer Programming), mostly because it is so very much cheaper in computer cycles to pick a uniform variate than a Gaussian one. (Not that this matters any more.) – Eric Towers Jun 19 '14 at 4:17
One might mention forbidding points very close to the origin (which actually have the highest probability density in the Gaussian distribution) since they may produce artifacts (or plain division by 0) during normalisation, which may perturb the uniform distribution on the sphere. – Marc van Leeuwen Jun 19 '14 at 9:03

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).

-
@String: I have provided a justification of Marsaglia's formula in my answer. I think it is a useful method since it requires an average of $8/\pi\doteq2.546479$ random numbers per point to be generated whereas the latter method requires an average of $18/\pi\doteq5.729578$ random numbers per point to be generated. – robjohn Jun 18 '14 at 20:40
@robjohn: Thank you for the insight! Your figures regarding the efficiency of each method are much more relevant than my percentages, indeed. In particular it shows that Marsaglia's method is even more efficient than the Gaussian approach which requires 3 random numbers per point. BTW the last method was added for its simplicity, not efficiency ;) – String Jun 18 '14 at 21:04
@robjohn: Why do you call this efficient? The method in David Z's and Jim Belk's answers (see comments under the OP) require exactly two uniformly distributed random numbers per point. And the resulting distribution is uniform on the sphere. – Jyrki Lahtonen Jun 19 '14 at 11:35
@JyrkiLahtonen: first, I never called this efficient; I was comparing the method in the Update with Marsaglia's method. However, transcendental functions (which can be sped up in hardware, but were cycle hogs when the formula was written) and random number generation (which are still usually library based) are the bottlenecks. David Z's answer (which I had not seen yet as it is linked from a comment to the question) calls arccos, then cos twice and sin thrice, though this can be optimized (see the end of my answer). – robjohn Jun 19 '14 at 14:36
Thanks, @robjohn. I did not consider HW efficiency, so that is fair point. Apologies. Back in the day when I had a need for this I was stuck with an integer valued PRNG, so I got away with building look-up-tables for the necessary transcendental functions. – Jyrki Lahtonen Jun 19 '14 at 15:23

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

-

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.

-

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.

-
The basic, remarkable, fact used here, is that rectangular maps (of shape $1:\pi$) of the sphere (Earth's surface) produced by projecting, from the axis joining the poles, perpendicularly outward to a cylinder wrapped around the equator is area preserving. Which is a reason to sometimes use this projection, or a stretched out version called Gall-Peters projection in cartography (even though it does deform the poles a lot). – Marc van Leeuwen Jun 19 '14 at 9:13
math.stackexchange.com/a/87243/312 – leonbloy Jun 19 '14 at 17:00
@MarcvanLeeuwen and others: I hope you will take the time to see the dynamic diagram I just added in my answer - showing the tranformation from the unit disc to the unit sphere given by Marsaglia's method. – String Jun 20 '14 at 9:40

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)$

-

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.

-

$\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;

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;
}