在球面上均匀分布n个点

Posted 2023-02-19

【中文标题】在球面上均匀分布n个点

【英文标题】：Evenly distributing n points on a sphere

【发布时间】：2012-03-24 22:34:53

【问题描述】：

我需要一种算法，它可以为我在一个球体周围提供 N 个点（可能小于 20 个）的位置，并将它们模糊地分散开来。不需要“完美”，但我只需要它，所以它们都不会聚集在一起。

This question 提供了很好的代码，但我找不到制作这种制服的方法，因为这似乎是 100% 随机的。 This blog post 推荐有两种方法允许输入球体上的点数，但 Saff and Kuijlaars 算法正是我可以转录的伪代码，而我发现的 code example 包含“节点 [k]”，我不能看不到解释并破坏了这种可能性。第二个博客示例是黄金分割螺旋，它给了我奇怪的、成束的结果，没有明确的方法来定义一个恒定的半径。 来自this question 的This algorithm 似乎可行，但我无法将该页面上的内容拼凑成伪代码或任何东西。

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

所以，我正在寻找的是简单的伪代码，用于在单位球体周围均匀分布 N 个点，以球坐标或笛卡尔坐标返回。如果它甚至可以稍微随机分布（想想围绕恒星的行星，分布得体，但有余地），那就更好了。

【问题讨论】：

你是什么意思“有点随机化”？您是指某种意义上的扰动吗？

OP 很困惑。 他正在寻找的是在一个球体上放置 n 个点，以使任意两点之间的最小距离尽可能大。 这将使这些点看起来“均匀分布”在整个球体。这与在球体上创建均匀随机分布完全无关，这是许多链接的意义所在，也是下面许多答案的意义所在。

如果您不希望它们看起来只是随机的，那么在球体上放置 20 个点并不是很多。

这是一种方法（它有代码示例）：pdfs.semanticscholar.org/97a6/…（看起来它使用排斥力计算）

当然，对于 4, 6, 8, 12, 20 中 N 上的值，存在精确解，其中每个点到（每个）它最近的邻居的距离对于所有人来说都是一个常数点和所有最近的邻居。

【参考方案1】：

取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，直到你环游世界一次，到那时你已经到了75,150，当你去90-3*dlat/2,(dlong/2)-180旁边。

显然，我是在球形地球表面上以度为单位工作的，通常的惯例是将 +/- 转换为 N/S 或 E/W。显然，这给了你一个完全非随机的分布，但它是均匀的，并且点没有聚集在一起。

要增加一定程度的随机性，您可以生成 2 个正态分布的（平均值为 0 和标准偏差为 dlat/3, dlong/3 视情况而定）并将它们添加到您的均匀分布点。

【讨论】：

如果你在 sin(lat) 而不是 lat 工作，那看起来会好很多。事实上，你会在两极附近聚集很多。

【参考方案2】：

在this example code 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

如果你绘制它，你会发现两极附近的垂直间距更大，因此每个点都位于大约相同的总区域空间中（靠近两极的空间更小”水平”，所以它给出更多的“垂直”）。

这与与邻居距离大致相同的所有点不同（这就是我认为您的链接所谈论的），但它可能足以满足您的需求，并且只需制作统一的纬度即可改进/lon 网格。

【讨论】：

很好，很高兴看到数学解决方案。我正在考虑使用螺旋和弧长分离。我仍然不确定如何获得最佳解决方案，这是一个有趣的问题。

您是否看到我编辑了答案以在顶部包含对 node[k] 的解释？我想这可能就是你所需要的......

太好了，感谢您的解释。稍后我会尝试一下，因为我目前没有时间，但非常感谢你帮助我。我会让你知道它最终是如何为我的目的工作的。 ^^

使用 Spiral 方法非常适合我的需求，非常感谢您的帮助和澄清。 :)

您的纬度转换为度数似乎不正确。你不应该也除以 pi 吗？

【参考方案3】：

你可以用少量的点来运行模拟：

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

【讨论】：

要改进我的答案，您应该将最接近的索引 = i 更改为最接近的索引 = randchoice(i,j)

【参考方案4】：

编辑：这并不能回答 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。

---

证明

现在我们不能统一分配 θ，否则我们会在两极处结块。我们希望分配与球形楔的表面积成正比的概率（该图中的 θ 实际上是 φ）：

赤道处的角位移 dφ 将导致 dφ*r 的位移。在任意方位角 θ 处的位移是多少？嗯，z 轴的半径是r*sin(θ)，因此与楔形相交的那个“纬度”的弧长是dφ * r*sin(θ)。因此，我们通过对从南极到北极的切片面积进行积分来计算要从中采样的区域的cumulative distribution。

