本文经AI新媒体量子位授权转载(公众号ID:QbitAI),转载请联系出处用于转载。随机、均匀的点图案在植物和动物中已经很常见。杨梅、草莓、荔枝、红毛丹等水果,表面有颗粒状或毛发状结构,随机均匀分布在水果表面:类似的图案也存在于动物身上,比如大家喜爱的毛肚torinse:同样,在计算机模拟下,有很多场景需要在空间上随机均匀地生成点。像生成动物毛发时毛孔的位置,多人游戏中玩家刷怪的位置,生成森林时树木的位置等等。这些场景的共同特点是任意两点之间的距离需要大于或等于等于一个下界(这个下界是预设的,改变它可以控制生成点之间的间隔)。但是,如果直接使用完全随机生成的点,很有可能会得到分布很不均匀的结果。有的地方“成簇”,有的地方稀疏:如果用这样的点来模拟头发等位置的生成,效果会很差。因此需要在生成点的过程中加入距离判断,剔除那些不符合要求的点。以前,用NumPy生成这样的效果往往需要70秒左右,非常不经济。现在,TaichiGraphics已经实现了基于Taichi的超快算法。同样的效果在单CPU线程上运行,生成这样一个pattern只需要0.7s,大约快了100倍。让我们看看他们是怎么做到的。使用Bridson算法实现以前有一种常见的算法叫做投飞镖(就像一个人被蒙上眼睛随机扔飞镖):每次在区域内随机选择一个点,并检查该点与所有获得的点之间的距离。是否有“冲突”。如果该点与一个已经得到的点的最小距离小于指定的下界,则丢弃该点,否则为合格点,将其加入到已有点集合中。重复此操作,直到获得足够的分数,或连续失败N次(N为设定的正整数)。但是这个算法效率很低。因为随着获取点数的增加,冲突的概率增加,获取新点所需的时间越来越长,而且每次比较当前点与所有已有点的距离也会降低效率。.相比之下,Bridson算法效率更高。该算法的原理来自RobertBridson2007年发表的论文《FastPoissonDiskSamplinginArbitraryDimensions》[1](论文很短,只有一页A4纸)。如果去掉标题和介绍,真正的算法内容只是一小段。一开始的动画演示了Bridson圆盘采样算法在400x400的网格区域上的运行效果。算法尝试获取100K个均匀散布点,实际生成约53.7K:此动画使用太极生成,单CPU线程运行,除去编译时间计算,耗时仅0.7s多一点,而将相同的代码翻译成NumPy大约需要70秒。[2]从上面的动画可以看出,Bridson算法非常类似于卷心菜的生长过程:我们从一个种子点开始,逐层添加新的点。每次我们添加一个新的点,它位于最外面的点周围,并尽可能地包裹最外层。为了避免每次都去检查所有存在点之间的距离,Taichi采用了所谓的网格技术:将整个空间划分为网格,对于一个需要检查的点,只要找到它所在的网格即可,然后检查它与相邻网格中的一个点之间的最小距离是否足够。只要此距离大于指定的下限,就不必再次检查其他点。这个技巧在图形和物理模拟中很常见。这个采样过程很难并行化,因为当一个线程“偷偷”增加一个新的点时,会改变所有其他线程对距离的判断。所以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:pipinstalltaichi==1.0.3然后,在代码的开头写上:importtaichiastiimporttaichi.mathastm这样就会导入Taichi和Taichi的数学模块。除了常用的数学函数,数学模块还提供了非常方便的向量运算。准备工作在泊松采样算法中,采样点之间的距离有一个下界r。我们假设整个区域是由N×N个大小相同的正方形组成的网格区域,使得每个小正方形的对角线长度正好为r,即网格的边长为r/√2。所以任何小正方形最多包含一个点。如下图所示:这就是我们前面提到的网格化方法,即对于任意一个点p,如果它所在的方格是D,那么必须定位到任意一个到p的距离小于等于r的点在D的中心,在一个由5×5的正方形组成的正方形区域中。在检查距离时,我们只需要计算这个子区域。我们用一个一维数组samples和一个N×N的二维数组grid来记录得到的采样点:samples保存了当前所有采样点的坐标,它的每个元素都是一个二维坐标(x,y)。grid[i,j]为整数,存储数组samples中第(i,j)个格子中采样点的下标。grid[i,j]=-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。最小采样间隔为每个小方块对角线的长度,即半径=sqrt(2)*dx。我们将目标采样点数设置为desired_examples=100000,这是一个视觉值,因为400x400的格子包含了160000个小方格,考虑到每个方格最多只有一个点,我们能得到的距离满足最大个数ofconstrainedpointsmustbelessthan160000.最初网格中没有点,因此需要将网格中的值设置为-1:grid.fill(-1)如何生成新点添加新随机时点,我们总是从现有点的附近开始随机选择一个位置,然后与已知点比较是否满足最小距离约束,如果是,则将其添加到现有点,否则丢弃并重新选择点.这里需要注意的是:当一个已经存在的点已经被填满后,我们后面添加一个新的点时就不用考虑它的邻近性了,我们可以用一个下标头来记录这个点。我们一致认为samples数组中靠近下标
=head的点。最开始head=0。samples是一个数组,长度为100K,并不代表我们真的可以得到这么多点,但是具体的个数无法提前确定,所以我们还需要用一个下标tail来记录目前为止获得的点数。最初tail=1因为我们将选择区域的中心作为第一个点。当然,这个初始点的位置可以是任意的。如前所述,当我们检查点p是否满足与现有点的最小距离约束时,不需要遍历检查所有点。查一下以p所在的正方形为中心的5×5正方形的正方形面积就可以了。我们编写一个单独的函数来检查一个点是否与现有点冲突:@ti.funcdefcheck_collision(p,index):x,y=indexcollision=Falseforiinrange(max(0,x-2),min(grid_n,x+3)):forjinrange(max(0,y-2),min(grid_n,y+3)):如果grid[i,j]!=-1:q=samples[grid[i,j]]if(q-p).norm()