Distributing Points on a Sphere

When working on algorithms dealing with real world objects, a common task is to view the object from a variety of points of view.  If we represent a point of view looking at the object as points on a sphere, we have the problem of distributing the points on the sphere evenly.  This prevents gaps (because points are too far apart in some places) and wasted computation (because points are too close in some places).

The naive approach would be to use something like latitude-latitude lines.  This is very simple:  rotate the point of view in longitude some number of degrees around the equator; increase latitude some number of degrees; repeat latitude.  There are a number of problems with this.  First, it assumes that the latitude and longitude degree increments add up to 360 degrees in longitude and +/- 90 degrees in latitude.  That means that the number of points on the sphere that can be used is very granular.  Second, and more importantly, the points are very unevenly distributed on the sphere. 

The points near the poles are very close to each other, and result in a great deal of uneven distribution and wasted computation.  

An alternative approach is to put the points on the sphere randomly and then adjust them to be more even.  A simple way to do this is treat them as connected to each other via springs and simulate their movement.  This is an interesting software engineering problem in itself, since you have to pick how many points each point connects to, what the mass of the points are, what the spring parameters are, friction or other non-linear forces, how to efficiently simulate, and then know when to stop.  In the most simple case, without friction, you can have waves of compression and expansion travel indefinitely over the sphere.  With too much friction or spring parameters too low, the simulation will take a very long time to settle down.  The biggest problem is that this whole approach is inefficient.

A clean, efficient solution, for any number of points, is to use a 3D version of a Golden Spiral.  In school, you may have learned about the distribution of seeds in a sunflower, how shells spiral out, and other natural growth phenomena form a series of Golden Spirals.  Here is a close up of the seeds in a sunflower:

The mathematical model for this uses Vogel’s method, discussed on this Golden Angle Stack Exchange page.  Basically, it is a logarithmic spiral with a very special value for the growth parameter.  Here are the equations, as a function of ‘n’, the point number in polar coordinates: and , where phi is the Golden Ratio.  And here is what they look like:

To get a better feel for how all this works, and the wonderful relationship to the Fibonacci series, watch these wonderful videos by Vihart:    Doodling in Math: Spirals, Fibonacci, and Being a Plant (Actually, it’s worthwhile seeing all the videos by Vihart.)  Conceptually, what is happening is that the angle between successive points is based on the Golden Angle, and the radius grows over time.  

If you watch the video, you will see that Vihart shows the spirals on a pine cone.  That’s actually pretty close to what we want for distributing points on a sphere.  The successive points on the sphere will still be based on the Golden Angle, so , as above with Vogel’s method.  In cartesian coordinates, x and y will be:   x = radius * cos(theta) and y = radius * sin(theta).   The z will be evenly distributed between the top of the sphere and the bottom of the sphere.  In order for this to sit on the unit sphere, this means that the radius value will be sqrt(1-z*z).  The resulting plot of points looks like this:

In the diagram, I’ve hidden the points on the back of the sphere.  If you go through the bother of looking at how far points are from their nearest neighbors, which can be used as a measure of how well distributed points are, you will see that they are nicely distributed.   

The (surprisingly simple) code in python looks like this:  

import math
import numpy as np

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

pts = numpy.zeros((n, 3))
for ii in range(0,n):
    z = (1 - (1/n)) - (2 * ii / n)
    r = math.sqrt(1 - z*z)

    d = gold_angle * ii
    x = r * math.cos(d)
    y = r * math.sin(d)
    pts[ii,0] = x
    pts[ii,1] = y
    pts[ii,2] = z

A reasonable question would be where the calculation for the golden_angle came from, since above it was defined as  2*pi / (phi+1).  If you plug in the formula for the golden ratio and simplify, you get pi * (3-sqrt(5)).  See Golden Angle .  By the way, I know the code is not very pythonic;  in practice, I’ve use linspace.