Répartir uniformément n points sur une sphère

121

J'ai besoin d'un algorithme qui puisse me donner des positions autour d'une sphère pour N points (moins de 20, probablement) qui les étale vaguement. Il n'y a pas besoin de "perfection", mais j'en ai juste besoin pour qu'aucune d'entre elles ne soit regroupée.

  • Cette question a fourni un bon code, mais je n'ai pas trouvé de moyen de créer cet uniforme, car cela semblait aléatoire à 100%.
  • Ce billet de blog recommandé avait deux façons permettant la saisie du nombre de points sur la sphère, mais l' algorithme Saff et Kuijlaars est exactement dans psuedocode que je pourrais transcrire, et l' exemple de code que j'ai trouvé contenait "node [k]", que je ne pouvais pas voir expliqué et ruiné cette possibilité. Le deuxième exemple de blog était la spirale de la section dorée, qui m'a donné des résultats étranges et groupés, sans moyen clair de définir un rayon constant.
  • Cet algorithme de cette question semble pouvoir fonctionner, mais je ne peux pas reconstituer ce qui se trouve sur cette page en psuedocode ou quoi que ce soit.

Quelques autres fils de questions que j'ai rencontrés parlaient de distribution uniforme aléatoire, ce qui ajoute un niveau de complexité qui ne me préoccupe pas. Je m'excuse que ce soit une question tellement idiote, mais je voulais montrer que j'ai vraiment regardé attentivement et que j'ai toujours échoué.

Donc, ce que je recherche, c'est un pseudocode simple pour répartir uniformément N points autour d'une sphère unitaire, qui retourne en coordonnées sphériques ou cartésiennes. Encore mieux s'il peut même distribuer avec un peu de randomisation (pensez aux planètes autour d'une étoile, décemment réparties, mais avec une marge de manœuvre).

