泊松盘采样（Poisson Disk Sampling）与Bridson算法

2024年，九位数学家团队发表了一份近 **1000页** 的巨著，证明了 **几何Langlands猜想**（geometric Langlands conjecture）。这是纯数学领域的巅峰成就之一，而我早已接受一个事实：即便他们证明的命题本身，我也永远无法理解，更不用说那份证明了。 而在光谱的另一端，2007年Robert Bridson发表了一篇仅 **一页** 的论文，被引用了近 **1000次**，且全文不到10分钟就能完全读懂。它针对计算机图形学和仿真中常见的一类问题——**随机放置物体，但彼此不能靠得太近**——给出了一个简洁优雅的解决方案。 比如，你想程序化生成一片森林，需要一种放置树木的方法。直接用均匀随机采样显然行不通： 有些树会重叠在一起！我们需要确保任意两棵树之间保持最小距离。满足这个约束的树木分布称为 **Poisson disk分布**。我们可以尝试一种朴素的拒绝采样方法：随机"扔飞镖"，如果某个点与其他已有点的距离小于最小阈值就丢弃。但如果没有更高效的数据结构，每次采样都需要线性时间检查碰撞，而拒绝率会迅速逼近100%。Bridson的算法则给出了一种高效的实现方式。

Bridson算法

假设点之间的最小距离为 $r$，我们在 $d$ 维空间中操作。Bridson算法步骤如下：

1. 将空间划分为边长为 $\frac{r}{\sqrt{d}}$ 的网格。这样可以保证每个网格单元内最多有一个点。
2. 初始化活跃列表 active，其中包含从空间中均匀随机选取的一个点。
3. 当 active 非空时，重复以下操作：
4. 从 active 中均匀随机选择一个点 $p$。
5. 在 $p$ 为中心、内半径 $r$、外半径 $2r$ 的环形区域内，最多尝试采样 $k$ 次。若找到一个合法的Poisson disk样本（利用网格进行高效的碰撞检测），则将其加入 active，并重新选择新的 $p$。若在 $k$ 次尝试内未找到合法点，则将 $p$ 从 active 中移除。Bridson建议取 $k = 30$。

均匀采样环形区域最简单的方法是：生成一个随机单位向量 $\vec{v} \in \mathbb{R}^d$，再在区间 $\left[\frac{1}{2^d}, 1\right)$ 上均匀抽取一个数 $x$，最终样本即为 $2rx^{1/d}\cdot\vec{v}$。几年前我制作了一个视频解释了这一做法的原理。在二维情况下，选取单位向量等价于选取角度 $\theta \in [0, 2\pi)$；在更高维度中，只需从正态分布中采样各分量然后归一化即可。

改进方案

我发现对 Bridson 算法有两处简单的改进，可以大幅减少生成同等数量点所需的迭代次数。第一处仅适用于二维，第二处则可推广至高维。

先从二维改进说起。考虑算法放置了点 $p$ 后，在其环形区域内采样得到新点 $q$ 的情形。我们将 $p$ 称为 $q$ 的父节点。这两点之间的关系蕴含着重要信息：当我们不可避免地采样以 $q$ 为中心的环形区域时，存在一整段角度范围可以跳过——落在其中的点会与 $p$ 过近。这段范围在下一张图中用虚线标出。

|p - q| = 1.50r

虽然直观理解不难，但要把它转化为公式就是一场繁琐的三角学练习。以下略去细节，直接断言虚线所围成的锥以角度 $\alpha$ 为中心、张角为 $2\beta$，其中

$$\begin{align*} \alpha &= \operatorname{atan2}(p_y - q_y, p_x - q_x), \\ \beta &= \min\left(\arccos\frac{|p-q|^2+3r^2}{4r \cdot |p-q|}, \arccos\frac{|p-q|}{2r}\right). \end{align*}$$

这个公式中唯一值得注意的部分是 $\beta$ 表达式里的取最小值操作。它考虑了这样一个事实：环形带的内圆还是外圆来界定圆锥，取决于 $p$ 与 $q$ 之间的距离。我们必须取最小值，才能保证圆锥与环形带的交点完全落在圆内。当两点间距离跨越 $\sqrt{3} \cdot r$ 时，圆锥的边界点会从外圆突变到内圆，这一点很值得留意。最小值中的第一项是与外圆的相交角，第二项则是与内圆的相交角。

实现这一优化只需记录每个点的父节点即可。之后就能算出圆锥的角度范围，并在圆锥之外的范围内生成下一个采样点的角度 $\theta$。

下图展示了在边长 $\ell = 100$、$r = 1$ 的网格上，Bridson 算法在有和无父节点优化的情况下分别生成了多少个点。每个数据点是 100 次试验的平均值。

这种优化很可能可以推广到更高维度，但这需要更多空间来存储环形带的接触向量，而且收益可能会递减，因为在高维空间中，环形带与球体的交集体积占比会趋于 insignificant。同样地，我们也可以在存储父节点的同时存储子节点，以消除更多的环形带区域，但这也会增加存储开销，并且会使选择 $\theta$ 的过程慢得多。

另一种改进不改变下一采样点的角度选取方式，而是改变距离的选取策略。考虑环形带中各点到其中心的距离分布，其累积分布函数（CDF）在区间 $[r, 2r]$ 上与 $x^d$ 成正比。关于这一点，我之前在视频中已有说明。但如果我们将指数改为某个不同于 $d$ 的常数 $c$ 呢？这样我们就可以让点更靠近或更远离中心。对于 $c \neq 0$，精确的 CDF 为

