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}{\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.

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 = \left(1 - \frac{1}{N} \right) \left(1 - \frac{2i}{N - 1} \right) \]

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.