Arriver à
la source
Que voulez-vous dire "avec un peu de randomisation"? Voulez-vous dire des perturbations dans un certain sens?
ninjagecko
32
OP est confus. Ce qu'il cherche, c'est de mettre n-points sur une sphère, de sorte que la distance minimale entre deux points soit aussi grande que possible. Cela donnera aux points l'apparence d'être «uniformément répartis» sur toute la sphère. Cela n'a aucun rapport avec la création d'une distribution aléatoire uniforme sur une sphère, ce qui est le sujet de bon nombre de ces liens et ce dont parlent de nombreuses réponses ci-dessous.
BlueRaja - Danny Pflughoeft
1
20, ce n'est pas beaucoup de points à placer sur une sphère si vous ne voulez pas qu'ils aient l'air juste au hasard.
John Alexiou
2
Voici un moyen de le faire (il a des exemples de code): pdfs.semanticscholar.org/97a6/… (on dirait qu'il utilise des calculs de force de répulsion)
trusktr
1
Bien sûr, pour les valeurs sur N dans {4, 6, 8, 12, 20}, il existe des solutions exactes dans lesquelles la distance de chaque point à (chacun de) ses plus proches voisins est une constante pour tous les points et tous les voisins les plus proches.
dmckee --- ex-moderator chaton

Réponses:

13

Dans cet exemple de code, il node[k] n'y a que le kième nœud. Vous générez un tableau N points et node[k]est le kième (de 0 à N-1). Si c'est tout ce qui vous déroute, j'espère que vous pourrez l'utiliser maintenant.

(en d'autres termes, kest un tableau de taille N qui est défini avant le début du fragment de code et qui contient une liste des points).

Alternativement , en s'appuyant sur l'autre réponse ici (et en utilisant 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 vous tracez cela, vous verrez que l'espacement vertical est plus grand près des pôles de sorte que chaque point est situé à peu près dans la même surface totale de l'espace (près des pôles, il y a moins d'espace "horizontalement", donc cela donne plus "verticalement" ).

Ce n'est pas la même chose que tous les points ayant à peu près la même distance par rapport à leurs voisins (ce dont je pense que vos liens parlent), mais cela peut être suffisant pour ce que vous voulez et améliore simplement la création d'une grille de latitude / longitude uniforme .

Andrew Cooke
la source
sympa, c'est bon de voir une solution mathématique. Je pensais utiliser une séparation en hélice et en longueur d'arc. Je ne sais toujours pas comment obtenir la solution optimale, ce qui est un problème intéressant.
robert king
avez-vous vu que j'ai modifié ma réponse pour inclure une explication du nœud [k] en haut? Je pense que c'est peut-être tout ce dont vous avez besoin ...
andrew cooke
Merveilleux, merci pour l'explication. Je vais l'essayer plus tard, car je n'ai pas le temps actuellement, mais merci beaucoup de m'avoir aidé. Je vous ferai savoir comment cela finit par fonctionner pour mes besoins. ^^
Befall
L'utilisation de la méthode Spiral répond parfaitement à mes besoins, merci beaucoup pour l'aide et la clarification. :)
Befall
13
Le lien semble mort.
Scheintod
140

L'algorithme de la sphère de Fibonacci est idéal pour cela. Il est rapide et donne des résultats qui, en un coup d'œil, tromperont facilement l'œil humain. Vous pouvez voir un exemple réalisé avec le traitement qui affichera le résultat au fil du temps à mesure que des points sont ajoutés. Voici un autre excellent exemple interactif réalisé par @gman. Et voici une implémentation 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 échantillons vous donnent ceci:

entrez la description de l'image ici

Fnord
la source
une variable n est appelée lors de la définition de phi: phi = ((i + rnd)% n) * incrément. Est-ce que n = échantillons?
Andrew Staroscik
@AndrewStaroscik oui! Quand j'ai écrit le code pour la première fois, j'ai utilisé "n" comme variable et j'ai changé le nom plus tard, mais je n'ai pas fait de diligence raisonnable. Merci d'avoir attrapé ça!
Fnord
4
@Xarbrough le code vous donne des points autour d'une sphère unitaire, donc multipliez simplement chaque point par le scalaire que vous voulez pour le rayon.
Fnord
2
@Fnord: Pouvons-nous faire cela pour des dimensions supérieures?
pikachuchameleon
108

La méthode de la spirale d'or

Vous avez dit que vous ne pouviez pas faire fonctionner la méthode de la spirale dorée et c'est dommage parce que c'est vraiment très bon. Je voudrais vous donner une compréhension complète de cela afin que vous puissiez peut-être comprendre comment éviter que cela ne soit «groupé».

Voici donc un moyen rapide et non aléatoire de créer un treillis qui est approximativement correct; comme indiqué ci-dessus, aucun réseau ne sera parfait, mais cela peut être suffisant. Il est comparé à d'autres méthodes, par exemple sur BendWavy.org, mais il a juste un joli et joli look ainsi qu'une garantie d'espacement uniforme dans la limite.

Primer: spirales de tournesol sur le disque de l'unité

Pour comprendre cet algorithme, je vous invite d'abord à vous pencher sur l'algorithme 2D de la spirale de tournesol. Ceci est basé sur le fait que le nombre le plus irrationnel est le nombre d'or (1 + sqrt(5))/2et si l'on émet des points par l'approche «se tenir au centre, tourner un nombre d'or de tours entiers, puis émettre un autre point dans cette direction», on construit naturellement un spirale qui, à mesure que vous atteignez des nombres de points de plus en plus élevés, refuse néanmoins d'avoir des «barres» bien définies sur lesquelles les points s'alignent. (Note 1.)

L'algorithme pour l'espacement pair sur un disque est,

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()

et il produit des résultats qui ressemblent à (n = 100 et n = 1000):

entrez la description de l'image ici

Espacement radial des points

La principale chose étrange est la formule r = sqrt(indices / num_pts); comment suis-je arrivé à celui-là? (Note 2.)

Eh bien, j'utilise la racine carrée ici parce que je veux que celles-ci aient un espacement de zone uniforme autour du disque. Cela revient à dire que dans la limite du grand N je veux qu'une petite région R ∈ ( r , r + d r ), Θ ∈ ( θ , θ + d θ ) contienne un nombre de points proportionnel à son aire, qui est r d r d θ . Maintenant, si nous prétendons parler ici d'une variable aléatoire, cela a une interprétation simple comme disant que la densité de probabilité conjointe pour ( R , Θ ) est juste crpour une certaine constante c . La normalisation sur le disque unité forcerait alors c = 1 / π.

Maintenant, laissez-moi vous présenter une astuce. Cela vient de la théorie des probabilités où il est connu sous le nom d' échantillonnage du CDF inverse : supposons que vous vouliez générer une variable aléatoire avec une densité de probabilité f ( z ) et que vous ayez une variable aléatoire U ~ Uniform (0, 1), tout comme sort de random()dans la plupart des langages de programmation. Comment est-ce que tu fais ça?

  1. Tout d'abord, transformez votre densité en une fonction de distribution cumulative ou CDF, que nous appellerons F ( z ). Un CDF, rappelez-vous, augmente de façon monotone de 0 à 1 avec le dérivé f ( z ).
  2. Calculez ensuite la fonction inverse F -1 ( z ) du CDF .
  3. Vous constaterez que Z = F -1 ( U ) est distribué en fonction de la densité cible. (Note 3).

Maintenant, l'astuce en spirale du nombre d'or espace les points dans un motif bien uniforme pour θ alors intégrons cela; pour le disque unitaire, on se retrouve avec F ( r ) = r 2 . Donc, la fonction inverse est F -1 ( u ) = u 1/2 , et donc nous générerions des points aléatoires sur le disque en coordonnées polaires avec r = sqrt(random()); theta = 2 * pi * random().

Maintenant, au lieu d' échantillonner au hasard cette fonction inverse, nous l' échantillonnons uniformément , et la bonne chose à propos de l'échantillonnage uniforme est que nos résultats sur la façon dont les points sont répartis dans la limite d'un grand N se comporteront comme si nous l'avions échantillonné au hasard. Cette combinaison est le truc. Au lieu de random()nous utiliser (arange(0, num_pts, dtype=float) + 0.5)/num_pts, de sorte que, disons, si nous voulons échantillonner 10 points, ils le sont r = 0.05, 0.15, 0.25, ... 0.95. Nous échantillonnons uniformément r pour obtenir un espacement de surface égale, et nous utilisons l'incrément de tournesol pour éviter de terribles «barres» de points dans la sortie.

Maintenant en train de faire le tournesol sur une sphère

Les changements que nous devons apporter pour doter la sphère de points impliquent simplement de remplacer les coordonnées polaires par des coordonnées sphériques. La coordonnée radiale n'entre bien sûr pas dans cela car nous sommes sur une sphère unitaire. Pour garder les choses un peu plus cohérentes ici, même si j'ai été formé en tant que physicien, j'utiliserai les coordonnées des mathématiciens où 0 ≤ φ ≤ π est la latitude descendant du pôle et 0 ≤ θ ≤ 2π est la longitude. Donc, la différence par rapport à ci-dessus est que nous remplaçons fondamentalement la variable r par φ .

Notre élément d'aire, qui était r d r d θ , devient maintenant le sin ( φ ) d φ d θ pas beaucoup plus compliqué . Donc, notre densité conjointe pour un espacement uniforme est sin ( φ ) / 4π. En intégrant θ , nous trouvons f ( φ ) = sin ( φ ) / 2, donc F ( φ ) = (1 - cos ( φ )) / 2. En inversant cela, nous pouvons voir qu'une variable aléatoire uniforme ressemblerait à acos (1 - 2 u ), mais nous échantillonnons uniformément plutôt que de manière aléatoire, donc nous utilisons à la place φ k = acos (1 - 2 ( k+ 0,5) / N ). Et le reste de l'algorithme projette simplement ceci sur les coordonnées x, y et 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()

Encore une fois pour n = 100 et n = 1000, les résultats ressemblent à: entrez la description de l'image ici entrez la description de l'image ici

De plus amples recherches

Je voulais crier au blog de Martin Roberts. Notez que ci-dessus j'ai créé un offset de mes indices en ajoutant 0,5 à chaque index. C'était juste visuellement attrayant pour moi, mais il s'avère que le choix du décalage compte beaucoup et n'est pas constant sur l'intervalle et peut signifier obtenir jusqu'à 8% de meilleure précision dans l'emballage s'il est choisi correctement. Il devrait également y avoir un moyen de faire en sorte que sa séquence R 2 couvre une sphère et il serait intéressant de voir si cela produisait également une belle couverture uniforme, peut-être telle quelle mais peut-être devant être, par exemple, prise à partir de seulement la moitié de l'unité carrée coupée en diagonale environ et étirée pour former un cercle.

Remarques

  1. Ces «barres» sont formées par des approximations rationnelles d'un nombre, et les meilleures approximations rationnelles d'un nombre proviennent de son expression de fraction continue, z + 1/(n_1 + 1/(n_2 + 1/(n_3 + ...)))zest un entier et n_1, n_2, n_3, ...est une séquence finie ou infinie d'entiers positifs:

    def continued_fraction(r):
        while r != 0:
            n = floor(r)
            yield n
            r = 1/(r - n)

    Comme la partie de fraction 1/(...)est toujours comprise entre zéro et un, un grand entier dans la fraction continue permet une approximation rationnelle particulièrement bonne: "un divisé par quelque chose entre 100 et 101" est meilleur que "un divisé par quelque chose entre 1 et 2." Le nombre le plus irrationnel est donc celui qui n'est 1 + 1/(1 + 1/(1 + ...))et n'a pas d'approximations rationnelles particulièrement bonnes; on peut résoudre φ = 1 + 1 / φ en multipliant par φ pour obtenir la formule du nombre d'or.

  2. Pour les gens qui ne sont pas si familiers avec NumPy - toutes les fonctions sont «vectorisées», donc c'est sqrt(array)la même chose que ce que d'autres langages pourraient écrire map(sqrt, array). Il s'agit donc d'une application composant par composant sqrt. Il en va de même pour la division par un scalaire ou l'addition avec des scalaires - ceux-ci s'appliquent à tous les composants en parallèle.

  3. La preuve est simple une fois que vous savez que c'est le résultat. Si vous demandez quelle est la probabilité que z < Z < z + d z , cela revient à demander quelle est la probabilité que z < F -1 ( U ) < z + d z , appliquez F aux trois expressions en notant que c'est une fonction monotone croissante, d'où F ( z ) < U < F ( z + d z ), étend le côté droit pour trouver F ( z ) + f( z ) d z , et comme U est uniforme, cette probabilité est juste f ( z ) d z comme promis.

CR Drost
la source
4
Je ne sais pas pourquoi c'est si bas, c'est de loin la meilleure méthode rapide pour le faire.
whn
2
@snb merci pour les aimables paroles! il est si bas en partie parce qu'il est beaucoup, beaucoup plus jeune que toutes les autres réponses ici. Je suis surpris qu'il se porte aussi bien qu'il l'a été.
CR Drost
Une question qui me reste est: combien de points ai-je besoin de distribuer pour une distance maximale donnée entre deux points?
Felix D.
1
@FelixD. Cela ressemble à une question qui pourrait devenir très compliquée très rapidement, surtout si vous commencez à utiliser, par exemple, des distances de grand cercle plutôt que des distances euclidiennes. Mais peut-être que je peux répondre à une question simple, si l'on convertit les points de la sphère en leur diagramme de Voronoi, on peut décrire chaque cellule de Voronoi comme ayant une aire d'environ 4π / N et on peut convertir cela en une distance caractéristique en prétendant que c'est plutôt un cercle qu'un losange, πr² = 4π / N. Alors r = 2 / √ (N).
CR Drost
2
L'utilisation du théorème d'échantillonnage avec une entrée réellement uniforme au lieu d'une entrée aléatoire uniforme est l'une de ces choses qui me fait dire "Eh bien, pourquoi le # $% & n'y ai-je pas pensé?" . Agréable.
dmckee --- ex-moderator chaton
86

C'est ce qu'on appelle des points d'emballage sur une sphère, et il n'y a pas de solution générale parfaite (connue). Cependant, il existe de nombreuses solutions imparfaites. Les trois plus populaires semblent être:

  1. Créez une simulation . Traitez chaque point comme un électron contraint à une sphère, puis exécutez une simulation pendant un certain nombre d'étapes. La répulsion des électrons tendra naturellement le système vers un état plus stable, où les points sont à peu près aussi éloignés les uns des autres que possible.
  2. Rejet de l'hypercube . Cette méthode fantaisie est en fait très simple: vous choisissez uniformément des points (beaucoup plus nqu'eux) à l' intérieur du cube entourant la sphère, puis vous rejetez les points à l'extérieur de la sphère. Traitez les points restants comme des vecteurs et normalisez-les. Ce sont vos "échantillons" - choisissez n-les en utilisant une méthode (aléatoire, gourmande, etc.).
  3. Approximations en spirale . Vous tracez une spirale autour d'une sphère et répartissez uniformément les points autour de la spirale. En raison des mathématiques impliquées, celles-ci sont plus compliquées à comprendre que la simulation, mais beaucoup plus rapides (et impliquent probablement moins de code). Le plus populaire semble être celui de Saff et al .

Un beaucoup plus d' informations sur ce problème se trouve ici

BlueRaja - Danny Pflughoeft
la source
Je vais examiner la tactique en spirale qu'Andrew Cooke a publiée ci-dessous, cependant, pourriez-vous s'il vous plaît clarifier la différence entre ce que je veux et ce qu'est la "distribution aléatoire uniforme"? Est-ce juste un placement aléatoire à 100% des points sur une sphère afin qu'ils soient placés uniformément? Merci pour l'aide. :)
Befall
4
@Befall: "distribution aléatoire uniforme" se réfère à la distribution de probabilité étant uniforme - cela signifie que, lors du choix d'un point aléatoire sur la sphère, chaque point a une probabilité égale d'être choisi. Cela n'a rien à voir avec la répartition spatiale finale des points, et n'a donc rien à voir avec votre question.
BlueRaja - Danny Pflughoeft
Ahhh, d'accord, merci beaucoup. La recherche de ma question a conduit à une tonne de réponses pour les deux, et je ne pouvais pas vraiment comprendre ce qui était inutile pour moi.
Befall
Pour être clair, chaque point a une probabilité nulle d'être choisi. Le rapport des probabilités que le point appartienne à deux zones quelconques de la surface de la sphère est égal au rapport des surfaces.
AturSams
2
Le dernier lien est maintenant mort
Felix D.
10

Ce que vous recherchez s'appelle un revêtement sphérique . Le problème de la couverture sphérique est très difficile et les solutions sont inconnues sauf pour un petit nombre de points. Une chose que l'on sait avec certitude est que, étant donné n points sur une sphère, il existe toujours deux points de distance d = (4-csc^2(\pi n/6(n-2)))^(1/2)ou plus proches.

Si vous voulez une méthode probabiliste pour générer des points uniformément répartis sur une sphère, c'est simple: générer des points dans l'espace uniformément par distribution gaussienne (elle est intégrée à Java, pas difficile de trouver le code pour les autres langages). Donc dans un espace tridimensionnel, vous avez besoin de quelque chose comme

Random r = new Random();
double[] p = { r.nextGaussian(), r.nextGaussian(), r.nextGaussian() };

Projetez ensuite le point sur la sphère en normalisant sa distance par rapport à l'origine

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 distribution gaussienne en n dimensions est sphérique symétrique donc la projection sur la sphère est uniforme.

Bien sûr, il n'y a aucune garantie que la distance entre deux points quelconques dans une collection de points générés uniformément sera limitée ci-dessous, vous pouvez donc utiliser le rejet pour appliquer toutes les conditions que vous pourriez avoir: il est probablement préférable de générer toute la collection, puis rejeter toute la collection si nécessaire. (Ou utilisez le "rejet anticipé" pour rejeter toute la collection que vous avez générée jusqu'à présent; ne gardez simplement pas certains points et n'en perdez pas d'autres.) Vous pouvez utiliser la formule dci-dessus, moins une certaine marge, pour déterminer la distance minimale entre points en dessous desquels vous rejetterez un ensemble de points. Vous devrez calculer n choisir 2 distances, et la probabilité de rejet dépendra du mou; il est difficile de dire comment, alors lancez une simulation pour avoir une idée des statistiques pertinentes.

Edward Doolittle
la source
Vote positif pour les expressions de distance maximale minimale. Utile pour limiter le nombre de points que vous souhaitez utiliser. Une référence à une source faisant autorité pour cela serait bien, cependant.
dmckee --- ex-moderator chaton
6

Cette réponse est basée sur la même `` théorie '' qui est bien décrite par cette réponse

J'ajoute cette réponse comme suit:
- Aucune des autres options ne correspond au besoin «d'uniformité» «sur place» (ou pas de manière évidente). (En notant que la planète ressemble à un comportement de distribution particulièrement recherché dans la demande initiale, vous rejetez simplement de la liste finie des k points créés de manière uniforme au hasard (aléatoire par rapport au nombre d'index dans les k éléments).)
--Le plus proche un autre impl vous a forcé à décider du `` N '' par `` axe angulaire '', contre juste `` une valeur de N '' sur les deux valeurs de l'axe angulaire (ce qui à de faibles comptes de N est très difficile à savoir ce qui peut ou peut ne pas avoir d'importance ( Par exemple, vous voulez «5» points - amusez-vous))
- De plus, il est très difficile de `` comprendre '' comment différencier les autres options sans aucune image, alors voici à quoi ressemble cette option (ci-dessous), et l'implémentation prête à l'emploi qui va avec.

