编辑自太极图 　　

　　

　　量子比特|微信官方账号QbitAI 　　

　　

　　由随机均匀的点组成的图案在动物和植物中很常见。 　　

　　

　　杨梅、草莓、荔枝、红毛丹等水果表面都有颗粒状或毛状结构，随机均匀地散布在水果表面： 　　

　　

　　 　　

　　

　　类似的模式在动物身上也有，比如人人都爱涮的毛肚： 　　

　　

　　 　　

　　

　　同样，在计算机模拟下，有许多场景需要在空间中随机均匀地生成点。 　　

　　

　　比如动物毛发生成时毛孔的位置，多人游戏中玩家出生的位置，森林生成时树木的位置等等。 　　

　　

　　这些场景的共同特点是任何两点之间的距离应该是大于等于一个下界(这个下限是预设的，可以通过改变它来控制生成点之间的间隔)。 　　

　　

　　但是，如果直接使用完全随机生成的点数，很有可能会得到非常不均匀的分布结果，有些地方“聚集”，有些地方稀疏： 　　

　　

　　 　　

　　

　　如果用这样的点来模拟头发等位置的生成，效果会很差。 　　

　　

　　因此，需要在生成点的过程中加入一个距离判断来剔除那些不符合要求的点。 　　

　　

　　以前，使用NumPy产生这样的效果经常需要70s,这是非常不经济的。 　　

　　

　　 　　

　　

　　现在，太极图形已经实现了一种基于Taichi,的超快速算法，在单个CPU线程上运行也有同样的效果。只需要0.7s就能产生这样的模式，其速度和100倍.差不多 　　

　　

　　让我们看看他们是怎么做的。 　　

　　

　　在Bridson算法实现之前，一种常见的算法是dart throwing(就像一个被蒙住眼睛的人随意扔飞镖): 　　

　　

　　在每个区域内随机选取一个点，检查该点与所有获得的点是否存在“冲突”。 　　

　　

　　如果这个点和一个获得的点之间的最小距离小于指定的下限，则丢弃这个点，否则它是一个合格的点，并将其添加到现有的点集。 　　

　　

　　重复此操作，直到获得足够的点数，或者连续失败n次(n为设定的正整数)。 　　

　　

　　但是这个算法的效率很低。 　　

　　

　　随着获得点数的增加，冲突的概率越来越大，获得新点数所需的时间也越来越长。每次比较当前点和所有已有点的距离也会降低效率。 　　

　　

　　相比之下，Bridson算法的效率更高。 　　

　　

　　这个算法的原理来自于Robert Bridson在2007年发表的论文《仲裁维中的快速泊松圆盘采样》1(论文很短，只有一张A4纸)。如果去掉标题和介绍，真正的算法内容只是一小段。 　　

　　

　　开始时的这张移动图片演示了Bridson disk采样算法在400x400.的一个网格区域上的运行效果。该算法试图获得100K个均匀分散的点，实际上产生了大约53.7K: 　　

　　

　　 　　

　　

　　这个动画由Taichi生成，运行在一个CPU线程上。除去编译的时间计算，只需要0.7s多一点，而同样的代码翻译成《NumPy. 2》需要70s左右 　　

　　

　　从上面的动画效果可以看出，Bridson算法很像白菜的生长过程：我们从一个种子点开始，一层一层的添加新点。 　　

　　

　　每次我们添加一个新点，它都位于最外层的点的周围，并尽可能覆盖最外层。 　　

　　

　　为了避免每次都检查所有现有点的距离，太极采用了所谓的网格化:技术 　　

　　

　　将整个空间划分成网格。对于一个需要检查的点，只需找到它的网格，然后检查它与相邻网格中的点之间的最小距离。 　　

　　

　　只要这个距离大于 　　

指定的下界，更远处的点就不必再检查了。这个技巧在图形学和物理仿真中是非常常用的。

　　

这个采样过程很难并行化，因为当一个线程“偷偷”加入一个新的点的时候，会改变其它所有线程对距离的判断。所以Taichi仅使用单个CPU线程来执行这个算法：

　　

ti.init(arch=ti.cpu)上面的代码中通过指定arch=ti.cpu来让程序运行在CPU上。

　　

你可能会想，既然是单线程+CPU，那为什么不直接写纯Python呢？别着急，我们的计算部分会放在ti.kernel函数中，这种函数并不运行在Python虚拟机中，而是会被Taichi编译执行，所以会比纯Python的实现快很多倍！

　　

在我们介绍Bridson算法的具体实现之前，你不妨猜猜这个Taichi程序包含多少行代码？

　　

