均匀地在球体上分布n个点

Question

我需要一个算法,可以给我一个球体周围的位置N点(可能小于20),模糊地将它们展开.没有必要"完美",但我只是需要它,所以没有一个被捆绑在一起.

• 这个问题提供了很好的代码,但是我找不到制作这种制服的方法,因为这似乎是100%随机化的.
• 这个博客文章推荐有两种方式允许输入球体上的点数,但Saff和Kuijlaars算法完全是我可以转录的伪代码,我发现的代码示例包含"node [k]",我不能看到解释并破坏了这种可能性.第二个博客的例子是黄金分割螺旋,它给了我奇怪的,褶皱的结果,没有明确的方法来定义恒定的半径.
• 这种算法从这个问题好像它可能工作,但我无法拼凑出那是什么网页上成伪代码或任何东西.

我遇到的一些其他问题主题是随机均匀分布,这增加了我不关心的复杂程度.我很抱歉这是一个如此愚蠢的问题,但我想表明我真的很努力,但仍然很短暂.

所以,我正在寻找的是简单的伪代码,可以在单位球体周围均匀分布N个点,这些点可以返回球形或笛卡尔坐标.如果它甚至可以通过一点随机分布来更好(想想围绕恒星的行星,分散得很好,但有余地的余地).

Answer 1

Fibonacci球体算法非常适合这种情况.它速度快,结果一目了然很容易愚弄人眼.您可以看到完成处理的示例,该示例将随着时间的推移显示添加点的结果.这是 @gman制作的另一个很棒的互动示例.这是一个带有简单随机选项的快速python版本:

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个样本给你这个:

Answer 2

这被称为球体上的包装点,并且没有(已知的)通用的完美解决方案.但是,有很多不完美的解决方案.最受欢迎的三个似乎是:

1. 创建一个模拟.将每个点视为约束到球体的电子,然后运行模拟一定数量的步骤.电子的排斥力自然会使系统处于更稳定的状态,在这种状态下,这些点可以得到的距离彼此相差很远.
2. 超立方体拒绝.这种奇特的方法实际上非常简单:你在球体周围的立方体内统一选择点(远远超过n它们),然后拒绝球体外的点.将剩余的点视为向量,并将它们标准化.这些是你的"样本" - n使用某种方法(随机,贪婪等)选择它们.
3. 螺旋近似.您在球体周围追踪螺旋线,并均匀分布螺旋周围的点.由于所涉及的数学,这些比模拟更复杂,但更快(并且可能涉及更少的代码).最受欢迎的似乎是Saff等人.

一个很多关于这个问题的更多信息,可以发现这里

Answer 3

金色螺旋方法

你说你无法使用金色螺旋方法工作,这是一种耻辱,因为它真的非常非常好.我想让你对它有一个完整的了解,这样你就可以理解如何避免它被"束缚".

所以这是一种快速,非随机的方法来创建一个近似正确的晶格; 如上所述,没有格子是完美的,但这可能"足够好".它与其他方法进行了比较,例如BendWavy.org,但它只是有一个漂亮漂亮的外观以及有限的均匀间距的保证.

底漆:单位圆盘上的向日葵螺旋

为了理解这个算法,我首先邀请您来看看2D向日葵螺旋算法.这是基于这样一个事实,即最无理的数字是黄金比例(1 + sqrt(5))/2,如果通过"站在中心,转动整个转弯的黄金比率,然后在该方向上发出另一个点"的方法发出点,一个人自然地构建了一个螺旋,当你得到越来越多的点数时,仍然拒绝使用明确定义的"条形"来表示这些点.^(注1)

磁盘上均匀间隔的算法是,

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

它产生的结果看起来像(n = 100和n = 1000):

径向间隔点

关键奇怪的是公式r = sqrt(indices / num_pts); 我怎么来那个？^(笔记2.)