avec N à 20:

entrez la description de l'image ici
puis N à 80: entrez la description de l'image ici


voici le code python3 prêt à l'emploi, où l'émulation est la même source: " http://web.archive.org/web/20120421191837/http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere " trouvé par d'autres . (Le tracé que j'ai inclus, qui se déclenche lorsqu'il est exécuté en tant que `` main '', est tiré 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

testé à de faibles comptes (N en 2, 5, 7, 13, etc.) et semble fonctionner `` bien ''

Matt S.
la source
5

Essayer:

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 fonction ci-dessus doit fonctionner en boucle avec N boucle totale et k itération de courant de boucle.

Il est basé sur un motif de graines de tournesol, sauf que les graines de tournesol sont incurvées en un demi-dôme, puis à nouveau dans une sphère.

Voici une image, sauf que je mets la caméra à mi-chemin à l'intérieur de la sphère pour qu'elle ressemble à 2d au lieu de 3D car la caméra est à la même distance de tous les points. http://3.bp.blogspot.com/-9lbPHLccQHA/USXf88_bvVI/AAAAAAAAADY/j7qhQsSZsA8/s640/sphere.jpg

alientiel
la source
2

Healpix résout un problème étroitement lié (pixellisation de la sphère avec des pixels de surface égale):

http://healpix.sourceforge.net/

