在球体上均匀分布n个点

我需要一个algorithm,可以让我围绕一个球体的位置N点(可能less于20),模糊地扩散出去。 没有必要“完美”,但我只是需要它,所以没有一个捆绑在一起。

  • 这个问题提供了很好的代码,但是我找不到一个统一的方法,因为这看起来是100%随机的。
  • 这个博客文章推荐有两种方法允许input球体上的点数,但是Saff和Kuijlaarsalgorithm正好在我能够转录的伪代码中 ,并且我find的代码示例包含了“node [k]”,我不能看到解释和破坏的可能性。 第二个博客的例子是黄金分割螺旋,它给了我奇怪的结果,没有明确的方法来定义一个恒定的半径。
  • 从这个问题的 algorithm似乎可能工作,但我不能拼凑在那个页面上的东西到psuedocode或任何东西。

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

所以,我正在寻找的是简单的伪代码来均匀分布单位球体周围的N个点,要么返回球面坐标或笛卡尔坐标。 甚至更好,如果它甚至可以散布一点随机(认为行星围绕一个明星,体面散开,但有余地余地)。

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

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

或者 ,在这里build立另一个答案(和使用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 

如果你绘制的话,你会发现在两极附近的垂直间距是较大的,所以每个点位于大约相同的总空间区域内(在极点附近空间“水平”的空间较小,所以它更“垂直” )。

这与所有与邻居有相同距离的点(这是我认为你的链接正在讨论的点)是不一样的,但是对于你想要和改进的方面来说,简单地做一个统一的纬度/经度网格。

这被称为一个球体上的包装点,并没有(已知的)一般的,完美的解决scheme。 但是,有很多不完善的解决scheme。 三个最受欢迎似乎是:

  1. 创build一个模拟 。 将每个点视为受限于一个球体的电子,然后对某个步骤进行模拟。 电子的排斥力自然会使系统趋于一个更稳定的状态,在这种状态下,各点之间的距离可以相差很远。
  2. 超立方体拒绝 。 这种奇特的方法实际上很简单:在球体周围的立方体内部统一select点(远多于n个) ,然后拒绝球体外的点。 将其余点作为vector,并将它们归一化。 这些是你的“样本” – 使用某种方法(随机,贪婪等)select其中的n
  3. 螺旋近似 。 您可以在一个球体周围绘制一个螺旋线,并将这些点均匀分布在螺旋线周围。 由于涉及math,这些比模拟更为复杂,但要快得多(可能涉及更less的代码)。 最stream行的似乎是由Saff等人 。

关于这个问题的更多信息可以在这里find

斐波那契球algorithm非常适合这一点。 它速度快,给出的结果一目了然容易欺骗人的眼睛。 您可以看到一个处理完成的例子,随着点数的增加, 这个例子会显示结果。 这是 @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个样本给你这个:

在这里输入图像描述

金色的螺旋方法

你说你不可能得到金色螺旋的方法去工作,这是一个耻辱,因为它真的很好。 我想给你一个完整的理解,以便你可以理解如何避免“聚集”。

所以这里有一个快速,非随机的方法来创build一个近似正确的点阵; 如上所述,没有格子是完美的,但这可能是“足够好”的。 它与其他方法比较,例如在BendWavy.org,但它只是有一个漂亮,漂亮的外观,以及有限的均匀间隔的保证。

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

为了理解这个algorithm,我首先邀请你看看2D向日葵螺旋algorithm。 这是基于最不合理的数字是黄金比例(1 + sqrt(5))/2的事实,并且如果通过“站在中心,转动整个转弯的黄金比例的方法发出点,然后发射另一个指向那个方向“,那么自然会形成一个螺旋,当你越来越高的点数时,它就拒绝有一个明确定义的”条“。 磁盘上均匀间隔的algorithm是,

 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) ; 我是怎么来的? 那么我希望这些在球体周围有均匀的空间 这与在大N的极限下我想要一个小区域R∈( rr + dr ), θ∈θθ + )包含许多与r d r成比例的点相同θ 。 现在假设这一个随机variables,我们这里正在讨论这个,这有一个直截了当的解释,说联合概率密度只是对一些常数c的 cr ; 单位盘上的归一化强制c = 1 /π。

这里有个技巧,它来自概率理论,它被称为抽样反CDF :假设你想产生一个概率密度为fz )的随机variables,并且有一个随机variablesU〜Uniform(0,1)就像大多数编程语言中的random() 。 你怎么做到这一点?

  1. 首先把你的密度变成一个累积分布函数 Fz ),它从0到1单调地与导数fz )相关。
  2. 然后计算CDF的反函数F -1z )。
  3. 你会发现Z = F -1U )是根据目标密度分配的。

(一旦你知道这是结果,certificate是简单的,如果你问什么是z < z < z + d z的概率,这跟问什么是z < F -1U )< z + d z ,将F应用到所有三个expression式中,注意到它是一个单调递增函数,因此Fz )< U < Fz + d z ),展开右边出来findFz )+ f z )d z ,并且由于U是一致的,所以这个概率就是fz )d z

现在,黄金比例的螺旋技巧将θ的点指向了θ的一个很好的均匀模式,所以让我们把它整合出来。 对于单位圆,我们剩下Fr )= r 2 。 所以反函数是F -1u )= u 1/2 ,因此我们将在极坐标系中在球体上产生随机点,其中r = sqrt(random()); theta = 2 * pi * random() r = sqrt(random()); theta = 2 * pi * random()

无论如何, 就是平方根的来源:而不是对这个反函数进行随机采样,我们对它进行均匀采样,关于均匀采样的好处在于,我们关于点在大N极限内展开的结果将performance为就好像我们已经随机抽样了一样。