好吧,我在这里使用平方根,因为我希望这些在球体周围有均匀的区域间距.这等于说,在大的限度Ñ我想一点区域ř ∈([R ,- [R + d - [R ),Θ ∈(θ,θ + d θ)来包含正比于它的面积,一个点的数量,其是- [R d - [R d θ.现在,如果我们假装我们在这里谈论一个随机变量,这有一个直截了当的解释,即(R,Θ)的联合概率密度对于某些常数c来说只是cr.然后单位磁盘上的归一化强制c = 1 /π.

现在让我介绍一个技巧.它来自概率论,它被称为逆CDF的采样:假设你想要生成一个概率密度为f(z)的随机变量,你有一个随机变量U ~Uniform(0,1),就像出来的那样random()在大多数编程语言中.你怎么做到这一点？

1. 首先,将您的密度转换为累积分布函数 F(z),记住,它使用导数f(z)从0单调增加到1 .
2. 然后计算CDF的反函数F ^-1(z).
3. 你会发现Z = F ^-1(U)是根据目标密度分布的.^(注3).

现在,黄金比例螺旋技巧将这些点以θ的漂亮均匀模式对齐,所以让我们将它整合出来; 对于单位圆,我们留下F(r)= r ².因此,反函数是F ^-1(u)= u ^1/2,因此我们将在极坐标中生成球体上的随机点r = sqrt(random()); theta = 2 * pi * random().

现在我们不是对这个反函数进行随机抽样,而是对它进行均匀采样,而均匀采样的优点是我们关于点如何在大N的极限范围内展开的结果将表现得好像我们随机采样它一样.这种组合就是诀窍.而不是random()我们使用(arange(0, num_pts, dtype=float) + 0.5)/num_pts,所以,如果我们想要取样10分r = 0.05, 0.15, 0.25, ... 0.95.我们统一采样r得到等面积间距,我们使用向日葵增量来避免输出中可怕的"条"点.

现在在球体上做向日葵

我们需要用点来点球体只需要切换球坐标的极坐标.当然,径向坐标不会进入,因为我们在单位球上.为了让事情变得更加一致这里,尽管我是一位训练有素的物理学家,我会用数学家的坐标,其中0≤ φ ≤π是纬度来从极0下来≤ θ ≤2π是经度.因此,从上述不同的是,我们基本上取代了变量[R与φ.

我们的区域元素,这是[R d [R d θ,现在变成了没有,备受更复杂的罪孽(φ)d φ d θ.所以我们的均匀间距的连接密度是sin(φ)/4π.积分θ,我们发现f(φ)= sin(φ)/ 2,因此F(φ)=(1-cos(φ))/ 2.反相此我们可以看到,一个均匀随机变量看起来像ACOS(1 - 2 ü),但我们采样均匀,而不是随机的,所以我们改为使用φ _ķ = ACOS(1 - 2(ķ + 0.5)/ Ñ).算法的其余部分只是将其投影到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()

同样,对于n = 100和n = 1000,结果如下:

笔记

1. 那些"条形"是由一个数的有理逼近形成的,一个数的最佳有理逼近来自它的连续分数表达式,z + 1/(n_1 + 1/(n_2 + 1/(n_3 + ...)))其中z是一个整数,并且n_1, n_2, n_3, ...是正整数的有限或无限序列:

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

   由于分数部分1/(...)总是在0和1之间,因此连续分数中的大整数允许特别好的有理逼近:"1除以100和101之间的某个"优于"1除以1和2之间的某个".因此,最无理的数字是1 + 1/(1 + 1/(1 + ...))没有特别好的理性近似的数字; 一个可以解决φ = 1 + 1/φ通过乘以φ得到为黄金比例的公式.

   2. 对于那些不熟悉NumPy的人来说 - 所有的功能都是"矢量化"的,所以这sqrt(array)和其他语言的写法一样map(sqrt, array).所以这是一个逐个组件的sqrt应用程序.同样也适用于标量划分或添加标量 - 这些适用于并行的所有组件.

   3. 一旦你知道这是结果,证明很简单.如果你问z < Z < z + d z的概率是多少,这与询问z < F ^-1(U)< z + d z的概率是什么相同,将F应用于所有三个表达式,注意它是单调递增函数,因此F(z)< U < F(z + d z),将右侧展开以找到F(z)+ f(z)d z,并且因为U是均匀的,这个概率就是f(z)d z如承诺的那样.

Answer 4

在这个例子中,代码 node[k]只是第k个节点.您正在生成数组N个点并且node[k]是第k个(从0到N-1).如果这一切让您感到困惑,希望您现在可以使用它.

(换句话说,k是一个大小为N的数组,它在代码片段开始之前定义,并包含一个点列表).

或者,在此处构建另一个答案(并使用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

如果你绘制它,你会看到在极点附近的垂直间距更大,这样每个点位于大约相同的空间总面积(靠近极点的地方,"水平"的空间更小,因此它提供了更多的"垂直" ).

这与所有与其邻居具有相同距离的点(这是我认为您的链接所讨论的内容)不同,但它可能足以满足您的需求,并且可以简单地制作统一的lat/lon网格.

Answer 5

您正在寻找的是一种球形覆盖物.球形覆盖问题非常困难,除了少量的点之外,解决方案是未知的.有一点可以肯定的是,给定球体上的n个点,总是存在两个距离d = (4-csc^2(\pi n/6(n-2)))^(1/2)或更近的点.

如果你想要一种概率方法来生成均匀分布在球体上的点,那么很容易:通过高斯分布统一生成空间中的点(它内置于Java中,不难找到其他语言的代码).所以在三维空间中,你需要类似的东西

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

然后通过标准化距离原点的距离将点投影到球体上

double norm = Math.sqrt( (p[0])^2 + (p[1])^2 + (p[2])^2 ); 
double[] sphereRandomPoint = { p[0]/norm, p[1]/norm, p[2]/norm };

n维的高斯分布是球对称的,因此球体上的投影是均匀的.

当然,不能保证统一生成的点集合中任意两点之间的距离将限制在下面,因此您可以使用拒绝来强制执行您可能具有的任何此类条件:可能最好生成整个集合然后如有必要,拒绝整个收藏.(或者使用"早期拒绝"来拒绝你到目前为止生成的整个集合;只是不要保留一些点并放弃其他点.)你可以使用d上面给出的公式,减去一些松弛,来确定它们之间的最小距离您将拒绝一组积分的点数.你必须计算n选择2个距离,拒绝的概率取决于松弛; 很难说是怎么做的,所以运行一个模拟来了解相关的统计数据.

Answer 6

这个答案是基于这个答案所概述的相同"理论"

我将这个答案添加为:
- 没有其他选项适合"均匀性"需要"点亮"(或者显然不是那么明显).(注意在原始问题中让行星看起来像分布看起来特别需要的行为,你只是随机地从k统一创建的点的有限列表中拒绝(随机地回到k项中的索引计数).)
- 最接近其他impl强迫你通过'角轴'决定'N',而不是两个角轴值上的'N的一个值'(在N的低计数处,非常难以知道什么可能,或者可能无关紧要(例如,你想要'5'点 - 玩得开心))
- 此外,很难"不知道"如何在没有任何图像的情况下区分其他选项,所以这里是这个选项的样子(下图),并且准备就绪 -与它一起运行的实现.

N为20:

然后N在80:

这是可立即运行的python3代码,其中仿真是相同的来源:" http://web.archive.org/web/20120421191837/http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere " .(我所包含的绘图,以'main'形式运行时的触发,取自: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

在低计数(N,2,5,7,13等)测试,似乎工作'很好'

Answer 7

尝试：

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); 
};

上面的函数应该在循环中运行，总共有 N 个循环和 k 个循环当前迭代。

它基于向日葵种子图案，除了向日葵种子弯曲成半圆顶，然后再弯曲成球体。

这是一张图片，除了我将相机放在球体内部的一半处，因此它看起来是 2d 而不是 3d，因为相机与所有点的距离相同。 http://3.bp.blogspot.com/-9lbPHLccQHA/USXf88_bvVI/AAAAAAAAADY/j7qhQsSZsA8/s640/sphere.jpg

Answer 8

Healpix 解决了一个密切相关的问题（用相等面积的像素对球体进行像素化）：

http://healpix.sourceforge.net/

这可能有点矫枉过正，但也许在看过它之后，您会意识到它的其他一些不错的特性对您来说很有趣。它不仅仅是一个输出点云的函数。

我降落在这里试图再次找到它；“healpix”这个名字并不完全让人想起球体......