Ways to generate random points on a sphere

In [1]:
#All necesary import
import sys
import time

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline

sys.path.append('../')

from ten.utils import generate_random_points_in_sphere
In [3]:
n = 3000
R = 5

theta = np.random.uniform(low = 0, high = 2*np.pi, size = n)
phi = np.random.uniform(low = 0, high = np.pi, size = n)

x = np.sin(phi)*np.cos(theta)
y = np.sin(phi)*np.sin(theta)
z = np.cos(phi)

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z)
Out[3]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa775223ba8>

Note the large number of dots in the poles, there is no uniform distriuciĆ³n.

The right way

In [5]:
n = 3000
theta = 2*np.pi*np.random.uniform(size=n)
phi = np.arccos(2*np.random.uniform(size=n) - 1)

x = np.sin(phi)*np.cos(theta)
y = np.sin(phi)*np.sin(theta)
z = np.cos(phi)

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z)
Out[5]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa7750cf4a8>

Now we have a uniform distribution.But, what if we want to generate inside the sphere?

Inside the sphere

The wrong way

The above code, we add that the radius is variable.

In [6]:
n = 3000
theta = 2*np.pi*np.random.uniform(size=n)
phi = np.arccos(2*np.random.uniform(size=n) - 1)

#theta = np.random.uniform(low=0, high=2*np.pi, size=n)
#phi = np.arcsin(2 * np.random.uniform(low=0, high=2*np.pi, size=n) - 1)
u = np.random.uniform(low=0, high=10, size=n)

x = np.sin(phi)*np.cos(theta)*u
y = np.sin(phi)*np.sin(theta)*u
z = np.cos(phi)*u

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z)
Out[6]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa77510c4a8>

Obviously there is not a uniform distribution

The trivial way, with a if...

In [7]:
n = 3000 
R = 5
a = 0

x_l = []
y_l = []
z_l = []

t = time.time()

while a <= n:
    x = np.random.uniform(low=-R, high=R, size=1)
    y = np.random.uniform(low=-R, high=R, size=1)
    z = np.random.uniform(low=-R, high=R, size=1)
    if x*x + y*y + z*z <= R*R:
        a += 1
        x_l.append(x)
        y_l.append(y)
        z_l.append(z)

print("time =",time.time() - t)

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x_l, y_l, z_l)
time = 0.30646824836730957
Out[7]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa774f50780>