C'est probablement exagéré, mais peut-être qu'après l'avoir regardé, vous vous rendrez compte que certaines de ses autres belles propriétés vous intéressent. C'est bien plus qu'une simple fonction qui génère un nuage de points.

J'ai atterri ici en essayant de le retrouver; le nom "healpix" n'évoque pas exactement des sphères ...

Andrew Wagner
la source
1

avec un petit nombre de points, vous pouvez exécuter une simulation:

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)]
Robert King
la source
pour améliorer ma réponse, vous devriez changer le plus proche_index = i en le plus proche_index = randchoice (i, j)
robert king
1

Prenez les deux plus grands facteurs de votre N, si N==20les deux plus grands facteurs sont {5,4}, ou, plus généralement {a,b}. Calculer

dlat  = 180/(a+1)
dlong = 360/(b+1})

Mettez votre premier point à {90-dlat/2,(dlong/2)-180}, votre deuxième à {90-dlat/2,(3*dlong/2)-180}, votre troisième à {90-dlat/2,(5*dlong/2)-180}, jusqu'à ce que vous ayez fait le tour du monde une fois, à quel moment vous devez savoir {75,150}quand vous allez à côté {90-3*dlat/2,(dlong/2)-180}.

Évidemment, je travaille cela en degrés sur la surface de la terre sphérique, avec les conventions habituelles pour traduire +/- en N / S ou E / W. Et évidemment, cela vous donne une distribution complètement non aléatoire, mais elle est uniforme et les points ne sont pas regroupés.

