Spreading points on a disc and on a sphere

When working with computer graphics, machine vision or simulations of various kind, two problems (among many others) which keep popping out are

• How to have N points approximately evenly distributed over a disc
• How to have N points approximately evenly distributed over a sphere

One way to build a solution for those two problems is to rely on a particle system. Each point is a particle, particles are repulsing each other with a $\frac{1}{r^2}$ force. Put the particles on the disc or the sphere surface, at random locations. Shake vigorously the particles, add some small drag force, let the particles reach a steady state, shake again ... rinse, repeat. It is not as easy as it sounds. How to setup the forces ? How to handle the boundary of the disc ? Which integrator to use to compute the particle's motion ? And for more than a few hundred particles, it will be terribly slow.

Spirals

There are much simpler and leaner algorithms to evenly distribute N points over a disc or a sphere, if one can live with just a good approximation. The main idea fits in one word : spiral ! Let's start with a disc. Imagine a spiral of dots starting from the center of the disc. In polar coordinates, the N points are produced by the sequence

$${\rho }_i=\theta i$$

$${\tau }_i=\sqrt{\frac{i}{N}}$$

${\rho }_i$ and ${\tau }_i$ are respectively the angle in radian and the radius of the i-th point. Why ${\tau }_i=\sqrt{\frac{i}{N}}$ ? Let's try to cut a unit disc in one disc and one ring of equal areas, by cutting the unit disc at radius $r_c$. The small disc and the ring are of equal areas, thus $\pi r_1^2=\pi -\pi r_1^2$. Thus $r_1=\sqrt{\frac{1}{2}}$. Now, let's try to cut the unit disc in one disc and 2 concentric rings, all of equal areas. A slightly more complex calculation (a system of 2 linear equations) would tell us that we should cut the unit disc at radius $r_1=\sqrt{\frac{1}{3}}$ and $r_2=\sqrt{\frac{2}{3}}$. The calculation is a more complex when cutting the unit disc in N equal areas rings, but you can guess the answer now : $r_1=\sqrt{\frac{1}{N}}$, $r_2=\sqrt{\frac{2}{N}}$, $r_i=\sqrt{\frac{i}{N}}$.

What about the ideal $\theta$ ? Brace yourself... This is the golden angle, $\pi (3-\sqrt{5})$ ! It's roughly equal to 137.508°, or about 2.39996 radians. The golden angle is the "most irrational" angle, defined as $2\pi (1-\frac{1}{\varphi })$ with $\varphi$ being the golden ratio. If $\theta$ is a rational number, we would obtain clusters of points aligned with the center of the disk. Thus $\theta$ have to be irrational. As it turns out, the golden ratio is the irrational number the hardest approximate with a continued fraction. Written as a continued fraction, the golden ratio is the irrational number with the slowest convergence of all the irrational numbers. The golden angle gives the best possible spread for the ${\rho }_i$ angles. This method to spread points over a disc is called Vogel's method.

Vogel's method is dead simple to implement. In pure, barebone Python, it can be done like this

import math

n = 256

golden_angle = math.pi * (3 - math.sqrt(5))

points = []
for i in xrange(n):
  theta = i * golden_angle
  r = math.sqrt(i) / math.sqrt(n)
  points.append((r * math.cos(theta), r * math.sin(theta)))

Using numpy, one can use vector operations like this

import numpy

n = 256

radius = numpy.sqrt(numpy.arange(n) / float(n))

golden_angle = numpy.pi * (3 - numpy.sqrt(5))
theta = golden_angle * numpy.arange(n)

points = radius[:,None] * numpy.stack([numpy.cos(theta), numpy.sin(theta)])

Spheres

What about the sphere ? We can reuse the golden angle spiral trick. In cylindrical coordinates, we would generate N points like this

$${\rho }_i=\theta i$$ $${\tau }_i=\sqrt{1-z_i^2}$$ $$z_i=(1-\frac{1}{N})(1-\frac{2i}{N-1})$$

And voila, an extremely cheap method to spread point evenly over a sphere's surface ! The code for this method is very close to the one for Vogel's method.

n = 256

golden_angle = numpy.pi * (3 - numpy.sqrt(5))
theta = golden_angle * numpy.arange(n)

z = numpy.linspace(1. - 1. / n, 1. / n - 1., n)
radius = numpy.sqrt(1. - z ** 2)

points = numpy.stack([radius * numpy.cos(theta),
                      radius * numpy.sin(theta),
                      z], axis = 1)

Note that you can also use this method to make a sphere mesh. The sphere will look good even with few polygons, but the indexing and computing the normals won't be especially fast. For this specific problem, making a sphere mesh, Iñigo Quilez have a very nice, efficient method introduced here.

I made ready-to-use Python functions on Github for the disc and the sphere.