python - Distribuyendo uniformemente n puntos en una esfera
algorithm math (14)
El método de la espiral dorada
Dijiste que no podías hacer funcionar el método de la espiral dorada y es una pena porque es realmente, muy bueno. Me gustaría darte una comprensión completa de eso para que tal vez puedas entender cómo evitar que se "agrupe".
Así que aquí hay una manera rápida y no aleatoria de crear un enrejado que sea aproximadamente correcto; como se discutió anteriormente, ningún enrejado será perfecto, pero esto puede ser "lo suficientemente bueno". Se compara con otros métodos, por ejemplo, en BendWavy.org pero tiene un aspecto bonito y bonito, así como una garantía sobre el espaciado uniforme en el límite.
Primer: espirales de girasol en el disco de la unidad
Para comprender este algoritmo, primero lo invito a mirar el algoritmo en 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 "pararse en el centro, gire una proporción áurea de giros enteros, luego emita otro señale en esa dirección, "uno naturalmente construye una espiral que, a medida que obtiene un número cada vez mayor de puntos, se niega a tener" barras "bien definidas sobre las que alinearse los puntos. (Nota 1.)
El algoritmo para el espaciado par 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 (n = 100 yn = 1000):
Espaciando los puntos radialmente
Lo más extraño 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 estos tengan espaciado de área uniforme alrededor de la esfera. Eso es lo mismo que decir que en el límite de N grande quiero que una pequeña región R ∈ ( r , r + d r ), Θ ∈ ( θ , θ + d θ ) contenga una cantidad de puntos proporcional a su área, lo cual es r d r d θ . Ahora bien, si pretendemos que estamos hablando de una variable aleatoria aquí, esto tiene una interpretación directa al decir que la densidad de probabilidad conjunta para ( R , Θ ) es solo cr para alguna constante c . La normalización en el disco de la unidad fuerza a c = 1 / π.
Ahora déjame presentarte un truco. Viene de la teoría de la probabilidad donde se conoce como muestreo de la CDF inversa : supongamos que quiere generar una variable aleatoria con una densidad de probabilidad f ( z ) y tiene una variable aleatoria U ~ Uniform (0, 1), tal como sale de random() en la mayoría de los lenguajes de programación. ¿Cómo haces esto?
- Primero, convierte tu densidad en una función de distribución acumulativa F ( z ), que, recuerda, aumenta monótonamente de 0 a 1 con la derivada f ( z ).
- Luego calcule la función inversa CDF F -1 ( z ).
- Encontrará que Z = F -1 ( U ) se distribuye de acuerdo con la densidad objetivo. (Nota 3).
Ahora el truco en espiral de la relación de oro separa los puntos en un patrón agradablemente uniforme para θ así que integremos eso; para el círculo unitario nos quedamos con F ( r ) = r 2 . Entonces la función inversa es F -1 ( u ) = u 1/2 , y por lo tanto generaríamos puntos aleatorios en la esfera en coordenadas polares con r = sqrt(random()); theta = 2 * pi * random() r = sqrt(random()); theta = 2 * pi * random() .
Ahora, en lugar de muestrear aleatoriamente esta función inversa, la estamos muestreando uniformemente , y lo bueno del muestreo uniforme es que nuestros resultados sobre cómo se distribuyen los puntos en el límite del N grande se comportarán 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 , por lo que, digamos, si queremos muestrear 10 puntos, son r = 0.05, 0.15, 0.25, ... 0.95 . De manera uniforme, se muestra r para obtener un espaciado de igual área, y utilizamos el incremento de girasol para evitar "barras" terribles de puntos en la salida.
Ahora haciendo el girasol en una esfera
Los cambios que necesitamos hacer para acotar la esfera con puntos simplemente implican cambiar las coordenadas polares para las coordenadas esféricas. La coordenada radial, por supuesto, no entra en esto porque estamos en una esfera de unidad. 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 bajando del polo y 0 ≤ θ ≤ 2π es la longitud. Entonces la diferencia de arriba es que básicamente estamos reemplazando la variable r con φ .
Nuestro elemento de área, que era r d r d θ , ahora se convierte en el pecado no mucho más complicado ( φ ) d φ d θ . Por lo tanto, nuestra densidad conjunta para espaciado uniforme es sin ( φ ) / 4π. Integrando θ , encontramos f ( φ ) = sin ( φ ) / 2, por lo tanto F ( φ ) = (1 - cos ( φ )) / 2. Invirtiendo esto, podemos ver que una variable aleatoria uniforme se vería como acos (1 - 2 u ), pero muestreamos uniformemente en lugar de aleatoriamente, por lo que usamos φ k = acos (1 - 2 ( k + 0.5) / N ). Y el resto del algoritmo solo está proyectando esto en las coordenadas x, y, yz:
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í:
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 de fracción continua,
z + 1/(n_1 + 1/(n_2 + 1/(n_3 + ...)))dondezes 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)Como la fracción parte
1/(...)está siempre entre cero y uno, un 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 lo 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", por lo que
sqrt(array)es lo mismo que otros lenguajes podrían escribirmap(sqrt, array). Entonces esta es una aplicaciónsqrtcomponente por componente. Lo mismo también se aplica a la división mediante un escalar o una adición 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 pregunta cuál es la probabilidad de que z < Z < z + d z , esto es lo mismo que preguntar cuál es la probabilidad de que z < F -1 ( U ) < z + d z , aplique F a las tres expresiones y observe que es una función monótonamente creciente, por lo tanto, F ( z ) < U < F ( z + d z ), expande el lado derecho para encontrar F ( z ) + f ( z ) d z , y como U es uniforme, esta probabilidad es f ( z ) d z como se prometió.
Necesito un algoritmo que me pueda dar posiciones alrededor de una esfera para N puntos (menos de 20, probablemente) que los difunde vagamente. No hay necesidad de "perfección", pero solo lo necesito para que ninguno esté agrupado.
- Esta pregunta proporcionó un buen código, pero no pude encontrar una manera de hacer este uniforme, ya que parecía ser 100% aleatorio.
- Esta publicación de blog recomendada tenía dos formas para ingresar el número de puntos en la esfera, pero el algoritmo de Saff y Kuijlaars está exactamente en psuedocode que pude transcribir, y el ejemplo de código que encontré contenía "nodo [k]", que no pude ver explicado y arruinó esa posibilidad. El segundo ejemplo de blog fue el Golden Section Spiral, que me dio resultados extraños y agrupados, sin una forma clara de definir un radio constante.
- Este algoritmo de esta pregunta parece que podría funcionar, pero no puedo reconstruir lo que hay en esa página en psuedocode ni nada.
Algunos otros hilos de preguntas que encontré hablaban de la distribución aleatoria uniforme, que agrega un nivel de complejidad que no me preocupa. Me disculpo por el hecho de que esta es una pregunta tan tonta, pero quería mostrar que realmente me he esforzado y aún me he quedado corto.
Entonces, lo que estoy buscando es un pseudocódigo simple para distribuir uniformemente N puntos alrededor de una esfera de unidad, que retorna en coordenadas esféricas o cartesianas. Incluso mejor si puede incluso distribuirse con un poco de aleatorización (piense en planetas alrededor de una estrella, decentemente extendidos, pero con espacio para el margen de maniobra).
El algoritmo de la esfera Fibonacci es ideal para esto. Es rápido y da resultados que a primera vista engañarán fácilmente al ojo humano. Puede ver un ejemplo hecho con procesamiento que mostrará el resultado a lo largo del tiempo a medida que se agreguen puntos. Aquí hay otro gran ejemplo interactivo hecho por @gman. Y aquí hay una versión rápida de Python con una opción de aleatorización simple:
import math, random
def fibonacci_sphere(samples=1,randomize=True):
rnd = 1.
if randomize:
rnd = random.random() * samples
points = []
offset = 2./samples
increment = math.pi * (3. - math.sqrt(5.));
for i in range(samples):
y = ((i * offset) - 1) + (offset / 2);
r = math.sqrt(1 - pow(y,2))
phi = ((i + rnd) % samples) * increment
x = math.cos(phi) * r
z = math.sin(phi) * r
points.append([x,y,z])
return points
1000 muestras te dan esto:
Esta respuesta se basa en la misma "teoría" que se resume bien en esta respuesta
Estoy agregando esta respuesta como:
- Ninguna de las otras opciones se ajusta a la ''uniformidad'' necesaria ''de manera directa'' (o no obviamente, claramente). (Observando que el planeta tiene un comportamiento de aspecto de distribución particularmente deseado en la pregunta original, simplemente rechaza de la lista finita de los k puntos creados de manera uniforme al azar (al azar, con el recuento de índice en los elementos k hacia atrás)).
--La otra impl más cercana te obligó a decidir la ''N'' por ''eje angular'', frente a solo ''un valor de N'' en ambos valores de eje angular (que en recuentos bajos de N es muy complicado saber qué puede ser, o puede no importar (por ejemplo, quiere ''5'' puntos - diviértase))
- Además, es muy difícil ''asimilar'' cómo diferenciar entre las otras opciones sin ninguna imagen, así que esta es la apariencia de esta opción (abajo), y la implementación lista para ejecutar que la acompaña.
con N a los 20:
y luego N a 80:
aquí está el código python3 listo para ejecutar, donde la emulación es la misma fuente: " here " encontrado por otros . (El trazado que he incluido, que se dispara cuando se ejecuta como ''principal'', está tomado de: http://www.scipy.org/Cookbook/Matplotlib/mplot3D )
from math import cos, sin, pi, sqrt
def GetPointsEquiAngularlyDistancedOnSphere(numberOfPoints=45):
""" each point you get will be of form ''x, y, z''; in cartesian coordinates
eg. the ''l2 distance'' from the origion [0., 0., 0.] for each point will be 1.0
------------
converted from: http://web.archive.org/web/20120421191837/http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere )
"""
dlong = pi*(3.0-sqrt(5.0)) # ~2.39996323
dz = 2.0/numberOfPoints
long = 0.0
z = 1.0 - dz/2.0
ptsOnSphere =[]
for k in range( 0, numberOfPoints):
r = sqrt(1.0-z*z)
ptNew = (cos(long)*r, sin(long)*r, z)
ptsOnSphere.append( ptNew )
z = z - dz
long = long + dlong
return ptsOnSphere
if __name__ == ''__main__'':
ptsOnSphere = GetPointsEquiAngularlyDistancedOnSphere( 80)
#toggle True/False to print them
if( True ):
for pt in ptsOnSphere: print( pt)
#toggle True/False to plot them
if(True):
from numpy import *
import pylab as p
import mpl_toolkits.mplot3d.axes3d as p3
fig=p.figure()
ax = p3.Axes3D(fig)
x_s=[];y_s=[]; z_s=[]
for pt in ptsOnSphere:
x_s.append( pt[0]); y_s.append( pt[1]); z_s.append( pt[2])
ax.scatter3D( array( x_s), array( y_s), array( z_s) )
ax.set_xlabel(''X''); ax.set_ylabel(''Y''); ax.set_zlabel(''Z'')
p.show()
#end
probado en conteos bajos (N en 2, 5, 7, 13, etc.) y parece funcionar "bien"
Esto funciona y es mortalmente simple. Tantos puntos como quieras:
private function moveTweets():void {
var newScale:Number=Scale(meshes.length,50,500,6,2);
trace("new scale:"+newScale);
var l:Number=this.meshes.length;
var tweetMeshInstance:TweetMesh;
var destx:Number;
var desty:Number;
var destz:Number;
for (var i:Number=0;i<this.meshes.length;i++){
tweetMeshInstance=meshes[i];
var phi:Number = Math.acos( -1 + ( 2 * i ) / l );
var theta:Number = Math.sqrt( l * Math.PI ) * phi;
tweetMeshInstance.origX = (sphereRadius+5) * Math.cos( theta ) * Math.sin( phi );
tweetMeshInstance.origY= (sphereRadius+5) * Math.sin( theta ) * Math.sin( phi );
tweetMeshInstance.origZ = (sphereRadius+5) * Math.cos( phi );
destx=sphereRadius * Math.cos( theta ) * Math.sin( phi );
desty=sphereRadius * Math.sin( theta ) * Math.sin( phi );
destz=sphereRadius * Math.cos( phi );
tweetMeshInstance.lookAt(new Vector3D());
TweenMax.to(tweetMeshInstance, 1, {scaleX:newScale,scaleY:newScale,x:destx,y:desty,z:destz,onUpdate:onLookAtTween, onUpdateParams:[tweetMeshInstance]});
}
}
private function onLookAtTween(theMesh:TweetMesh):void {
theMesh.lookAt(new Vector3D());
}
Esto se conoce como puntos de empaque en una esfera, y no hay (conocida) solución general perfecta. 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 naturalmente tenderá al sistema a un estado más estable, donde los puntos están lo más alejados posible entre ellos.
- Rechazo de Hypercube . Este método de sonido sofisticado es realmente muy simple: usted elige puntos de manera uniforme (mucho más que
nde ellos) dentro del cubo que rodea la esfera, y luego rechaza los puntos fuera de la esfera. Trate los puntos restantes como vectores y normalícelos. Estas son sus "muestras": elijande ellas usando algún método (aleatorio, codicioso, etc.). - Aproximaciones en espiral Trazas una espiral alrededor de una esfera y distribuyes uniformemente los puntos alrededor de la espiral. Debido a las matemáticas involucradas, estas son más complicadas de entender que la simulación, pero mucho más rápidas (y probablemente involucren menos código). El más popular parece ser Saff, et al .
Se puede encontrar mucha más información sobre este problema here
Healpix resuelve un problema estrechamente relacionado (pixelación de la esfera con píxeles de igual área):
http://healpix.sourceforge.net/
Probablemente sea excesivo, pero tal vez después de verlo, te darás cuenta de que algunas de sus otras propiedades son interesantes para ti. Es mucho más que una función que genera una nube de puntos.
Aterricé aquí tratando de encontrarlo de nuevo; el nombre "healpix" no evoca exactamente esferas ...
Lo que estás buscando se llama cobertura esférica . El problema de cobertura esférica es muy difícil y las soluciones son desconocidas a excepción de un pequeño número de puntos. Una cosa que se sabe con certeza es que dados n puntos en una esfera, siempre existen dos puntos de distancia d = (4-csc^2(/pi n/6(n-2)))^(1/2) o más cerca.
Si desea un método probabilístico para generar puntos distribuidos uniformemente en una esfera, es fácil: generar puntos en el espacio uniformemente por distribución gaussiana (está integrado en Java, no es difícil encontrar el código para otros idiomas). Entonces en el espacio tridimensional, necesitas algo como
Random r = new Random();
double[] p = { r.nextGaussian(), r.nextGaussian(), r.nextGaussian() };
Luego proyecta el punto en la esfera normalizando su distancia desde el origen
double norm = Math.sqrt( (p[0])^2 + (p[1])^2 + (p[2])^2 );
double[] sphereRandomPoint = { p[0]/norm, p[1]/norm, p[2]/norm };
La distribución gaussiana en n dimensiones es esféricamente simétrica por lo que la proyección sobre la esfera es uniforme.
Por supuesto, no hay garantía de que la distancia entre dos puntos en una colección de puntos generados uniformemente se limite a continuación, por lo que puede usar el rechazo para imponer las condiciones que pueda tener: probablemente sea mejor generar toda la colección y luego rechazar toda la colección si es necesario. (O utilice el "rechazo temprano" para rechazar toda la colección que ha generado hasta ahora, simplemente no conserve algunos puntos y descarte otros). Puede usar la fórmula para d dada anteriormente, menos cierta holgura, para determinar la distancia mínima entre los puntos debajo de los cuales rechazarás un conjunto de puntos. Tendrá que calcular n elegir 2 distancias, y la probabilidad de rechazo dependerá de la holgura; es difícil decir cómo, así que ejecuta una simulación para tener una idea de las estadísticas relevantes.
O ... para colocar 20 puntos, calcule los centros de las caras icosaedroicas. Para 12 puntos, encuentra los vértices del icosaedro. Por 30 puntos, el punto medio de los bordes del icosaedro. puede hacer lo mismo con el tetraedro, el cubo, el dodecaedro y los octaedros: un conjunto de puntos está en los vértices, otro en el centro de la cara y otro en el centro de los bordes. Sin embargo, no se pueden mezclar.
Tome los dos factores más grandes de su N , si N==20 entonces los dos factores más grandes son {5,4} , o, más generalmente, {a,b} . Calcular
dlat = 180/(a+1)
dlong = 360/(b+1})
Pon tu primer punto en {90-dlat/2,(dlong/2)-180} , tu segundo en {90-dlat/2,(3*dlong/2)-180} , tu 3 ° en {90-dlat/2,(5*dlong/2)-180} , hasta que hayas tropezado alrededor del mundo una vez, para cuando tienes aproximadamente {75,150} cuando vayas a {90-3*dlat/2,(dlong/2)-180} .
Obviamente estoy trabajando esto en grados en la superficie de la tierra esférica, con las convenciones usuales para traducir +/- a N / S o E / W. Y obviamente esto te da una distribución completamente no aleatoria, pero es uniforme y los puntos no están agrupados.
Para agregar un cierto grado de aleatoriedad, podría generar 2 normalmente distribuidos (con media 0 y dev. Estándar de {dlat / 3, dlong / 3} según corresponda) y agréguelos a sus puntos distribuidos uniformemente.
Tratar:
function sphere ( N:float,k:int):Vector3 {
var inc = Mathf.PI * (3 - Mathf.Sqrt(5));
var off = 2 / N;
var y = k * off - 1 + (off / 2);
var r = Mathf.Sqrt(1 - y*y);
var phi = k * inc;
return Vector3((Mathf.Cos(phi)*r), y, Mathf.Sin(phi)*r);
};
La función anterior debe ejecutarse en bucle con N loop total e k loop current iteration.
Se basa en un patrón de semillas de girasol, excepto que las semillas de girasol están curvadas en una media cúpula y de nuevo en una esfera.
Aquí hay una imagen, excepto que puse la cámara a medio camino dentro de la esfera para que parezca 2d en lugar de 3d porque la cámara está a la misma distancia de todos los puntos. http://3.bp.blogspot.com/-9lbPHLccQHA/USXf88_bvVI/AAAAAAAAADY/j7qhQsSZsA8/s640/sphere.jpg
con pequeños números de puntos, podrías ejecutar una simulación:
from random import random,randint
r = 10
n = 20
best_closest_d = 0
best_points = []
points = [(r,0,0) for i in range(n)]
for simulation in range(10000):
x = random()*r
y = random()*r
z = r-(x**2+y**2)**0.5
if randint(0,1):
x = -x
if randint(0,1):
y = -y
if randint(0,1):
z = -z
closest_dist = (2*r)**2
closest_index = None
for i in range(n):
for j in range(n):
if i==j:
continue
p1,p2 = points[i],points[j]
x1,y1,z1 = p1
x2,y2,z2 = p2
d = (x1-x2)**2+(y1-y2)**2+(z1-z2)**2
if d < closest_dist:
closest_dist = d
closest_index = i
if simulation % 100 == 0:
print simulation,closest_dist
if closest_dist > best_closest_d:
best_closest_d = closest_dist
best_points = points[:]
points[closest_index]=(x,y,z)
print best_points
>>> best_points
[(9.921692138442777, -9.930808529773849, 4.037839326088124),
(5.141893371460546, 1.7274947332807744, -4.575674650522637),
(-4.917695758662436, -1.090127967097737, -4.9629263893193745),
(3.6164803265540666, 7.004158551438312, -2.1172868271109184),
(-9.550655088997003, -9.580386054762917, 3.5277052594769422),
(-0.062238110294250415, 6.803105171979587, 3.1966101417463655),
(-9.600996012203195, 9.488067284474834, -3.498242301168819),
(-8.601522086624803, 4.519484132245867, -0.2834204048792728),
(-1.1198210500791472, -2.2916581379035694, 7.44937337008726),
(7.981831370440529, 8.539378431788634, 1.6889099589074377),
(0.513546008372332, -2.974333486904779, -6.981657873262494),
(-4.13615438946178, -6.707488383678717, 2.1197605651446807),
(2.2859494919024326, -8.14336582650039, 1.5418694699275672),
(-7.241410895247996, 9.907335206038226, 2.271647103735541),
(-9.433349952523232, -7.999106443463781, -2.3682575660694347),
(3.704772125650199, 1.0526567864085812, 6.148581714099761),
(-3.5710511242327048, 5.512552040316693, -3.4318468250897647),
(-7.483466337225052, -1.506434920354559, 2.36641535124918),
(7.73363824231576, -8.460241422163824, -1.4623228616326003),
(10, 0, 0)]
editar: Esto no responde la pregunta que el OP quería hacer, dejándolo aquí en caso de que las personas lo encuentren útil de alguna manera.
Usamos la regla de la multiplicación de la probabilidad, combinada con infinitessimals. Esto da como resultado 2 líneas de código para lograr el resultado deseado:
longitude: φ = uniform([0,2pi))
azimuth: θ = -arcsin(1 - 2*uniform([0,1]))
(definido en el siguiente sistema de coordenadas :)
Generalmente, tu idioma tiene un primitivo de número aleatorio uniforme. Por ejemplo, en python puedes usar random.random() para devolver un número en el rango [0,1) . Puedes multiplicar este número por k para obtener un número aleatorio en el rango [0,k) . Por lo tanto, en python, uniform([0,2pi)) significaría random.random()*2*math.pi
Prueba
Ahora no podemos asignar θ uniformemente, de lo contrario nos agruparíamos en los polos. Deseamos asignar probabilidades proporcionales al área de superficie de la cuña esférica (el θ en este diagrama es en realidad φ):
Un desplazamiento angular dφ en el ecuador dará como resultado un desplazamiento de dφ * r. ¿Cuál será ese desplazamiento en un azimut arbitrario θ? Bueno, el radio del eje z es r*sin(θ) , por lo que la longitud de arco de esa "latitud" que se cruza con la cuña es dφ * r*sin(θ) . Por lo tanto, calculamos la distribución acumulada del área a muestrear a partir de ella, integrando el área del corte desde el polo sur al polo norte.
(donde cosas = dφ*r )
Ahora intentaremos obtener el inverso de la CDF para obtener muestras de ella: http://en.wikipedia.org/wiki/Inverse_transform_sampling
Primero, normalizamos dividiendo nuestro casi CDF por su valor máximo. Esto tiene el efecto secundario de cancelar el dφ y r.
azimuthalCDF: cumProb = (sin(θ)+1)/2 from -pi/2 to pi/2
inverseCDF: θ = -sin^(-1)(1 - 2*cumProb)
Así:
let x by a random float in range [0,1]
θ = -arcsin(1-2*x)
En este ejemplo, el node[k] código node[k] es solo el k-nodo. Está generando una matriz N puntos y el node[k] es la k (de 0 a N-1). Si eso es todo lo que te está confundiendo, ojalá puedas usar eso ahora.
(en otras palabras, k es una matriz de tamaño N que se define antes de que comience el fragmento de código, y que contiene una lista de los puntos).
Alternativamente , construyendo sobre la otra respuesta aquí (y usando Python):
> cat ll.py
from math import asin
nx = 4; ny = 5
for x in range(nx):
lon = 360 * ((x+0.5) / nx)
for y in range(ny):
midpt = (y+0.5) / ny
lat = 180 * asin(2*((y+0.5)/ny-0.5))
print lon,lat
> python2.7 ll.py
45.0 -166.91313924
45.0 -74.0730322921
45.0 0.0
45.0 74.0730322921
45.0 166.91313924
135.0 -166.91313924
135.0 -74.0730322921
135.0 0.0
135.0 74.0730322921
135.0 166.91313924
225.0 -166.91313924
225.0 -74.0730322921
225.0 0.0
225.0 74.0730322921
225.0 166.91313924
315.0 -166.91313924
315.0 -74.0730322921
315.0 0.0
315.0 74.0730322921
315.0 166.91313924
Si trazas eso, verás que el espaciado vertical es más grande cerca de los polos, de modo que cada punto se encuentra aproximadamente en la misma área total de espacio (cerca de los polos hay menos espacio "horizontalmente", por lo que da más "verticalmente") )
Esto no es lo mismo que todos los puntos que tienen aproximadamente la misma distancia con sus vecinos (que es de lo que creo que hablan sus enlaces), pero puede ser suficiente para lo que desea y mejora simplemente haciendo una grilla uniforme lat / lon .
# create uniform spiral grid
numOfPoints = varargin[0]
vxyz = zeros((numOfPoints,3),dtype=float)
sq0 = 0.00033333333**2
sq2 = 0.9999998**2
sumsq = 2*sq0 + sq2
vxyz[numOfPoints -1] = array([(sqrt(sq0/sumsq)),
(sqrt(sq0/sumsq)),
(-sqrt(sq2/sumsq))])
vxyz[0] = -vxyz[numOfPoints -1]
phi2 = sqrt(5)*0.5 + 2.5
rootCnt = sqrt(numOfPoints)
prevLongitude = 0
for index in arange(1, (numOfPoints -1), 1, dtype=float):
zInc = (2*index)/(numOfPoints) -1
radius = sqrt(1-zInc**2)
longitude = phi2/(rootCnt*radius)
longitude = longitude + prevLongitude
while (longitude > 2*pi):
longitude = longitude - 2*pi
prevLongitude = longitude
if (longitude > pi):
longitude = longitude - 2*pi
latitude = arccos(zInc) - pi/2
vxyz[index] = array([ (cos(latitude) * cos(longitude)) ,
(cos(latitude) * sin(longitude)),
sin(latitude)])