安装和导入Taichi首先推荐大家使用最新的Taichi发布版本，这样可以使用更丰富的功能，在不同平台上的支持也更稳定。截止本文写作时最新版本是1.0.3:

　　

pip install taichi==1.0.3然后，在代码开头写上：

　　

import taichi as tiimport taichi.math as tm这样会导入Taichi以及Taichi的math模块。math模块除了包含常用的数学函数之外，还提供了非常方便的向量运算。

　　

准备工作在泊松采样算法中，采样点之间的距离有一个下界r。

　　

我们假设整个区域是由N×N个同样大小的方格组成的网格区域，使得每个小方格的对角线长度正好是r，即网格的边长是r/√2。

　　

于是任何小方格中至多包含一个点。如下图所示：

　　

　　

这就是我们前面提到的网格化方法，即对于任何一个点p，设它所在的方格是D，则任何与p的距离小于等于r的点必然位于以D中心的、由5×5个方格组成的正方形区域中。

　　

在检查距离时，我们只要针对这个子区域进行计算即可。

　　

我们用一个一维数组samples和一个N×N的二维数组grid来记录已经得到的采样点：

　　

samples保存当前所有已经采样点的坐标，它的每个元素是一个二维坐标（x,y）。grid是一个整数，它存储的是第（i, j）个方格中采样点在数组samples中的下标。grid = -1表示这个方格中没有采样点。于是我们的初始设置可以这样写：

　　

grid_n = 400res = (grid_n, grid_n)dx = 1 / res<0>inv_dx = res<0>radius = dx * ti.sqrt(2)desired_samples = 100000grid = ti.field(dtype=int, shape=res)samples = ti.Vector.field(2, float, shape=desired_samples)这里网格大小设置为400x400，它占据的平面区域是<0,1>×<0,1>，所以网格的步长是dx = 1/400。采样的最小间隔是每个小方格对角线的长度，即radius = sqrt(2)*dx。

　　

我们把采样点的目标个数设置为desired_examples=100000，这是一个目测值，因为400x400的网格包含160000个小方格，考虑到每个方格中至多只有一个点，我们能得到的满足距离约束的点的最大数目肯定少于160000。

　　

初始时网格中没有任何点，所以需要将grid中的值都置为-1:

　　

grid.fill(-1)如何生成新的点在加入新的随机点时，我们总是从已有点的附近随机选择一个位置，然后比较它和已知点是否满足最小距离约束，是的话就将其加入已有点，否则就将其抛弃然后重新选择点。

　　

这里需要注意的是：

　　

当一个已有点附近已经被填满时，我们后面再加入新的点时就不必考虑它的附近了，可以用一个下标head来记录这一点。我们约定samples数组中下标< head的点附近都已经被填满，从而不必再考虑它们，只考虑下标>= head的点即可。初始时head = 0。samples是一个长度为100K的数组，这不代表我们真的能取到这么多点，但具体个数是多少无法事先确定，所以我们还需要用一个下标tail来记录目前已经获得的点的个数。初始时tail = 1，因为我们将选择区域中心作为第一个点。当然这个初始点的位置可以是任意的。正如前面提到的，当我们检查一个点p是否与已有点满足最小距离约束时，没有必要遍历检查所有的点。只要检查以p所在方格为中心，由5×5个方格组成的正方形区域即可。检查一个点是否和已有点冲突的逻辑我们单独写成一个函数：

　　

@ti.funcdef check_collision(p, index): x, y = index collision = False for i in range(max(0, x - 2), min(grid_n, x + 3)): for j in range(max(0, y - 2), min(grid_n, y + 3)): if grid != -1: q = samples> if (q - p).norm() < radius - 1e-6: collision = True return collision其中p是需要检查点的坐标，index=(x, y)是p所在的方格的下标。

　　

我们遍历所有满足x-2 <= i <= x+2和y-2 <= j <= y+2的下标(i, j)，检查方格(i, j)中是否已经有点，即 grid是否等于-1。有的话它与p的距离是否小于radius，然后返回对应的判断。

　　

完成了准备工作，我们可以开始正式的循环了：

　　

@ti.kerneldef poisson_disk_sample(desired_samples: int) -> int: samples<0> = tm.vec2(0.5) grid = 0 head, tail = 0, 1 while head < tail and head < desired_samples: source_x = samples head += 1 for _ in range(100): theta = ti.random() * 2 * tm.pi offset = tm.vec2(tm.cos(theta), tm.sin(theta)) * (1 + ti.random()) * radius new_x = source_x + offset new_index = int(new_x * inv_dx) if 0 <= new_x<0> < 1 and 0 <= new_x<1> < 1: collision = check_collision(new_x, new_index) if not collision and tail < desired_samples: samples = new_x grid = tail tail += 1 return tail首先我们把区域的中心，即坐标为(0.5,0.5)的点选择为初始点，让它作为“种子”将随机点逐渐扩散到整个区域。

　　

