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

Ответ 1

На основе последнего подхода на этой странице вы можете просто сгенерировать вектор, состоящий из независимых выборок из трех стандартных нормальных распределений, а затем нормализовать вектор такой, что его величина равна 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) или на поверхностях гипермерных больших размерных гиперсфер.

Ответ 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 имеет фиксированное значение.

Ответ 3

После некоторого обсуждения с @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 xrange(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 xrange(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

Ответ 4

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

Выберите 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

Ответ 5

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

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

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

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

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

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

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

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

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