1

I am looking for an algorithm to uniformly sample the surface of a unit sphere.

I have been reading in https://mathworld.wolfram.com/SpherePointPicking.html.

Their approach appears to translate as:

def randomPoints(n):
    u = np.random.uniform(0, 1, size=n)
    v = np.random.uniform(0, 1, size=n)
    theta = 2*np.pi*u
    phi = np.arccos(2*v - 1)
    return theta, phi

Where theta = [0, 2pi) and phi = [0, pi]

I would like to know if anyone can tell me if this is the best approach to achieve a uniform distribution around the surface of a unit sphere. Or if there is some other algorithm which will reach an isotropic state faster.

EDIT: I have also been trying with a method by George Marsaglia 1972, The Annals of Mathematical Statistics, vol. 43, No. 2, 645 - 646.

r = np.sqrt(np.random.uniform(0.0, 1.0, n))
t = np.random.uniform(0, 2*np.pi, n)
v1, v2 = r*np.cos(t), r*np.sin(t)
s = v1*v1 + v2*v2 

a = 2*v1*(np.sqrt(1-s))
b = 2*v2*(np.sqrt(1-s))
c = 1-2*s

theta = np.arctan2(b, a) 
phi = np.arccos(c/(np.sqrt(a**2 + b**2 + c**2)))

Best Regards

mystic snake
  • 59
  • 5

1 Answers1

2

First I explain the algorithm: The function f(ϕ, θ)= sin(ϕ)/4π is the pdf of the uniform distribution on the unit sphere (because the area element on the sphere is dσ=sin(ϕ)dϕdθ, and the sphere area is 4π). The marginal densities are g(ϕ)=sin(ϕ)/2, and h(θ)=1/2π, i.e. the two random variables, ϕ and θ, are independent, and θ is uniform distributed on [0,2π]. The cummulative distribution function for ϕ is G(ϕ)=(1-cos(ϕ))/2, and its inverse G^(-1)(u)=acos(1-2u). Hence for sampling uniformly n points on the unit sphere we define:

import numpy as np
import plotly.graph_objects as go
from numpy import sin, cos, pi, arccos, sqrt

n= 2000
u = np.random.rand(n) 
phi = arccos(1 - 2*u) #inverse transform sampling 
theta = 2*pi*np.random.rand(n) #uniform sampling 
x=cos(theta)*sin(phi)
y=sin(theta)*sin(phi)
z=cos(phi)

fig=go.Figure(go.Scatter3d(x=x, y=y, z=z, mode="markers",marker_size=2, marker_color="red"))
fig.update_layout(width=500, height=500)

Another algorithm:

n=2000
z1  = np.random.normal(size=n)
z2=np.random.normal(size=n)
z3=np.random.normal(size=n)
r=sqrt(z1**2+z2**2+z3**2)
x, y, z=z1/r, z2/r, z3/r
fig=go.Figure(go.Scatter3d(x=x, y=y, z=z, mode="markers",marker_size=2, marker_color="red"))
fig.update_layout(width=500, height=500)

The two samplings are illustrated here (at left with 1st algorithm, at right with the 2nd one). enter image description here From my experience I cannot say that one is better than another. You should experiment for your use case and decide which is the most suitable to achieve the task

xecafe
  • 740
  • 6
  • 12