（其中的东西=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)

【讨论】：

这不等同于他丢弃的“100%随机”选项吗？我的理解是，他希望它们比均匀随机分布更均匀。

@BlueRaja-DannyPflughoeft：嗯，很公平。我想我没有仔细阅读这个问题。无论如何，我把它留在这里，以防其他人觉得它有用。感谢您指出这一点。

【参考方案5】：

这被称为球体上的打包点，并且没有（已知的）通用、完美的解决方案。但是，有很多不完美的解决方案。最受欢迎的三个似乎是：

创建模拟。将每个点视为限制在球体上的电子，然后运行一定数量的步骤的模拟。电子的排斥力自然会使系统趋于更稳定的状态，在这种状态下，这些点彼此之间的距离尽可能远。 超立方体拒绝。这种听起来很花哨的方法实际上非常简单：你统一选择球体周围的立方体内部的点（远远超过n），然后拒绝球体外部的点。将剩余的点视为向量，并将它们归一化。这些是您的“样本” - 使用某种方法（随机、贪婪等）从中选择 n。 螺旋近似。您围绕球体追踪螺旋，然后将点均匀分布在螺旋周围。由于涉及数学，这些比模拟更复杂，但速度更快（并且可能涉及更少的代码）。最受欢迎的似乎是Saff, et al。

很多关于这个问题的更多信息可以找到here

【讨论】：

我将研究安德鲁·库克在下面发布的螺旋策略，但是，您能否澄清一下我想要的和“均匀随机分布”之间的区别？那只是在球体上100％随机放置点以便均匀放置吗？谢谢您的帮助。 :)

@Befall：“均匀随机分布”是指概率分布是均匀的——这意味着，当在球体上选择一个随机点时，每个点都有相同的可能性被选中。它与点的最终空间分布无关，因此与您的问题无关。

啊，好的，非常感谢。搜索我的问题会为这两个问题带来大量答案，我无法真正理解哪个对我来说毫无意义。

明确一点，每个点被选中的概率为零。该点属于球体表面上任意两个区域的概率之比等于曲面的比值。

最后一个链接现在失效了

【参考方案6】：

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

【讨论】：

如果您写一些文字来解释这是什么意思，那将会很有帮助，因此 OP 不必相信它会起作用。

【参考方案7】：

这很有效，而且非常简单。想要多少分就多少：

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

【参考方案8】：

此答案基于this answer 很好概述的相同“理论”

我将此答案添加为： - 没有其他选项适合“一致性”需要“现场”（或不明显 - 显然如此）。 （注意要在原始问题中获得特别想要的类似分布的行为，您只需从随机均匀创建的 k 个点的有限列表中拒绝（随机返回 k 个项目中的索引计数）。） - 最接近的其他 impl 迫使您通过“角轴”来决定“N”，而不是在两个角轴值上仅“N 的一个值”（在 N 的低计数下很难知道可能发生什么，或者可能无关紧要（例如，您想要“5”分——玩得开心）） --此外，很难“理解”如何在没有任何图像的情况下区分其他选项，所以这就是这个选项的样子（下图），以及随之而来的准备运行的实现。

N 为 20：

然后是 80 时的 N：

这是可以运行的 python3 代码，其中的仿真是同一来源：“http://web.archive.org/web/20120421191837/http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere”由其他人找到。 （我包含的绘图在作为“主”运行时触发，取自：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 in 2, 5, 7, 13 等）下进行测试，似乎工作得“不错”

【参考方案9】：

试试：

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

【参考方案10】：

OR... 放置 20 个点，计算二十面体面的中心。对于 12 个点，找到二十面体的顶点。对于 30 点，二十面体边缘的中点。您可以对四面体、立方体、十二面体和八面体做同样的事情：一组点位于顶点，另一组位于面的中心，另一组位于边缘的中心。但是，它们不能混合使用。

【讨论】：

一个好主意，但它只适用于 4、6、8、12、20、24 或 30 点。

如果你想作弊，你可以使用面和顶点的中心。它们将不是等间距，而是一个不错的近似值。这很好，因为它具有确定性。

【参考方案11】：

斐波那契球体算法非常适合这一点。它速度很快，并且给出的结果一目了然，很容易欺骗人眼。 You can see an example done with processing 将随着时间的推移显示添加点的结果。 Here's another great interactive example 由 @gman 制作。这是一个简单的python实现。

import math


def fibonacci_sphere(samples=1000):

    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 个样本为您提供：

【讨论】：

在定义 phi 时调用了一个变量 n：phi = ((i + rnd) % n) * increment。 n = 样本吗？

这太棒了！谢谢！！！ Here's something random I made using it, warning contains sound and uses WebGL

