Random Points on a Sphere (Generating Random Sequences III)

Last week, we discussed how we could generate uniform random points in a triangle using a (tiny) bit of linear algebra, mostly the parallelogram rule, and a random variable uniform on [0,1]. It required a tiny bit of math and was computationally very simple.

This time, let us see how we can generate uniformly distributed random points on a sphere.

At first, you may think choosing independently latitude and longitude (in polar coordinates) would do the trick, but you quickly realize, after plotting it the points, that they are inordinately numerous near the poles and comparatively sparse at the equator:

So what went wrong? Choosing both coordinates independently fails to take into account that not all latitudes have the same circumference and, accordingly, their differential area varies. Simply put: there are more points to be put on a latitude near the equator than near a pole.

So a method of mapping randomly, but uniformly spaced, points on a sphere must take this basic observation into account. Let us simplify the problem by considering we are working on a unit sphere, that is, a sphere with a radius of 1. Let us suppose that the South-North axis is the Y axis. Then the sphere goes from -1 (the South Pole) to +1 (the North Pole).

To distribute points equidistantly on the sphere, we have to know the differential surface length between two given latitude, and to do so, we will need to compute a surface of revolution (that’s calculus 201 or 301).

At height y the radius of the disk perpendicular to y that reaches the surface of the sphere is given by x=\sqrt{1-y^2}, which is just Pythagoras’ formula with for a radius of 1. The differential radius is:

\displaystyle\frac{\partial\:x}{\partial\:y}=\frac{-y}{\sqrt{1-y^2}}

and

\displaystyle\left(\frac{\partial\:x}{\partial\:y}\right)^2=\frac{y^2}{1-y^2}

The integral of the surface of revolution is between latitudes a
and b:

\displaystyle{}S(a,b)=2\pi\int_a^b x\sqrt{1+\left(\frac{\partial\:x}{\partial\:y}\right)^2}\partial\:y

\displaystyle=2\pi\int_a^b \sqrt{1-y^2}\sqrt{1+\frac{y^2}{1-y^2}}\partial\:y

\displaystyle=2\pi\int_a^b\frac{1}{\sqrt{1-y^2}}\partial\:y

\displaystyle=\sin^{-1}(y)\Big|_a^b

The indefinite integral is therefore only the familiar \sin^{-1} (or “arcsin”). \sin^{-1}(y) maps -1 to 1 onto [-\pi/2,\pi/2] with equal density, which is exactly what we need (and just computed). Finally:

\theta=\mathcal{U}(-1,1)

\phi=\mathcal{U}(0,2\pi)

where, here also, \mathcal{U}(a,b) is a uniform random number on the interval [a,b]. We only have to compute the polar to Cartesian transformation (if we need to) for display:

x=\cos\theta\cos\phi

y=\sin\phi

z=\cos\theta\sin\phi

And indeed, the points are uniformly distributed:

That’s a lot of work, but it was completely worth the effort.

*
* *

Calculus, and math in general, is quite a powerful ally of the good programmer—even more so of the computer scientist—but it seems that many C.S. programs do not take math seriously enough. When doing my B.Sc. (maybe I wasn’t the best student ever, but that’s beside the point here) we only had 5 courses on a total of 30-something. We had a calculus course which had more or less the same contents as the two cégep courses and as such wasn’t that useful, a course on linear algebra that was also pretty dismal (I have quite a few stories about the lecturer), a course on probability and statistics, a course on statistics, and finally a course in discrete mathematics (the one I enjoyed the most, by far). Overall, it all seems a bit thin for people who will manipulate information for the rest of their lives, even if only to generate random numbers.

Let me end this digression by enjoining my readers, all of you, to question the quality and depth of your mathematical education and, if you’re still studying, to ask yourself what better math could do for you, and what means you should take to better your skills. I do that too, all the time, so it’s not like I’m saying I’m much better than you.

