黄金螺旋法
你说你无法让黄金螺旋方法发挥作用,这很遗憾,因为它真的非常非常好。我想让你对它有一个完整的了解,这样也许你就能明白如何避免它被“捆绑”。</p>
因此,这是一种快速、非随机的方法来创建近似正确的晶格;如上所述,没有晶格是完美的,但这可能已经足够了。它与其他方法进行了比较,例如在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 + d r ), Θ ∈ ( θ , θ + d θ ) 包含与其面积成比例的点数,这是r d r d θ。现在,如果我们假设我们在这里谈论的是一个随机变量,这有一个简单的解释,即 ( R , Θ ) 的联合概率密度就是cr对于一些常数c。单位圆盘上的归一化将迫使c = 1/π。
现在让我介绍一个技巧。它来自概率论,它被称为对逆 CDF 进行采样:假设您想生成一个概率密度为f ( z ) 的随机变量,并且您有一个随机变量U ~ Uniform(0, 1),就像从random()
在大多数编程语言中。你怎么做到这一点?
- 首先,将您的密度转换为累积分布函数或 CDF,我们将其称为F ( z )。记住,CDF 通过导数f ( z )从 0 单调增加到 1 。
- 然后计算 CDF 的反函数F -1 ( z )。
- 您会发现Z = F -1 ( U ) 是根据目标密度分布的。(注 3)。
现在黄金比例螺旋技巧将这些点以一个非常均匀的模式隔开,让我们把它整合起来;对于单位圆盘,我们剩下F ( r ) = r 2。所以反函数是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 θ,现在变成了不太复杂的 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 创建了索引的偏移量。这只是在视觉上吸引了我,但事实证明,偏移量的选择很重要,并且在间隔内不是恒定的,如果选择正确,这意味着打包的准确度可以提高 8%。还应该有一种方法可以让他的 R 2序列覆盖一个球体,看看这是否也产生了一个很好的均匀覆盖,也许是原样但可能需要,比如说,只取自一半单位正方形对角线左右切割并拉伸得到一个圆圈。
笔记
这些“条”是由一个数字的有理逼近形成的,一个数字的最佳有理逼近来自其连分数表达式,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 + 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所承诺的。