Solución:
El algoritmo de la esfera de Fibonacci es excelente para esto. Es rápido y da resultados que de un vistazo engañarán fácilmente al ojo humano. Puede ver un ejemplo realizado con procesamiento que mostrará el resultado a lo largo del tiempo a medida que se agregan puntos. Aquí hay otro gran ejemplo interactivo creado por @gman. Y aquí hay una implementación simple en Python.
import math
def fibonacci_sphere(samples=1):
points = []
phi = math.pi * (3. - math.sqrt(5.)) # golden angle in radians
for i in range(samples):
y = 1 - (i / float(samples - 1)) * 2 # y goes from 1 to -1
radius = math.sqrt(1 - y * y) # radius at y
theta = phi * i # golden angle increment
x = math.cos(theta) * radius
z = math.sin(theta) * radius
points.append((x, y, z))
return points
1000 muestras te dan esto:
El método de la espiral dorada
Dijiste que no podías hacer funcionar el método de la espiral dorada y eso es una pena porque es realmente bueno. Me gustaría darte una comprensión completa de esto para que tal vez puedas entender cómo evitar que esto se “acumule”.
Así que aquí hay una forma rápida y no aleatoria de crear una celosía que sea aproximadamente correcta; como se discutió anteriormente, ninguna celosía será perfecta, pero esto puede ser lo suficientemente bueno. Se compara con otros métodos, por ejemplo, en BendWavy.org, pero tiene un aspecto agradable y bonito, así como una garantía sobre el espaciado uniforme en el límite.
Imprimación: espirales de girasol en el disco unitario
Para comprender este algoritmo, primero lo invito a mirar el algoritmo de espiral de girasol 2D. Esto se basa en el hecho de que el número más irracional es la proporción áurea. (1 + sqrt(5))/2
y si uno emite puntos por el enfoque “párese en el centro, gire una proporción áurea de giros completos, luego emita otro punto en esa dirección”, uno naturalmente construye una espiral que, a medida que llega a un número cada vez mayor de puntos, no obstante se niega a tener “barras” bien definidas en las que los puntos se alineen.(Nota 1.)
El algoritmo para el espaciado uniforme en un disco es,
from numpy import pi, cos, sin, sqrt, arange
import matplotlib.pyplot as pp
num_pts = 100
indices = arange(0, num_pts, dtype=float) + 0.5
r = sqrt(indices/num_pts)
theta = pi * (1 + 5**0.5) * indices
pp.scatter(r*cos(theta), r*sin(theta))
pp.show()
y produce resultados que se parecen a (n = 100 yn = 1000):
Espaciar los puntos radialmente
La clave extraña es la fórmula r = sqrt(indices / num_pts)
; ¿Cómo llegué a ese? (Nota 2.)
Bueno, estoy usando la raíz cuadrada aquí porque quiero que tengan un espacio uniforme alrededor del disco. Eso es lo mismo que decir que en el límite de grandes norte Quiero una pequeña región R ∈ (r, r + dr), Θ ∈ (θ, θ + dθ) para contener un número de puntos proporcional a su área, que es r Dr Dθ. Ahora bien, si pretendemos que estamos hablando de una variable aleatoria aquí, esto tiene una interpretación sencilla como decir que la densidad de probabilidad conjunta para (R, Θ) es solo cr por alguna constante C. La normalización en el disco de la unidad forzaría C = 1 / π.
Ahora permítanme presentarles un truco. Proviene de la teoría de la probabilidad, donde se conoce como muestreo de la CDF inversa: suponga que desea generar una variable aleatoria con una densidad de probabilidad F(z) y tienes una variable aleatoria U ~ Uniforme (0, 1), como sale de random()
en la mayoría de los lenguajes de programación. ¿Cómo haces esto?
- Primero, convierta su densidad en una función de distribución acumulativa o CDF, que llamaremos F(z). Un CDF, recuerde, aumenta monótonamente de 0 a 1 con la derivada F(z).
- Luego calcula la función inversa de la CDF F-1(z).
- Encontraras eso Z = F-1(U) se distribuye de acuerdo con la densidad objetivo. (Nota 3).
Ahora, el truco de la espiral de la proporción áurea espacia los puntos en un patrón uniforme para θ así que integremos eso; para el disco unitario nos queda F(r) = r2. Entonces la función inversa es F-1(tu) = tu1/2, y por lo tanto generaríamos puntos aleatorios en el disco en coordenadas polares con r = sqrt(random()); theta = 2 * pi * random()
.
Ahora en lugar de al azar muestreando esta función inversa estamos uniformemente muestreándolo, y lo bueno del muestreo uniforme es que nuestros resultados sobre cómo se distribuyen los puntos en el límite de grandes norte se comportará como si lo hubiéramos muestreado al azar. Esta combinación es el truco. En lugar de random()
usamos (arange(0, num_pts, dtype=float) + 0.5)/num_pts
, de modo que, digamos, si queremos muestrear 10 puntos, son r = 0.05, 0.15, 0.25, ... 0.95
. Muestreamos uniformemente r para obtener un espaciado de áreas iguales, y usamos el incremento de girasol para evitar terribles “barras” de puntos en la salida.
Ahora haciendo el girasol en una esfera
Los cambios que debemos hacer para salpicar la esfera con puntos simplemente implican cambiar las coordenadas polares por coordenadas esféricas. La coordenada radial, por supuesto, no entra en esto porque estamos en una esfera unitaria. Para mantener las cosas un poco más consistentes aquí, aunque fui entrenado como físico, usaré las coordenadas de los matemáticos donde 0 ≤ φ ≤ π es la latitud que desciende del polo y 0 ≤ θ ≤ 2π es la longitud. Entonces, la diferencia con respecto a lo anterior es que básicamente estamos reemplazando la variable r con φ.
Nuestro elemento de área, que fue r Dr Dθ, ahora se convierte en el pecado no mucho más complicado (φ) Dφ Dθ. Entonces, nuestra densidad conjunta para un espaciado uniforme es sin (φ) / 4π. Integrando θ, encontramos F(φ) = pecado (φ) / 2, por lo tanto F(φ) = (1 – cos (φ)) / 2. Al invertir esto, podemos ver que una variable aleatoria uniforme se vería como acos (1 – 2 tu), pero tomamos muestras de manera uniforme en lugar de aleatoria, por lo que en su lugar usamos φk = acos (1-2 (k + 0,5) /norte). Y el resto del algoritmo solo proyecta esto en las coordenadas x, y, z:
from numpy import pi, cos, sin, arccos, arange
import mpl_toolkits.mplot3d
import matplotlib.pyplot as pp
num_pts = 1000
indices = arange(0, num_pts, dtype=float) + 0.5
phi = arccos(1 - 2*indices/num_pts)
theta = pi * (1 + 5**0.5) * indices
x, y, z = cos(theta) * sin(phi), sin(theta) * sin(phi), cos(phi);
pp.figure().add_subplot(111, projection='3d').scatter(x, y, z);
pp.show()
Nuevamente, para n = 100 yn = 1000, los resultados se ven así:
Más investigación
Quería dar un saludo al blog de Martin Roberts. Tenga en cuenta que antes creé un desplazamiento de mis índices agregando 0.5 a cada índice. Esto fue visualmente atractivo para mí, pero resulta que la elección de la compensación es muy importante y no es constante durante el intervalo y puede significar obtener hasta un 8% más de precisión en el empaque si se elige correctamente. También debería haber una forma de obtener su R2 secuencia para cubrir una esfera y sería interesante ver si esto también produce una cubierta uniforme agradable, tal vez tal como está, pero tal vez necesita ser, digamos, tomado de solo la mitad del cuadrado unitario cortado en diagonal más o menos y estirado a conseguir un círculo.
Notas
-
Esas “barras” están formadas por aproximaciones racionales a un número, y las mejores aproximaciones racionales a un número provienen de su expresión fraccionaria continua,
z + 1/(n_1 + 1/(n_2 + 1/(n_3 + ...)))
dóndez
es un número entero yn_1, n_2, n_3, ...
es una secuencia finita o infinita de enteros positivos:def continued_fraction(r): while r != 0: n = floor(r) yield n r = 1/(r - n)
Dado que la parte de la fracción
1/(...)
está siempre entre cero y uno, un número entero grande en la fracción continua permite una aproximación racional particularmente buena: “uno dividido por algo entre 100 y 101” es mejor que “uno dividido por algo entre 1 y 2”. El número más irracional es, por tanto, el que es1 + 1/(1 + 1/(1 + ...))
y no tiene aproximaciones racionales particularmente buenas; uno puede resolver φ = 1 + 1 /φ multiplicando por φ para obtener la fórmula de la proporción áurea. -
Para las personas que no están tan familiarizadas con NumPy, todas las funciones están “vectorizadas”, de modo que
sqrt(array)
es lo mismo que otros idiomas podrían escribirmap(sqrt, array)
. Así que este es un componente por componente.sqrt
solicitud. Lo mismo también se aplica a la división por un escalar o la suma con escalares, que se aplican a todos los componentes en paralelo. -
La prueba es simple una vez que sabes que este es el resultado. Si preguntas cuál es la probabilidad de que z < Z < z + dz, esto es lo mismo que preguntar cuál es la probabilidad de que z < F-1(U) z + dz, solicitar F a las tres expresiones señalando que es una función que aumenta monótonamente, por lo tanto F(z) U < F(z + dz), expanda el lado derecho hacia afuera para encontrar F(z) + F(z) Dz, y desde U es uniforme esta probabilidad es solo F(z) Dz como fue prometido.
Esto se conoce como puntos de empaquetamiento en una esfera, y no existe una solución general perfecta (conocida). Sin embargo, hay muchas soluciones imperfectas. Los tres más populares parecen ser:
- Crea una simulación. Trate cada punto como un electrón restringido a una esfera, luego ejecute una simulación para un cierto número de pasos. La repulsión de los electrones tenderá naturalmente al sistema a un estado más estable, donde los puntos están lo más lejos posible unos de otros.
-
Rechazo de hipercubo. Este método que suena elegante es realmente simple: eliges puntos de manera uniforme (mucho mas que
n
de ellos) dentro del cubo que rodea la esfera, luego rechace los puntos fuera de la esfera. Trate los puntos restantes como vectores y normalícelos. Estas son sus “muestras” – elijan
de ellos usando algún método (aleatorio, codicioso, etc.). - Aproximaciones en espiral. Traza una espiral alrededor de una esfera y distribuye uniformemente los puntos alrededor de la espiral. Debido a las matemáticas involucradas, estos son más complicados de entender que la simulación, pero mucho más rápidos (y probablemente involucran menos código). El más popular parece ser el de Saff, et al.
A lote se puede encontrar más información sobre este problema aquí