
Random Points on a Sphere - santaclaus
https://www.jasondavies.com/maps/random-points/
======
davmre
The easiest way to generate random points on a (hyper)sphere, at least with
respect to most sampling APIs, is to sample from an isotropic normal
distribution of the appropriate dimension (for a three-dimensional sphere,
just call randn() three times), then normalize the result to be on the surface
of the sphere.

This works because the normal distribution is rotationally symmetric, so each
sample is essentially a random unit vector scaled by a random (normally
distributed) radius. Once you remove the variation in radius you just get back
the random unit vector, i.e. a point on the unit sphere.

~~~
josephg
This is not quite correct. Generating a random vector using three calls to
rand() is the equivalent of picking a random point inside a cube. If you
normalize the points, you're effectively squishing the corners of the cube in
toward the center to make a sphere. Points near the corners of the containing
cube will appear more frequently than the top, bottom and 'sides'.

@ginko's comment is correct - you can fix the algorithm by throwing out any
points that lie outside the sphere before normalizing.

~~~
davmre
As others have pointed out, the distinction between randn() and rand() is
crucial here - the former gives you points from a (spherically symmetric)
normal distribution, the latter gives you uniform points from the unit cube.

One of the advantages of the normal-sampling route over @ginko's rejection-
based method is that in high dimensions _almost all_ of the volume of the unit
cube is situated outside of the unit sphere (the volume of the unit cube is
always 1, whereas the volume of the unit sphere decreases exponentially with
dimension). So the rejection method becomes exponentially slow in high
dimensions, while the Gaussian method still works just fine.

~~~
Retr0spectrum
Surely there would still be a bias towards the corners?

~~~
darkmighty
Interestingly, there isn't for Gaussians (see my comment below).

Other simple distributions tend to give biases towards the corners or axis.
Perhaps the Gaussian is unique in this regard? I'm not sure.

~~~
OscarCunningham
Yes, the fact that i.i.d. Gaussians give a rotationally symmetric joint
distribution is a defining property of the Gaussian.

~~~
darkmighty
Ah I suspected as much. Can you provide proof?

I think I found one, but I'm not sure:

Without loss of generality, take f(xi) = k1 * exp(-g(xi)) [1], for some g.
Then we need the joint pdf to satisfy f(x1,...,xn) = k2 * exp(-h(R^2)),
R=sum(xi^2)^1/2 (the R^2 and h(.) is w.l.g. too). So we get
g(x1)+g(x2)=h(x1^2+x2^2). Then assuming the functions g and h analytic we end
up needing g(x)= k * x^2, otherwise we get cross terms in the Taylor expansion
that can't be cancelled out for all xi. Sounds good?

[1] The function f trivially needs to be symmetric, justifying no loss of
generality.

~~~
madcaptenor
Perhaps coincidentally, someone just asked about this at math.stackexchange:
[http://math.stackexchange.com/questions/1255637/joint-pdf-
of...](http://math.stackexchange.com/questions/1255637/joint-pdf-of-n-1-i-i-d-
random-variables-isotropic-if-and-only-if-they-are-cen)

~~~
darkmighty
Yes, that was me :)

I wanted to make sure I got a proof, since I didn't really find this
elsewhere.

------
dietrichepp
> Note that the points are no longer independent of each other, hence they are
> no longer uniformly distributed.

This is false. Random variables can be uniform even if they are not
independent.

For example, let X be uniform(0, 1), and let Y = X + 0.1 mod 1. Clearly they
are not independent. Clearly X is uniform. When you realize that you could
redefine the variables as Y = uniform(0, 1) and X = Y - 0.1 mod 1, then it
becomes obvious that Y is also uniform.

It is a common problem to know a variable's distribution without knowing
whether it is independent.

~~~
davmre
Good catch. Probably what the article meant to say is that the lack of
independence means the points are not _conditionally_ uniform. That is,
although X and Y follow uniform distributions individually, the dependence
between them means that once you observe X, your resulting belief about Y is
_no longer_ a uniform distribution. Compare this to the independent case where
knowing X would tell you nothing about Y, so your belief about Y conditioned
on X would still just be uniform.

------
sjolsen
>Although we’ve successfully generated uniformly distributed points on a
sphere, it feels messy. Some points seem too close together, and some seem too
far apart.

This is how randomness _should_ look. If you randomly generate (using a high-
quality random engine) uniformly distributed points on a square plane or even
a line, you'll see the same apparent clustering.

Put another way, if you plot every point on a coordinate grid with integer
coordinates, you'll have a uniform distribution, but it will hardly be
"random."

------
ggchappell
FTA:

> One solution is to pick λ ∈ [-180°, 180°) as before and then set φ =
> cos^(-1)(2x - 1), where x is uniformly distributed and x ∈ [0, 1).

There is a fascinating fact hidden in the above sentence: if a point (x,y,z)
is uniformly distributed on the surface of a sphere -- say, of radius 1,
centered at (0,0,0) -- then each of x, y, and z is uniformly distributed over
the interval [-1, 1].

Honestly, that's amazing.

~~~
fsloth
Are you sure? If z is distributed over the surface of the sphere there is far
less surface at the extremums -1 and 1. Doesn't z being uniform in the range
(-1,1) lead to points being actually more densly concentrated near the poles?
Or did I misunderstand what you meant?