@Xarbrough 代码为您提供了单位球体周围的点，因此只需将每个点乘以您想要的半径标量即可。

@Fnord：我们可以为更高维度做这个吗？

真的很酷！！！您使用什么工具来生成该渲染？

【参考方案12】：

您正在寻找的东西称为球形覆盖层。球形覆盖问题非常困难，除了少量点外，解决方案未知。可以确定的一件事是，给定球体上的 n 个点，总是存在两个距离 d = (4-csc^2(\pi n/6(n-2)))^(1/2) 或更近的点。

如果您想要一种概率方法来生成均匀分布在球体上的点，这很容易：通过高斯分布在空间中均匀生成点（它内置在 Java 中，不难找到其他语言的代码）。所以在 3 维空间中，你需要类似

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 个距离，拒绝的概率将取决于 slack；很难说怎么做，所以运行一个模拟来感受一下相关的统计数据。

【讨论】：

赞成最小最大距离表达式。对于限制要使用的点数很有用。不过，参考权威来源会很好。

【参考方案13】：

黄金螺旋法

你说你无法让黄金螺旋方法发挥作用，这很遗憾，因为它真的非常非常好。我想让你对它有一个完整的了解，这样也许你就能明白如何避免它被“捆绑”。

所以这里有一个快速、非随机的方法来创建一个近似正确的格子；如上所述，没有晶格是完美的，但这可能已经足够了。它与其他方法进行比较，例如在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）

好吧，我在这里使用平方根，因为我希望它们在磁盘周围具有均匀的区域间距。就是说在大N的限制下我想要一个小区域R ∈ (r, r em> + dr), Θ ∈ (θ, θ + dθ) 包含与其面积成正比的点数，即 r dr dθ。现在，如果我们假设我们在这里谈论的是一个随机变量，这有一个简单的解释，即 (R, Θ) 的联合概率密度只是 cr 表示一些常量 c。单位圆盘上的归一化将强制 c = 1/π。

现在让我介绍一个技巧。它来自被称为sampling the inverse CDF 的概率论：假设你想生成一个概率密度为f(z)的随机变量你有一个随机变量 U ~ Uniform(0, 1)，就像大多数编程语言中的 random() 一样。你是怎么做到的？

首先，将您的密度转换为cumulative distribution function 或 CDF，我们将其称为 F(z)。请记住，CDF 通过导数 f(z) 从 0 单调增加到 1。 然后计算CDF的反函数F^-1(z)。 你会发现Z = F^-1(U)是按照目标密度分布的. ^{（注 3）。}

现在黄金比例螺旋技巧以一个非常均匀的模式将点分隔开来，θ 所以让我们把它整合起来；对于单位圆盘，我们剩下 F(r) = r²。所以反函数是F^-1(u) = u^1/2，因此我们将在磁盘上生成带有r = sqrt(random()); theta = 2 * pi * random() 的极坐标随机点。

现在不是随机对这个反函数进行采样，而是均匀对它进行采样，而均匀采样的好处在于我们的结果是关于点如何分布在大 N 的限制将表现得就像我们随机采样它一样。这种组合就是诀窍。我们使用(arange(0, num_pts, dtype=float) + 0.5)/num_pts 而不是random()，因此，如果我们要采样10 个点，它们是r = 0.05, 0.15, 0.25, ... 0.95。我们统一采样 r 以获得等面积间距，我们使用向日葵增量来避免输出中出现可怕的点“条”。

现在在球体上做向日葵

我们需要做的改变是用点来点缀球体，只需要将极坐标换成球坐标。径向坐标当然不包括在内，因为我们在一个单位球面上。为了在这里保持一致，即使我是受过物理学训练的，我也会使用数学家的坐标，其中 0 ≤ φ ≤ π 是从极点向下的纬度，0 ≤ θ ≤ 2π 为经度。所以与上面的不同之处在于，我们基本上是将变量 r 替换为 φ。

我们的面积元素，以前是 r dr dθ，现在变成了不太复杂的 sin(φ) dφ dθ。所以我们均匀间距的联合密度是 sin(φ)/4π。积分出θ，我们发现f(φ) = sin(φ)/2，因此 F(φ) = (1 − cos(φ))/2。反过来，我们可以看到均匀随机变量看起来像 acos(1 - 2 u)，但是我们均匀采样而不是随机采样，所以我们使用 φ_k = acos(1 − 2 (k + 0.5)/N)。算法的其余部分只是将其投影到 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，结果如下所示：

进一步研究

