Uniform distribution of points on the surface of a sphere?

Question

I am looking for an algorithm to optimally distribute N points on the surface of a sphere (not randomly, I speak about optimal distribution in the sense of the Tammes problem). I want to do that in C++, but the problem is not a language problem, it is an algorithmic problem. Here is the algorithm I currently use (in python):

# Function to distribute N points on the surface of a sphere 
# (source: http://www.softimageblog.com/archives/115)
def uniform_spherical_distribution(N): 
    pts = []   
    inc = math.pi * (3 - math.sqrt(5)) 
    off = 2 / float(N) 
    for k in range(0, int(N)): 
        y = k * off - 1 + (off / 2) 
        r = math.sqrt(1 - y*y) 
        phi = k * inc 
        pts.append([math.cos(phi)*r, y, math.sin(phi)*r])   
    return pts

For a large number of points the visual result seems quite good but for a low number of points (between 2 and 10), the result does not seem good at all (more points in a hemisphere). Is there a way to enhance the algorithm to be also ok for a low number of points... (the answers can be in C, C++ or Python).

Answer 1

Brute force probably works if all you need is the low number of points. By brute force I mean precomputing the solutions for 2..10 and looking them up in an array if N<=10.

For small Ns like that, it should be possible to solve the mathematical optimization problem distributing your points optimally. You need multistart but for small Ns it should not be a problem. Here is my attempt in AMPL:

param N;

var x{1..N};
var y{1..N};
var z{1..N};

var d;

param xbest{1..N};
param ybest{1..N};
param zbest{1..N};
param dbest;

maximize obj: d;

all_points_on_the_sphere{i in 1..N}:
  x[i]^2 + y[i]^2 + z[i]^2 = 1;

all_pairs{i in 1..N, j in 1..N : i<j}:
  (x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2 >= d;

fix_first_x:  x[1] = 1;
fix_first_y:  y[1] = 0;
fix_first_z:  z[1] = 0;  
fix_second_z: z[2] = 0;

#############################################

option show_stats 1;
option presolve 10;
option substout 1;
option var_bounds 2;
#option nl_comments 0;
#option nl_permute 0;
option display_precision 0;

option solver "/home/ali/ampl/ipopt";

for {k in 2..10} {
    let N := k;  
    solexpand _con;
    let dbest := -1.0;
    # Multistart
    for {1..2000} {
        for {j in 1.._snvars}
            let _svar[j] := Uniform(-1, 1);
        let d := N;     

        solve;

        if (solve_result_num < 200 and d > dbest) then {
            let dbest := d;
            for {i in 1..N} {
                let xbest[i] := x[i];
                let ybest[i] := y[i];
                let zbest[i] := z[i];
            }
        }
    }

    print "@@@";
    printf "N=%d, min distance %6f\n", N, sqrt(dbest);

    for {i in 1..N}
        printf "(%9.6f, %9.6f, %9.6f)\n", xbest[i], ybest[i], zbest[i];
}

It took 5 minutes to run this script on my machine and the solutions are:

N=2, min distance 2.000000
( 1.000000,  0.000000,  0.000000)
(-1.000000,  0.000000,  0.000000)

N=3, min distance 1.732051
( 1.000000,  0.000000,  0.000000)
(-0.500000,  0.866025,  0.000000)
(-0.500000, -0.866025,  0.000000)

N=4, min distance 1.632993
( 1.000000,  0.000000,  0.000000)
(-0.333333, -0.942809,  0.000000)
(-0.333333,  0.471405, -0.816497)
(-0.333333,  0.471405,  0.816497)

N=5, min distance 1.414214
( 1.000000,  0.000000,  0.000000)
(-0.208840,  0.977950,  0.000000)
(-0.000000,  0.000000,  1.000000)
(-0.212683, -0.977121,  0.000000)
( 0.000000,  0.000000, -1.000000)

N=6, min distance 1.414214
( 1.000000,  0.000000,  0.000000)
(-1.000000,  0.000000,  0.000000)
( 0.000000, -0.752754, -0.658302)
( 0.000000,  0.752754,  0.658302)
( 0.000000,  0.658302, -0.752754)
( 0.000000, -0.658302,  0.752754)

N=7, min distance 1.256870
( 1.000000,  0.000000,  0.000000)
(-0.688059, -0.725655,  0.000000)
( 0.210138, -0.488836, -0.846689)
( 0.210138, -0.488836,  0.846688)
(-0.688059,  0.362827,  0.628435)
(-0.688059,  0.362827, -0.628435)
( 0.210138,  0.977672,  0.000000)

N=8, min distance 1.215563
( 1.000000,  0.000000,  0.000000)
( 0.261204, -0.965284,  0.000000)
( 0.261204,  0.565450,  0.782329)
(-0.783612, -0.482642, -0.391165)
( 0.261204, -0.199917, -0.944355)
( 0.261204,  0.882475, -0.391165)
(-0.783612,  0.599750,  0.162026)
(-0.477592, -0.399834,  0.782329)

N=9, min distance 1.154701
( 1.000000,  0.000000,  0.000000)
(-0.500000, -0.866025,  0.000000)
( 0.333333, -0.577350, -0.745356)
(-0.500000,  0.866025,  0.000000)
(-0.666667, -0.000000,  0.745356)
(-0.666667,  0.000000, -0.745356)
( 0.333333, -0.577350,  0.745356)
( 0.333333,  0.577350, -0.745356)
( 0.333333,  0.577350,  0.745356)

N=10, min distance 1.091426
( 1.000000,  0.000000,  0.000000)
(-0.605995,  0.795469,  0.000000)
( 0.404394,  0.816443,  0.412172)
(-0.664045, -0.746251, -0.046407)
( 0.404394, -0.363508, -0.839242)
(-0.664045,  0.002497,  0.747688)
(-0.605995,  0.046721, -0.794096)
( 0.404394, -0.908368,  0.106452)
( 0.255933,  0.703344, -0.663179)
( 0.404394, -0.159620,  0.900548)

By looking at the numbers, it is obvious that some of the solutions could be computed analytically (I recognize sqrt(2) and sqrt(3), etc). I believe for N=2, 4 and 6 the solutions are straight line ([-1,0,0], [1,0,0]), tetrahedron, octahedron.

There is no strong guarantee that the above are the best possible distribution of the points. The nonlinear solver could have gotten stuck in a local optimum; as N grows, the number of local optima also grows.

You can just put the above solutions in an array and use them in Python, C++ or whatever language it is that you use.