# Choosing evenly distributed points from a million points on a sphere

I have $$\10^5\$$ to $$\10^6\$$ points on a sphere, and want to choose some points from them which are as close as uniformly distributed as possible. For that reason, I do the following: at each step I add the point in the data which is farthest away from any point which I have already chosen. In order to do this, I find the closest distance from each point in the data to any of the points that are already chosen. Then I pick the point from the data which has the largest minimum distance to any of the points that are already chosen. The algorithm is the following:

1. Choose a random point from all the $$\10^6\$$ points, add it to the chosen set of points

2. For each of the $$\10^6\$$ points:

2.1 Find the distances to all the points which are chosen already

2.2 Find the min distance to a point which is chosen already

2.3 Find the point which is furthest away from any point which is chosen already

2.4 Add this point to the chosen points

3. Repeat 2 until I have enough points (say 1000 or 10000 or 100000).

data is $$\10^6\$$ by 3 array of points on a sphere, e.g.

[[ 0.26750522 -0.92342735  0.27517791]
[-0.26753053  0.9228284  -0.27715548]
[ 0.34058837  0.35873472  0.86908513]
[-0.33563178 -0.35794416 -0.8713365 ]
[-0.8945523  -0.15682272  0.41854847]
[ 0.906739    0.15103321 -0.39371734]
[ 0.49138659 -0.64830133 -0.581588  ]
[-0.87922161  0.24492971 -0.40863039]
[ 0.18062012  0.39852148 -0.89919798]
[-0.49103509  0.65872966  0.57005243]
[-0.55615839 -0.82248645 -0.11918007]
[ 0.85802915 -0.25207255  0.44748789]
[-0.16990087 -0.40390349  0.89888579]]


I have implemented this in Python:

import random
import numpy

N = 1000

starting_point_ind = random.SystemRandom().randint (1 , len(data))

points_from_data = numpy.array(numpy.zeros((N, 3)), dtype=float, order='C')
points_from_data = data[starting_point_ind]

for i in range(1, N):
distancesToClosestPoints = numpy.array(numpy.zeros(len(data)), dtype=float)

for j, point in enumerate(data):
distancesToEstablishedPoints = numpy.linalg.norm (point - points_from_data[0:i] , axis=1)

distancesToClosestPoints[j] = distancesToEstablishedPoints.min()

k = distancesToClosestPoints.argmax()

points_from_data[i] = data[k]


This doesn't run very efficiently for large N, so I also tried to optimise it by doing

import random
from sklearn.metrics import pairwise_distances
import time

t0 = time.time()

N = 1000

starting_point_ind = random.SystemRandom().randint (1 , len(data))

points_from_data = numpy.array([data[starting_point_ind]])

for i in range(1, N):
pairwise_distances_to_data = pairwise_distances (data, Y=points_from_data, metric='euclidean', n_jobs=-1)
pairwise_distances_to_data = numpy.array(pairwise_distances_to_data)

min_distances_to_data = numpy.amin(pairwise_distances_to_data, axis=1)

k = min_distances_to_data.argmax()

points_from_data = numpy.append(points_from_data, [data[k]], axis=0)

print(points_from_data)


But this still has the issue that the pairwise distances are computed each time, and has obvious memory issues for large N. I would like to hear suggestions how I can improve the performance of the algorithm.

• You might be interested in this article. Make a mesh of necessary size, and pick one point per triangle.
– vnp
Jun 6, 2019 at 17:13
• @Gordon I think you've got a good question in your hands, you just need to clarify what your code really does and what you really want your code to do, then ask yourself if you've achieved your goal. Jun 6, 2019 at 18:32
• If I compute the norm of two of the vectors you gave in the data example, they are not the same length, so points in data aren't part of a sphere. Is the example data not good or is there something that doesn't work? Jun 6, 2019 at 18:49
• @IEatBagels I realise my description of what I am trying to achieve was loose at first, but I think our discussion narrowed it down. At any rate, the results produced by the 2 chunks of code are equivalent and are exactly what I need. I just want to know if I can do it in a better way and faster. The vectors in data should have a norm of around 1.0+/-eps. I think the discrepancy comes from not copying enough significant figures when I pasted it from Jupyter. You can easily generate 10^6 points on the sphere, though.
– John
Jun 6, 2019 at 18:53

Since we can assume all your points are on a sphere, there is a better way to compute distance than using the L2 Norm (Euclidean distance, which you're using).

If we transform your data points to Spherical coordinates, we could use a simple Manhattan distance on the theta and rho coordinates, which is less expensive to compute. Since you compute this distance about $$\Nm\$$ times (6 million points * 1000 sample points), I'm sure it would make a good performance difference.

Small note on the Spherical Coordinates :

Instead of having x,y,z coordinates, we use r, theta, rho, which correspond to :

• r : the length of the vector
• theta : angle on the XY plane
• phi : angle of inclination starting from theta

So, to change data to spherical coordinates, we'd do the following (I assume data is arranged in a x,y,z manner). The conversion is very fast, but it could be further optimized if needed :

def convert_to_spherical(data):
spherical_data = np.zeros_like(data)
for i,d in enumerate(data):
r = math.sqrt(d**2 + d**2 + d**2)
theta = math.acos(d / r)
phi = math.atan(d / d)

spherical_data[i] = np.array([r, theta, phi])

return spherical_data


Then, you could use this instead of data and change the metric used by pairwise_distances to manhattan. You'll need to dig a little bit to understand why this works, but consider that since the parameter r is always the same, you can simply compute the difference between the angles.

Since my laptop is pretty bad, I couldn't generate data to sample 1000 points out of 1000000. These times were generated sampling 300 points from 100000 points, I started timing after generating the points but I included the coordinates transformation in the times :

Time using your optimized solution : 213.1844208240509s
Time using the spherical coordinates : 140.83173203468323s


Considering this pretty big time difference, I'm pretty sure performance would be much better on your data sample.

If afterwards, you need to get your coordinates back to the cartesian plane, the formula is pretty easy and you can find it online easily.

I've generated the points using thanks to this answer :

def sample_spherical(npoints, ndim=3):
vec = numpy.random.randn(ndim, npoints)
vec /= numpy.linalg.norm(vec, axis=0)
return numpy.transpose(vec)

• Thank you for the answer. I have a question about using the Manhattan distances. Imagine I am comparing 2 points with ((\theta_1, \phi_1) = (0, 0.1)) and ((\theta_2, \phi_2) = (0, 2\pi - 0.1)). This would imply they are a distance (2\pi - 0.2) apart, when in reality they are a distance of 0.2 apart. Wouldn't this mess up the algorithm?
– John
Jun 6, 2019 at 21:01
• I'm pretty damn sure convert_to_spherical can also be computed using numpy instead of looping over the data rows as you do in your current implementation. Jun 6, 2019 at 21:58
• So how would I compare the distances then?
– John
Jun 7, 2019 at 0:12
• @Gordon I misread your question I'm sorry. I understand your problem and I didn't think about it at first, because if we're talking about radians 2*pi = , but the algorithm won't figure that out. You'd need to create a custom distance metric based on the manhattan distance that makes a modulo 2*pi. I'm pretty confident it'll still be faster, but I'm sorry I didn't think about it first. Jun 7, 2019 at 13:42
• @Gordon you could also check this out : en.wikipedia.org/wiki/Great-circle_distance Jun 7, 2019 at 13:44