我想对 Martin Roberts 的博客大声疾呼。请注意，上面我通过向每个索引添加 0.5 创建了索引的偏移量。这对我来说只是视觉上的吸引力，但it turns out that the choice of offset matters a lot 在整个间隔内并不是恒定的，如果选择正确，这意味着包装的准确度可以提高 8%。还应该有一种方法可以让his R₂ sequence 覆盖一个球体，看看这是否也产生了一个很好的均匀覆盖，也许是原样但可能需要，比如说，只取自一半单位正方形沿对角线左右切割并拉伸得到一个圆形。

注意事项

那些“条”是由一个数的有理逼近形成的，一个数的最佳有理逼近来自它的连分数表达式，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 之间，因此连分数中的大整数可以提供特别好的有理近似值：“一个除以 100 到 101 之间的某物”比“一个除以某物”好1到2之间。”因此，最无理数是1 + 1/(1 + 1/(1 + ...)) 并且没有特别好的有理近似值；可以通过乘以φ得到φ = 1 + 1/φ，得到黄金比例的公式。

对于不太熟悉 NumPy 的人来说——所有的函数都是“矢量化的”，所以 sqrt(array) 与其他语言可能写的 map(sqrt, array) 相同。所以这是一个逐个组件的sqrt 应用程序。这同样适用于标量除法或标量加法 - 并行适用于所有组件。

一旦你知道这是结果，证明就很简单了。如果你问 z Z z + dz 的概率是多少，这和问是一样的z F^-1(U) z + 的概率是多少dz，将 F 应用于所有三个表达式，注意它是一个单调递增函数，因此 F(z ) U F(z + dz)，右手边向外展开找到F(z) + f(z) dz，并且由于 U 是一致的，这个概率只是 f(z) dz 正如所承诺的。

【讨论】：

我不知道为什么这么差，这是迄今为止最好的快速方法。

@snb 谢谢你的好话！它之所以如此落后，部分原因是它比这里的所有其他答案都年轻得多。我很惊讶它的表现甚至和以前一样好。

@FelixD。这听起来像是一个很快就会变得非常复杂的问题，特别是如果你开始使用大圆距离而不是欧几里得距离。但也许我可以回答一个简单的问题，如果将球体上的点转换为他们的 Voronoi 图，则可以将每个 Voronoi 单元描述为具有大约 4π/N 的面积，并且可以通过假装它是一个圆来将其转换为特征距离。比菱形，πr² = 4π/N。那么r=2/√(N)。

使用具有实际均匀输入而不是随机均匀输入的采样定理是让我说 “好吧，为什么 #$%& 我没有想到这一点？ ".不错。

好问题！我相信我的答案更接近“它起作用的原因”，而 Martin 的答案更加精确。因此，根据定义，黄金比例满足 φ² = φ + 1，它重新排列为 φ – 1 = 1/φ，乘以 2 π，前导数字 1 刚刚被三角函数击中。所以在浮点数中，只需减去一个就会用 0 填充第 53 位，其中 1 会更正确。

【参考方案14】：

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

http://healpix.sourceforge.net/

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

我降落在这里试图再次找到它； “healpix”这个名字并不能完全唤起球体……

【参考方案15】：

@robert king 这是一个非常好的解决方案，但其中有一些草率的错误。我知道它对我有很大帮助，所以不要介意马虎。 :) 这是一个清理后的版本....

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

【参考方案16】：

根据 fnord 的回答，这是一个 Unity3D 版本，增加了范围：

代码：

// golden angle in radians
static float Phi = Mathf.PI * ( 3f - Mathf.Sqrt( 5f ) );
static float Pi2 = Mathf.PI * 2;

public static Vector3 Point( float radius , int index , int total , float min = 0f, float max = 1f , float angleStartDeg = 0f, float angleRangeDeg = 360 )

    // y goes from min (-) to max (+)
    var y = ( ( index / ( total - 1f ) ) * ( max - min ) + min ) * 2f - 1f;

    // golden angle increment
    var theta = Phi * index ; 
        
    if( angleStartDeg != 0 || angleRangeDeg != 360 )
    
        theta = ( theta % ( Pi2 ) ) ;
        theta = theta < 0 ? theta + Pi2 : theta ;
            
        var a1 = angleStartDeg * Mathf.Deg2Rad;
        var a2 = angleRangeDeg * Mathf.Deg2Rad;
            
        theta = theta * a2 / Pi2 + a1;
    

    // https://***.com/a/26127012/2496170
    
    // radius at y
    var rY = Mathf.Sqrt( 1 - y * y ); 
    
    var x = Mathf.Cos( theta ) * rY;
    var z = Mathf.Sin( theta ) * rY;

    return  new Vector3( x, y, z ) * radius;

要点：https://gist.github.com/nukadelic/7449f0872f708065bc1afeb19df666f7/edit

预览：