现在做一个球体上的向日葵

我们需要做的点变化,只需要将极坐标转换为球坐标。 当然径向坐标不会进入,因为我们在单位球面上。 为了让事情在这里保持一致,即使我被训练成物理学家,我将使用math家的坐标,其中0≤φ≤π是从极点向下的纬度,0≤θ≤2π是经度。 所以与上面的区别在于我们基本上用φ代替variablesr

我们的区域元素从r d r dθ变化到不太复杂的sin( φdφdθ ,所以我们的均匀间距的联合密度是sin( φ )/4π。 对θ进行积分,得到fφ )= sin( φ )/ 2,因此Fφ )=(1-cos( φ ))/ 2。 反过来,我们可以看到一个统一的随机variables看起来像acos(1 – 2 u ),但是我们一致而不是随机地采样,所以我们改用φk = acos(1 – 2( k + 0.5)/ N )。 algorithm的其余部分只是将其投影到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,结果如下: 在这里输入图像描述 在这里输入图像描述

这个答案是基于这个答案概述的相同的“理论”

我将这个答案添加为:
– 没有其他的select符合“统一”需要“点”(或不明显)。 (注意让原始问题中的分布寻找行为在特定情况下是需要的,那么你只要从k个随机产生的点的有限列表中排除(随机selectk个项目中的索引数)。)
– 最接近的其他impl迫使你用'angular度轴'来决定'N',而在两个angular度轴上只是'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 

在低计数(2,5,7,13等N),并似乎工作“很好”

你正在寻找什么叫做球形覆盖物 。 球形覆盖问题是非常困难的,解决scheme除了less数点以外是未知的。 有一点可以肯定的是,给定一个球体上的n个点,总是存在两个距离点d = (4-csc^2(\pi n/6(n-2)))^(1/2)或更近。

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

 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维上的高斯分布是球对称的,所以在球体上的投影是均匀的。

当然,不能保证统一生成的点集合中的任何两点之间的距离将被限制在下面,因此您可以使用拒绝来强制执行任何您可能具有的条件:可能最好生成整个集合,然后必要时拒绝整个collections。 (或者使用“早期拒绝”来拒绝你迄今为止产生的整个集合;只是不要保留某些点而放弃其他点)。可以使用上面给出的d的公式,减去一些松弛,来确定最小距离之间的点之间,你会拒绝一组点。 你必须计算nselect2个距离,拒绝的可能性取决于松弛; 很难说如何,所以运行一个模拟来感受相关的统计数据。

尝试:

 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

与less量的点你可以运行一个模拟:

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

N的两个最大的因子,如果N==20那么两个最大的因子是{5,4} ,或者更一般地{a,b} 。 计算

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

把你的第一点放在{90-dlat/2,(dlong/2)-180} ,你的第二点在{90-dlat/2,(3*dlong/2)-180} {90-dlat/2,(5*dlong/2)-180} ,直到你绊倒了世界一次,那么当你走到{90-3*dlat/2,(dlong/2)-180}

很明显,我正在用球形地球的表面度数,按照惯例将+/-转换成N / S或E / W。 显然,这给你一个完全非随机的分布,但是它是统一的,并不是聚集在一起的。

为了增加某种程度的随机性,可以生成2个正态分布的(平均值为0和标准偏差为{dlat / 3,dlong / 3}),并将它们添加到你的均匀分布的点上。

编辑:这不回答OP的意思问的问题,留在这里,以防人们发现它有用。

我们使用概率的乘法规则,结合无穷小。 这导致了2行代码来实现你想要的结果:

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

(在以下坐标系中定义:)

在这里输入图像描述

你的语言通常有一个统一的随机数字原语。 例如在python中,你可以使用random.random()返回范围[0,1)中的一个数字。 你可以把这个数乘以k得到一个在[0,k)范围内的随机数。 因此在Python中, uniform([0,2pi))意味着random.random()*2*math.pi


certificate

现在我们不能统一指定θ,否则我们会在极点上聚集。 我们希望将概率与球面楔的表面面积成正比(这个图中的θ实际上是φ):

在这里输入图像描述

赤道的angular位移dφ将导致dφ* r的位移。 这个位移是在什么方位θ? 那么从Z轴开始的半径是r*sin(θ) ,所以这个“纬度”与楔相交的长度是dφ * r*sin(θ) 。 因此,我们通过积分从南极到北极的切片面积来计算从其采样的区域的累积分布 。

在这里输入图像描述 (其中stuff = dφ*r

我们现在尝试从CDF中获取相反的结果: http : //en.wikipedia.org/wiki/Inverse_transform_sampling

首先,我们通过将我们的近似CDF除以最大值来归一化。 这有消除dφ和r的副作用。

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

从而:

 let x by a random float in range [0,1] θ = -arcsin(1-2*x) 
 # 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)]) 

这工作,它是致命的简单。 尽可能多的点:

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

或…放置20点,计算二十面体的中心。 对于12点,find二十面体的顶点。 30点,二十面体的边缘的中点。 你可以用四面体,立方体,十二面体和八面体做同样的事情:一组点在顶点上,另一个在面的中心,另一个在边的中心。 但是,他们不能混合使用。

Healpix解决了一个密切相关的问题(像素等于面积像素的球体):

http://healpix.sourceforge.net/

这可能是矫枉过正的,但也许看了之后,你会意识到其他一些不错的属性是有趣的。 这不仅仅是一个输出点云的函数。

我在这里登陆试图再次find它; 名字“healpix”不完全唤起领域…