~~~
kasdkashdkas
Near the poles the surface becomes very flat and therefore there is as much
surface area in a slice of constant height than anywhere else. Archimedes was
a clever fella :-)

~~~
fsloth
Oh, I was totally confused, of course the numbers in every dimension are
evenly distributed, silly me - I was thinking of great circles. Thanks :)

------
ginko
A very simple solution is to generate a 3D point in [-1,1]^3, regenerate if
it's length is >1 and project it on to the sphere otherwise.

~~~
reitzensteinm
Why would you regenerate it if the length exceeds 1? Can't you just normalize
the vector regardless of its length (excluding 0)?

~~~
dietrichepp
That would lead to the distribution being non-uniform: there would be more
points in the corners of the cube.

~~~
reitzensteinm
Ah, of course. To visualize it, the distance between the center and the
surface of the shape would be proportional to the probability of that vector
being selected after normalization. Thanks!

------
kayamon
The best and simplest results can be obtained by using a Hammersley Point
Set[1].

[1]
[https://www.cse.cuhk.edu.hk/~ttwong/papers/udpoint/udpoint.p...](https://www.cse.cuhk.edu.hk/~ttwong/papers/udpoint/udpoint.pdf)

------
thaumasiotes
related:
[http://en.wikipedia.org/wiki/Bertrand_paradox_%28probability...](http://en.wikipedia.org/wiki/Bertrand_paradox_%28probability%29)

("choose a random chord intersecting a circle")

If someone is familiar with this kind of thing, I'd like to know a little more
about the Jaynes solution described. It seems to me that if the circle were
translated off to some distant portion of the plane, it would be very unlikely
that any method of choosing chords would have the same distributional
properties on the original circle and on the translated circle, because the
chords are all originally restricted to intersect the first circle and, given
two circles in a plane, there are many lines that intersect one but not both
circles. Will Jaynes' method really shade the translated circle the same way
it will shade the original? My intuition says the translated one should be
shaded more heavily in a sort of "band" parallel to the line connecting the
origins of the two circles, or possibly in a crescent pattern with more
shading closer to the original.

~~~
kgwgk
It's easy to see that the second method can be applied so it works
independently of the position and size of the circle. Let's start by
reformulating it slightly: instead of a "random radius" we will use a "random
diameter".

> Choose a diameter of the circle, choose a point on the diameter and
> construct the chord through this point and perpendicular to the diameter.

In the general case, choose a line going through the origin.

=> In the circle, there is a diameter parallel to this line.

Then choose a point on the line, and construct the perpendicular line at that
point. (To be able to choose a point uniformly, you have to consider a finite
segment but this is not a problem as long as the location/size of the circle
is bounded).

=> This perpendicular line might not intersect the circle, but when it does it
determines a chord that goes through a point chosen uniformly on the diameter.

This generalized method can be used to cover uniformly multiple circles at the
same time, but obviously not every line will intersect every circle.

------
cturhan
If you have random N points and you want to distribute them uniformly on a
sphere, this method won't work. I needed this in two different side projects
and used this[1] which worked well.

[http://dx.doi.org/10.1016/j.jocs.2011.06.007](http://dx.doi.org/10.1016/j.jocs.2011.06.007)

------
krazydad
I stumbled upon an attractive & simple to implement non-random method that
uses the golden angle (an angle commonly observed in plant phyllotaxy) --
drawing a fibonacci spiral on the sphere while maintaining a constant surface
area. It is demonstrated here, on the open processing website.

[http://www.openprocessing.org/sketch/41142](http://www.openprocessing.org/sketch/41142)

    
    
      n = number of points
      phi = (sqrt(5)+1)/2 - 1
      ga = phi * 2 * PI
    
      for each point i (1..n)
        longitude = ga*i
        latitude = asin(-1 + 2*i/n)

------
logicallee
am I the only one who thought the points on the second sphere look _way_ more
randomly distributed than the points on the third sphere?

~~~
spacemanmatt
They are not random on the third sphere.

~~~
logicallee
so what's the point? WHo wants a "random" sequence like
0100101101011011011001001001010110111010010101010 - that's not random at all!

~~~
kgwgk
Sometimes you want "quasi-random" sequences:

[http://en.wikipedia.org/wiki/Low-
discrepancy_sequence](http://en.wikipedia.org/wiki/Low-discrepancy_sequence)

------
cafebeen
The statistical way of approaching this problem is to sample from a Von-Mises
Fisher distribution with k=0:

[https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distr...](https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution)

------
tantalor
The div.sphere div.wrap elements on this page are lovely, and they are created
by the client as well!

------
starmole
This post seems to barely scratch the surface of the problem. First you really
need to specify by what metric you want the distribution to be uniform.
Euclidean 3d distance? Arcs? The way I would do it would be to generate a lot
of random points by some easy scheme, like long/lat, then feed them into a
kdtree (k=3), find all neighbors for every point with an overestimating
metric, then check them against the real metric and remove the too close ones.

~~~
j2kun
It should be clear that the obvious notion of uniformity is the intended one.
I.e., it should be uniform with respect to the usual Lebesgue measure on the
sphere.

And your method is extremely complicated. You can solve this problem with 5
lines of code and no data structures.