Pour ajouter un certain degré d'aléatoire, vous pouvez générer 2 normalement distribués (avec une moyenne de 0 et un dev standard de {dlat / 3, dlong / 3} selon le cas) et les ajouter à vos points uniformément distribués.

Marque haute performance
la source
5
cela aurait l'air beaucoup mieux si vous travailliez en sin (lat) plutôt qu'en lat. en l'état, vous obtiendrez beaucoup de regroupement près des poteaux.
andrew cooke
1

edit: Cela ne répond pas à la question que l'OP avait l'intention de poser, la laissant ici au cas où les gens le trouveraient utile d'une manière ou d'une autre.

Nous utilisons la règle de multiplication de la probabilité, combinée à des nombres infinis. Il en résulte 2 lignes de code pour obtenir le résultat souhaité:

longitude: φ = uniform([0,2pi))
azimuth:   θ = -arcsin(1 - 2*uniform([0,1]))

(défini dans le système de coordonnées suivant :)

entrez la description de l'image ici

Votre langue a généralement une primitive de nombre aléatoire uniforme. Par exemple, en python, vous pouvez utiliser random.random()pour renvoyer un nombre dans la plage [0,1). Vous pouvez multiplier ce nombre par k pour obtenir un nombre aléatoire dans la plage [0,k). Ainsi en python, uniform([0,2pi))signifierait random.random()*2*math.pi.