$$F_c(x) = \begin{cases} 0 & x \leq r, \\ \frac{x^c - r^c}{(2^c - 1)r^c} & r < x \leq 2r, \\ 1 & x > 2r. \end{cases}$$

下图展示了随着 $c$ 变化，圆环区域内的 CDF 和 500 个随机样本呈现怎样的形态。记住，$c = 2$ 对应均匀分布。

c = 2.00

你可能会注意到滑块允许设置 $c = 0$，尽管 $F_0(x)$ 因除以零而无定义。为此，我们定义 $F_0(x) = \lim_{c \to 0} F_c(x)$，即

$$F_0(x) = \begin{cases} 0 & x \leq r, \\ \log_2{x} - \log_2{r} & r < x \leq 2r, \\ 1 & x > 2r. \end{cases}$$

为了用任意 $c$ 值对半径进行采样，我们可以采用逆变换采样法。当 $c = 0$ 时，半径取 $r \cdot 2^x$，其中 $x$ 为 $[0, 1)$ 上的均匀随机变量。否则，半径取 $2ry^{1/c}$，其中 $y$ 为 1 到 $\frac{1}{2^c}$ 之间的均匀随机变量。区间的上下界会因 $c$ 的正负而互换。

下面的热力图展示了 $c$ 对 Bridson 算法生成点数量的影响。实验设置与之前相同：在边长 $\ell = 100$、$r = 1$ 的网格上进行 100 次试验。

这似乎暗示我们应该将 $c$ 设为一个绝对值很大的负数，甚至趋近负无穷，使每个点到其父节点的距离恰好为 $r$。虽然这确实能最大化生成的点数、形成更紧密的排列，但代价是分布的"随机感"大打折扣。极端情况下（$c = -\infty$），会出现许多瑕疵，如点连成长串、距离限制无法覆盖的空隙等。事后也很容易还原出点是如何逐个生成的树结构。

因此，我们需要在最大化点密度与保留随机性之间取得平衡。当 $k$ 固定于 $[15, 40]$ 区间内时，经验表明配合父节点优化使用 $c = -1.4 - \frac{17}{\sqrt{k}}$ 效果最佳。该公式旨在使生成的点数大致匹配均匀最大泊松盘采样器的预期输出（下文详述）。对于 15 至 40 之间的每个 $k$ 值，我都通过二分搜索确定了 $c$，使得以每个点为中心、半径为 $r/2$ 所画的圆，总面积恰好占据整体的 54.7%。这一比例正是 随机序贯吸附 模型下圆形盘的饱和覆盖度。当然，该结论仅适用于二维空间，高维场景下的 $c$ 参数需另行调优。

点绘技术

此前我们一直将 $r$ 视为常量，但这并非必须。我们可以依据函数 $r\colon \mathbb{R}^d \to \mathbb{R}$ 动态设定点间的最小距离。因此，每放置一个点 $p$ 后，便在内径为 $r(p)$ 的环形区域内进行采样。一个有趣的应用是将 $r$ 映射为图像各像素的亮度值，从而生成点绘效果。左图展示了黑白图像的处理结果，右图则叠加了三组泊松盘采样样本（分别对应三个颜色通道）。 Birdson 算法本质上是串行的，但市面上也有专为并行设计的替代算法，可带来巨大的性能提升。其中我最喜欢的是 PixelPie，它能完全在 GPU 上运行。我曾利用该算法开发出了 Poisson Cam——一款基于泊松盘采样实现实时视频点绘的工具。这是我最喜欢的开发项目之一，因为它让我掌握了着色器编程、Rust、流压缩算法，以及 PixelPie 算法本身。

最大性与均匀性

2022 年，Scott A. Mitchell 发表了一篇 非常精彩的论文，提出了一种在二维空间生成泊松盘采样点的全新方法。据我所知，该方法自发表以来并未引起太多关注。本文的最后一部分将专门介绍这一成果。 在介绍算法之前，我想先说明三点，这也是本算法优于 Bridson 算法的地方：

1. **极大性**。算法结束后，保证无法在不违反泊松盘性质的前提下再放入任何点。
2. **均匀性**。算法对所有极大泊松盘样本集合均匀采样。
3. **确定性**。算法不依赖拒绝采样，因此不存在放置失败的尝试。

极大性与均匀性已被多种算法实现，其中最受欢迎的是**层次飞镖投掷法**，但 Mitchell 的算法是第一个在不依赖任何拒绝采样的情况下做到这一点的，同时性能也很出色。我虽然没有做过严谨的基准测试，但发现 [Mitchell 的实现](https://github.com/sandialabs/DeterministicMPS/tree/main) 运行速度与我在 $k = 20$ 和 $c = -5.2$ 下实现的 Bridson 算法（含父节点优化）大致相当，生成的点数也相近。据我所知，该算法唯一的缺点是复杂度显著更高。以下是大幅简化后的算法摘要：

1. 将空间划分为边长为 $\frac{r}{\sqrt{d}}$ 的网格。
2. 当仍有空间可放置新点时，循环执行：
3. 按剩余面积加权随机选取网格中的一个单元 $c$。
4. 将 $c$ 分解为互不相交的三角形与楔形区域。¹
5. 按面积加权随机选取一个三角形或楔形区域 $t$。
6. 在 $t$ 上均匀采样一点 $p$，并将其加入最终集合。
7. 从网格中挖去以 $p$ 为圆心、半径为 $r$ 的圆形区域。

原始论文已对该算法做了非常详尽的说明，因此我最能提供的帮助就是下面的可视化示意：

脚注

1. 楔形区域（chock）是一种由圆弧、径向射线与切线围成的三边形状。↩