黄金螺旋法
你说你无法让黄金螺旋方法发挥作用，这很遗憾，因为它真的非常非常好。我想让你对它有一个完整的了解，这样也许你就能明白如何避免它被“捆绑”。
所以这里有一个快速、非随机的方法来创建一个近似正确的格子；如上所述，没有晶格是完美的，但这可能已经足够了。它与其他方法进行比较，例如在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) = r2。所以反函数是F-1(u) = u1/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 R2 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 正如所承诺的。