Preuve

Maintenant, nous ne pouvons pas assigner θ uniformément, sinon nous aurions une agglutination aux pôles. Nous souhaitons attribuer des probabilités proportionnelles à la surface du coin sphérique (le θ dans ce diagramme est en fait φ):

entrez la description de l'image ici

Un déplacement angulaire dφ à l'équateur se traduira par un déplacement de dφ * r. Quel sera ce déplacement à un azimut arbitraire θ? Eh bien, le rayon de l'axe z est r*sin(θ), donc la longueur d'arc de cette "latitude" coupant le coin est dφ * r*sin(θ). Ainsi nous calculons la distribution cumulative de la zone à échantillonner, en intégrant la zone de la tranche du pôle sud au pôle nord.

entrez la description de l'image ici(où truc = dφ*r)

Nous allons maintenant essayer d'obtenir l'inverse du CDF pour en échantillonner: http://en.wikipedia.org/wiki/Inverse_transform_sampling

Tout d'abord, nous normalisons en divisant notre quasi-CDF par sa valeur maximale. Ceci a pour effet secondaire d'annuler les dφ et r.

azimuthalCDF: cumProb = (sin(θ)+1)/2 from -pi/2 to pi/2

inverseCDF: θ = -sin^(-1)(1 - 2*cumProb)