接下来的while循环是算法的主体，这个循环是串行执行的，只占用一个线程。

　　

我们每次找到第一个需要考虑的点samples，然后在以它为中心，半径为的圆环中随机选择100个点，逐个检查这100个点是否不超出<0,1>×<0,1>的区域范围，以及是否和已有点冲突。

　　

如果都满足的话它就是一个合格的点，我们将它的坐标和方格下标更新到samples和grid中，并将已有点的个数tail增加1。在这100个点都检查完后，可能有多个点会被加入已有点的集合。

　　

注意在半径为的圆环中采样可以让我们得到的点在满足最小距离约束的同时距离已有点也不会太远。

　　

当这100个点都检查完毕后，我们可以认为samples这个点的周围已经没有空白区域可以放置新的点，所以将head增加1，并重新检查下一个samples 附近的区域。

　　

当所有的点周围的空间都已经被填满，即head = tail时，或者我们已经获得了desired_samples个点，即tail = desired_samples时循环结束。这时samples中下标在0~tail-1范围内的点就是全部的已有点。

　　

展示动画效果我们可以只用几行代码，就把整个采样的过程用动画的形式显示出来：

　　

num_samples = poisson_disk_sample(desired_samples)gui = ti.GUI("Poisson Disk Sampling", res=800, background_color=0xFFFFFF)count = 0speed = 300while gui.running: gui.circles(samples.to_numpy()<:min(count * speed, num_samples)>, color=0x000000, radius=1.5) count += 1 gui.show()这里我们控制动画的速度为每生成300个点就绘制一帧。

　　

至此我们已经介绍完了程序的所有要点，把各部分组合起来：

　　

import taichi as tiimport taichi.math as tmti.init(arch=ti.cpu)grid_n = 400res = (grid_n, grid_n)dx = 1 / res<0>inv_dx = res<0>radius = dx * ti.sqrt(2)desired_samples = 100000grid = ti.field(dtype=int, shape=res)samples = ti.Vector.field(2, float, shape=desired_samples)grid.fill(-1)@ti.funcdef check_collision(p, index): x, y = index collision = False for i in range(max(0, x - 2), min(grid_n, x + 3)): for j in range(max(0, y - 2), min(grid_n, y + 3)): if grid != -1: q = samples> if (q - p).norm() < radius - 1e-6: collision = True return collision@ti.kerneldef poisson_disk_sample(desired_samples: int) -> int: samples<0> = tm.vec2(0.5) grid = 0 head, tail = 0, 1 while head < tail and head < desired_samples: source_x = samples head += 1 for _ in range(100): theta = ti.random() * 2 * tm.pi offset = tm.vec2(tm.cos(theta), tm.sin(theta)) * (1 + ti.random()) * radius new_x = source_x + offset new_index = int(new_x * inv_dx) if 0 <= new_x<0> < 1 and 0 <= new_x<1> < 1: collision = check_collision(new_x, new_index) if not collision and tail < desired_samples: samples = new_x grid = tail tail += 1 return tailnum_samples = poisson_disk_sample(desired_samples)gui = ti.GUI("Poisson Disk Sampling", res=800, background_color=0xFFFFFF)count = 0speed = 300while gui.running: gui.circles(samples.to_numpy()<:min(count * speed, num_samples)>, color=0x000000, radius=1.5) count += 1 gui.show()代码总行数：60。

　　

One More Thing具体来说，这篇代码实现了两个操作：

　　

在60行代码中实现了一个完整的泊松采样动画。在一个400x400的网格中采集了53k个点，但耗时不到1秒。相关代码可以在文末的原文链接中找到。

　　

严格来说，本文实现的算法和Bridson论文里描述的算法有一点点不一样的地方（更简单一些），但是效果却差不多。

　　

你能看出是哪里不一样吗？（TIP：可以关注一下原论文Step 2中“active list”的处理方式）

　　

项目地址：
https://github.com/taichi-dev/poisson-sampling-homework

　　

参考资料：
<1>Robert Bridson的原论文见Fast Poisson Disk Sampling in Arbitrary Dimensions.
<2>Poisson采样用Taichi, Numpy, Numba实现的benchmark比较见GitHub

　　

― 完 ―

　　

量子位 QbitAI 头条号签约

　　

关注我们，第一时间获知前沿科技动态