17 Responses to Random Points on a Sphere (Generating Random Sequences III)

  1. Ofek Shilon says:

    An alternative approach is –
    (1) Independently generate 3 coordinates from a *Normal* distribution,
    (2) The resulting vector distribution is dependent only on the radius (immediate check), so it suffices to normalize it.
    This approach has the advantage that it is easily generalizable to any dimension.

    • Steven Pigeon says:

      And it’s also quite easier, normalisation is very simple. Using a normal (0,1) distribution (and normalization) I get:

      From a Normal Distribution

      But if we use uniform, it is not too bad until we reach a certain density where artefacts (uneveness of distribution) become obvious:

      From Uniform Random

      Well, we have to rotate/scale to see the clumps:

      From Uniform Random

      (use right-click “view image” to see full resolution)

      • AM says:

        The “clumps” are the result of projecting a cubic volume onto a sphere. If you reject the points outside the sphere (mag>1) as you generate them, you get an even distribution.

  2. Ofek Shilon says:

    Here is a mention of another nice idea, by Marsaglia: http://mathworld.wolfram.com/SpherePointPicking.html

    A harder riddle is – how do you generate a uniform distribution in a ball volume? :)

    • Steven Pigeon says:

      My immediate idea would be to use a variation of the method I presented in the post, but generalizing it to spherical shells (depending on radius) rather than circumferences (depending on lattitudes).

      • Simon Read says:

        Easy!!!

        1) Generate point in cube – reject points outside sphere.

        2) Using some mathematics:

        First, generate a point on the surface of a sphere, THEN choose a radius.

        Volume relates to radius

        volume (of a shell) ~ radius^2 . thickness (factor involving pi)

        dV = r^2 dr

        ergo, integrate to get:

        V = r^3 (with factor, etc.)

        So… you want something which is uniformly distributed in volume… so

        r = V ^ (1/3) with V distributed uniformly

        or

        radius = ( U(0,1) ) ^(1/3)

        • The rejection method is very bad. The ratio of the volume of the sphere to the volume of the bounding (hyper)cube goes to zero somewhat fast as the number of dimensions increases (see comments below). In other words, as you increase the number of dimensions, the number of retries increases as well.

          As for 2), I’m not sure that’s the way I would derive the equations.

  3. Ofek Shilon says:

    I’m having trouble with the HTML comment formatting – I wish comments had a ‘preview’ feature. here goes another attempt

    While this may be possible, it would be very hard to generalize to an arbitrary dimensions. When I needed a solution some years ago, I came across Roger Stafford’s incredible solution. Here are some (too few) words of explanation.

    Note that Matlab’s gammainc function is commonly referred to as the regularized-lower-incomplete gamma function, so the code may seem misleading at first sight.

    • Steven Pigeon says:

      Would’nt a method based on rejection be efficient? Do the corners of a cube bounding a sphere grow or shrink as we add dimensions?

      • Ofek Shilon says:

        The ratio n-ball-volume/n-bounding-cube-volume goes to zero as n approaches infinity, so rejection grows more and more inefficient.

        • Steven Pigeon says:

          Quite so. The volume of a unit hypersphere is given by:

          \displaystyle{V_n=\frac{\pi^\frac{n}{2}}{\Gamma(\frac{n}{2}+1)}}

          while the volume of the bounding hypercube is H_n=2^n. And so

          \displaystyle{\lim_{n\to\infty} \frac{V_n}{H_n}=0}

          (the volume of a hypersphere of radius r is

          \displaystyle{V_n=\frac{\pi^\frac{n}{2}r^n}{\Gamma(\frac{n}{2}+1)}}

          )

  4. Ian Cutress says:

    I believe there are a couple of issues with your surface of resolution mathematics. The integral should reduce to 1.

  5. […] 2D points uniformly and randomly distributed a triangle, discussing the method of rejection. In a third post, I’ve discussed how to generate points on a […]

  6. […] series, I’ve discussed permutations, generating uniform random points in a triangle, on a sphere, the inversion method, even recycling expensive random bits. In this entry, I will use the […]

  7. Simon Read says:

    You put theta = U(-1,1)

    Surely it should be (following your above analysis)
    theta = arcsin ( U(-1,1) )

    ??

    There is an easier way to derive the above.

    A = area of sphere.

    dA= r^2 cos(theta) d.theta d.phi

    then, immediately:

    A = r^2 sin(theta) d.phi

    (phi is the longitude)

    so theta = arcsin A, maybe with a scaling factor to put it in the range -1 to 1.

    BUT WAIT…

    I don’t like the arcsin, as you lose some resolution, i.e. d(arcsin(x))/d(x) goes to infinity as x reaches +-1.

    That means the values of theta generated will fall into certain values with gaps between.

    There’s a totally different way of generating points on the surface of a sphere. This would also work with higher-dimensional spheres.

    Start with a vector of a point on the surface, and rotate it randomly.
    Just choose a matrix:
    [ cos(theta) sin(theta) 0 ]
    [ -sin(theta) cos(theta) 0 ]
    [ 0 0 1 ]

    Where theta is uniformly distributed in the range -pi to pi. In fact, the way this works, it can be uniformly distributed in ANY range, but a small range means you have to do more work before the coordinates are sufficiently random.

    1) Swap the coordinates of the vector around randomly (effectively some reflections in planes going through the origin)
    2) New vector = matrix . old vector

    Do this a couple of times and your point (vector) is randomly distributed on the surface of the sphere.

    There are variations, involving randomly choosing which elements of the matrix to put the cos(theta) and sin(theta) in, ensuring that it’s still a rotation matrix.

    Rounding error might give points marginally off the surface, so you can normalise the vector if you want.

    Essentially, the above generates a series of vectors, each one a random rotation of the previous vector. The individual transformations won’t necessarily be uniformly distributed, but the series will.

    This will free you from dependence on latitude and longitude, not particularly nice coordinates as they don’t behave gracefully, i.e.
    1) they have a finite range and just stop
    2) a small movement of the point occasionally means a massive change in theta and phi.

  8. […] another post, I gave the formula for the unit […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: