Сгенерируйте случайную выборку точек, распределенных на поверхности единичной сферы

Основываясь на последнем подходе на этой странице, вы можете просто сгенерировать вектор, состоящий из независимых выборок из трех стандартных нормальных распределений, затем нормализуйте вектор так, чтобы его величина была равна 1:

import numpy as np

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

Например:

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d

phi = np.linspace(0, np.pi, 20)
theta = np.linspace(0, 2 * np.pi, 40)
x = np.outer(np.sin(theta), np.cos(phi))
y = np.outer(np.sin(theta), np.sin(phi))
z = np.outer(np.cos(theta), np.ones_like(phi))

xi, yi, zi = sample_spherical(100)

fig, ax = plt.subplots(1, 1, subplot_kw={'projection':'3d', 'aspect':'equal'})
ax.plot_wireframe(x, y, z, color='k', rstride=1, cstride=1)
ax.scatter(xi, yi, zi, s=100, c='r', zorder=10)

Тот же метод также распространяется на выбор равномерно распределенных точек на единичной окружности (ndim=2) или на поверхностях единичных гиперсфер более высоких измерений.

Точки на поверхности сферы можно выразить с помощью двух сферических координат theta и phi с 0 < theta < 2pi и 0 < phi < pi.

Формула перевода в декартовы x, y, z координаты:

x = r * cos(theta) * sin(phi)
y = r * sin(theta) * sin(phi)
z = r * cos(phi)

где r - радиус сферы.

Таким образом, программа могла произвольно выбирать theta и phi в их диапазонах с равномерным распределением и генерировать из них декартовы координаты.

Но тогда точки распределяются более плотно на полюсах сферы. Чтобы точки были равномерно распределены на поверхности сферы, phi необходимо выбрать как phi = acos(a), где -1 < a < 1 выбрано для равномерного распределения.

Для кода Numpy он будет таким же, как в Выборка равномерно распределенных случайных точек внутри сферического объема , за исключением того, что переменная radius имеет фиксированное значение.

Другой способ, который в зависимости от оборудования может быть намного быстрее.

Выберите a, b, c, чтобы получить три случайных числа от -1 до 1 каждое.

Рассчитать r2 = a^2 + b^2 + c^2

Если r2> 1.0 (= точка не в сфере) или r2 ‹0,00001 (= точка слишком близко к центру, у нас будет деление на ноль при проецировании на поверхность сферы), вы отбрасываете значения и выберите другой набор случайных a, b, c

В противном случае у вас есть случайная точка (относительно центра сферы):

ir = R / sqrt(r2)
x = a * ir
y = b * ir
z = c * ir

После некоторого обсуждения с @Soonts мне стало интересно узнать о производительности трех подходов, используемых в ответах: один с генерацией случайных углов, один с использованием нормально распределенных координат и один, отклоняющий равномерно распределенные точки.

Вот моя попытка сравнения:

import numpy as np

def sample_trig(npoints):
    theta = 2*np.pi*np.random.rand(npoints)
    phi = np.arccos(2*np.random.rand(npoints)-1)
    x = np.cos(theta) * np.sin(phi)
    y = np.sin(theta) * np.sin(phi)
    z = np.cos(phi)
    return np.array([x,y,z])

def sample_normals(npoints):
    vec = np.random.randn(3, npoints)
    vec /= np.linalg.norm(vec, axis=0)
    return vec

def sample_reject(npoints):
    vec = np.zeros((3,npoints))
    abc = 2*np.random.rand(3,npoints)-1
    norms = np.linalg.norm(abc,axis=0) 
    mymask = norms<=1
    abc = abc[:,mymask]/norms[mymask]
    k = abc.shape[1]
    vec[:,0:k] = abc
    while k<npoints:
       abc = 2*np.random.rand(3)-1
       norm = np.linalg.norm(abc)
       if 1e-5 <= norm <= 1:  
           vec[:,k] = abc/norm
           k = k+1
    return vec

Затем за 1000 очков

In [449]: timeit sample_trig(1000)
1000 loops, best of 3: 236 µs per loop

In [450]: timeit sample_normals(1000)
10000 loops, best of 3: 172 µs per loop

In [451]: timeit sample_reject(1000)
100 loops, best of 3: 13.7 ms per loop

Обратите внимание, что в реализации, основанной на отклонении, я сначала сгенерировал npoints выборки и выбросил плохие, а для генерации остальных точек я использовал только цикл. Казалось, что прямое пошаговое отклонение занимает больше времени. Я также удалил проверку деления на ноль, чтобы было чище сравнение с случаем sample_normals.

Удаление векторизации из двух прямых методов ставит их в один ряд:

def sample_trig_loop(npoints):
    x = np.zeros(npoints)
    y = np.zeros(npoints)
    z = np.zeros(npoints)
    for k in range(npoints):
        theta = 2*np.pi*np.random.rand()
        phi = np.arccos(2*np.random.rand()-1)
        x[k] = np.cos(theta) * np.sin(phi)
        y[k] = np.sin(theta) * np.sin(phi)
        z[k] = np.cos(phi)
    return np.array([x,y,z])

def sample_normals_loop(npoints):
    vec = np.zeros((3,npoints))
    for k in range(npoints):
      tvec = np.random.randn(3)
      vec[:,k] = tvec/np.linalg.norm(tvec)
    return vec

In [464]: timeit sample_trig(1000)
1000 loops, best of 3: 236 µs per loop

In [465]: timeit sample_normals(1000)
10000 loops, best of 3: 173 µs per loop

In [466]: timeit sample_reject(1000)
100 loops, best of 3: 14 ms per loop

In [467]: timeit sample_trig_loop(1000)
100 loops, best of 3: 7.92 ms per loop

In [468]: timeit sample_normals_loop(1000)
100 loops, best of 3: 10.9 ms per loop

(отредактировано с учетом исправлений, внесенных в комментарии)

В 2004 году я исследовал несколько подходов к этой проблеме с использованием постоянного времени.

предполагая, что вы работаете в сферических координатах, где theta - угол вокруг вертикальной оси (например, долгота), а phi - угол, поднятый вверх от экватора (например, широта), тогда для получения равномерного распределения случайных точек в полушарии к северу от по экватору вы делаете это:

1. выберите theta = rand (0, 360).
2. выберите phi = 90 * (1 - sqrt (rand (0, 1))).

чтобы получить точки на сфере, а не на полусфере, тогда просто инвертируйте phi в 50% случаев.

для любопытных, аналогичный подход имеет место для создания равномерно распределенных точек на единичном диске:

1. выберите theta = rand (0, 360).
2. выберите radius = sqrt (rand (0, 1)).

У меня нет доказательств правильности этих подходов, но я использовал их с большим успехом за последнее десятилетие или около того, и я убежден в их правильности.

некоторые иллюстрации (с 2004 г.) различных подходов находятся здесь, включая визуализация подхода к выбору точек на поверхности куба и их нормализации на сфере.