Donc:

let x by a random float in range [0,1]
θ = -arcsin(1-2*x)
ninjagecko
la source
n'est-ce pas équivalent à l'option qu'il a écartée comme étant «100% aléatoire»? je crois comprendre qu'il veut qu'ils soient espacés plus uniformément qu'une distribution aléatoire uniforme.
andrew cooke
@ BlueRaja-DannyPflughoeft: Hmm, c'est juste. Je suppose que je n'ai pas lu la question aussi attentivement que j'aurais dû. Je laisse cela ici de toute façon au cas où d'autres le trouveraient utile. Merci de l'avoir signalé.
ninjagecko
1

OU ... pour placer 20 points, calculez les centres des faces icosaèdres. Pour 12 points, trouvez les sommets de l'icosaèdre. Pour 30 points, le point médian des arêtes de l'icosaèdre. vous pouvez faire la même chose avec le tétraèdre, le cube, le dodécaèdre et les octaèdres: un ensemble de points se trouve sur les sommets, un autre au centre de la face et un autre au centre des arêtes. Ils ne peuvent cependant pas être mélangés.

user19371
la source
Une bonne idée, mais cela ne fonctionne que pour 4, 6, 8, 12, 20, 24 ou 30 points.
Le gars au chapeau
Si vous voulez tricher, vous pouvez utiliser le centre des faces et des vertices. Ils ne seront pas équi-espacés mais une approximation décente. C'est bien parce que c'est déterministe.
chessofnerd
0
# 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)])
ksmith
la source
4
Il serait utile que vous écriviez un texte expliquant ce que cela est censé faire, pour que le PO n'ait pas à croire que cela fonctionnera.
hcarver du
0

