Erzeugen Sie eine zufällige Auswahl von Punkten, die auf der Oberfläche einer Einheitskugel verteilt sind

Basierend auf dem letzten Ansatz auf dieser Seite können Sie einfach einen Vektor aus unabhängigen Stichproben aus drei Standardnormalverteilungen generieren und den Vektor so normalisieren, dass seine Größe 1 ist:

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

Zum Beispiel:

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)

Dieselbe Methode verallgemeinert auch das Aufnehmen von gleichmäßig verteilten Punkten auf dem Einheitskreis (ndim=2) oder auf den Oberflächen von höherdimensionalen Einheitshyperkugeln.

Punkte auf der Oberfläche einer Kugel können mit zwei Kugelkoordinaten theta und phi mit 0 < theta < 2pi und 0 < phi < pi ausgedrückt werden.

Umrechnungsformel in kartesische x, y, z-Koordinaten:

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

wobei r der Radius der Kugel ist.

Das Programm könnte also theta und phi in gleichmäßigen Verteilungen zufällig abtasten und daraus die kartesischen Koordinaten generieren.

Dann werden die Punkte dichter auf den Polen der Kugel verteilt. Damit Punkte gleichmäßig auf der Kugeloberfläche verteilt werden, muss phi als phi = acos(a) ausgewählt werden, wobei -1 < a < 1 für eine gleichmäßige Verteilung ausgewählt wird.

Für den Numpy-Code wäre dies das Gleiche wie bei Sampling gleichmäßig verteilter zufälliger Punkte innerhalb eines Kugelvolumens , mit der Ausnahme, dass die Variable radius einen festen Wert hat.

Nach einigen Diskussionen mit @Soonts wurde ich neugierig auf die Leistung der drei Ansätze, die in den Antworten verwendet wurden: einer, bei dem zufällige Winkel erzeugt werden, einer, der normalverteilte Koordinaten verwendet, und einer, der gleichmäßig verteilte Punkte ablehnt.

Hier ist mein versuchter Vergleich:

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

Dann für 1000 Punkte

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

Beachten Sie, dass ich in der ablehnungsbasierten Implementierung zuerst npoints-Samples generiert und die schlechten Samples verworfen habe. Ich habe nur eine Schleife zum Generieren der restlichen Punkte verwendet. Es schien so zu sein, dass die direkte schrittweise Ablehnung länger dauert. Ich habe auch die Prüfung auf Division durch Null entfernt, um einen saubereren Vergleich mit dem sample_normals-Fall zu erhalten.

Durch das Entfernen der Vektorisierung aus den beiden direkten Methoden werden sie in den gleichen Bereich versetzt:

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

Ein anderer Weg, dass abhängig von der Hardware viel schneller sein könnte.

Wählen Sie a, b, c als drei Zufallszahlen zwischen -1 und 1

r2 = a^2 + b^2 + c^2 berechnen

Wenn r2> 1.0 (= der Punkt liegt nicht in der Kugel) oder r2 <0.00001 (= der Punkt liegt zu nahe am Zentrum, werden wir durch die Projektion auf die Oberfläche der Kugel durch Null geteilt), verwerfen Sie die Werte und wähle einen anderen Satz von zufälligen a, b, c

Ansonsten haben Sie Ihren zufälligen Punkt (relativ zum Mittelpunkt der Kugel):

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

(bearbeitet, um Korrekturen aus Kommentaren wiederzugeben)

ich habe 2004 einige zeitlich konstante Ansätze für dieses Problem untersucht.

angenommen, Sie arbeiten in sphärischen Koordinaten, wobei theta der Winkel um die vertikale Achse (z. B. der Längengrad) und phi der Winkel ist, der vom Äquator nach oben angehoben wird (z. B. der Breitengrad) Gleichmäßige Verteilung von zufälligen Punkten auf der Hemisphäre nördlich des Äquators

1. wähle theta = Rand (0, 360).
2. wähle phi = 90 * (1 - sqrt (Rand (0, 1))).

um Punkte auf einer Kugel anstatt auf einer Halbkugel zu erhalten, negieren Sie phi in 50% der Fälle.

für die Neugierigen gilt ein ähnlicher Ansatz für die Erzeugung gleichmäßig verteilter Punkte auf einer Einheitsscheibe:

1. wähle theta = Rand (0, 360).
2. wähle radius = sqrt (Rand (0, 1)).

ich habe keine Beweise für die Korrektheit dieser Ansätze, aber ich habe sie in den letzten zehn Jahren sehr erfolgreich angewendet und bin von ihrer Korrektheit überzeugt.

einige Beispiele (ab 2004) für die verschiedenen Ansätze sind hier , einschließlich einer Visualisierung des Ansatzes, Punkte auf der Oberfläche eines Würfels auszuwählen und sie auf der Kugel zu normalisieren.