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

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$ 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 bit 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}{\phi})$ with $\phi$ 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. Effect of the angle step parameter, on 256 points. The 2 leftmost points layout are for rational values of the angle step. The top center one is for square root of 2, the 3 others are for others irrational values.

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

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

points = numpy.zeros((n, 2))
points[:,0] = numpy.cos(theta)
points[:,1] = numpy.sin(theta)


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

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.0 / n, 1.0 / n - 1, n)
radius = numpy.sqrt(1 - z * z)

points = numpy.zeros((n, 3))