@robert king C'est une solution vraiment sympa mais qui contient quelques bugs bâclés. Je sais que cela m'a beaucoup aidé, alors ne faites pas attention à la négligence. :) Voici une version nettoyée ....

from math import pi, asin, sin, degrees
halfpi, twopi = .5 * pi, 2 * pi
sphere_area = lambda R=1.0: 4 * pi * R ** 2

lat_dist = lambda lat, R=1.0: R*(1-sin(lat))

#A = 2*pi*R^2(1-sin(lat))
def sphere_latarea(lat, R=1.0):
    if -halfpi > lat or lat > halfpi:
        raise ValueError("lat must be between -halfpi and halfpi")
    return 2 * pi * R ** 2 * (1-sin(lat))

sphere_lonarea = lambda lon, R=1.0: \
        4 * pi * R ** 2 * lon / twopi

#A = 2*pi*R^2 |sin(lat1)-sin(lat2)| |lon1-lon2|/360
#    = (pi/180)R^2 |sin(lat1)-sin(lat2)| |lon1-lon2|
sphere_rectarea = lambda lat0, lat1, lon0, lon1, R=1.0: \
        (sphere_latarea(lat0, R)-sphere_latarea(lat1, R)) * (lon1-lon0) / twopi


def test_sphere(n_lats=10, n_lons=19, radius=540.0):
    total_area = 0.0
    for i_lons in range(n_lons):
        lon0 = twopi * float(i_lons) / n_lons
        lon1 = twopi * float(i_lons+1) / n_lons
        for i_lats in range(n_lats):
            lat0 = asin(2 * float(i_lats) / n_lats - 1)
            lat1 = asin(2 * float(i_lats+1)/n_lats - 1)
            area = sphere_rectarea(lat0, lat1, lon0, lon1, radius)
            print("{:} {:}: {:9.4f} to  {:9.4f}, {:9.4f} to  {:9.4f} => area {:10.4f}"
                    .format(i_lats, i_lons
                    , degrees(lat0), degrees(lat1)
                    , degrees(lon0), degrees(lon1)
                    , area))
            total_area += area
    print("total_area = {:10.4f} (difference of {:10.4f})"
            .format(total_area, abs(total_area) - sphere_area(radius)))

test_sphere()
Ismael Harun
la source
-1

Cela fonctionne et c'est d'une simplicité mortelle. Autant de points que vous le souhaitez:

    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());
    }
Arthur Fleur
la source