关于 Fast Poisson-Disk Sample Generation
Chaos Chiao
一、随便说说
从MSRA回来后仔细地想了很多事情，看着自己一度觉得很了不起的渲染器，觉得有点可笑。大而全的东西，注定是非常容易失败的，尤其对于单干，把东西铺得太开，就很危险了。尽管我再做GS的时候已经有过这样的教训了。
感觉这段时间有非常强烈变强的欲望，但知道欲速则不达，所以还是停下了项目，重新找出来了很多布满尘土的数学书，翻着做题。不知为什么，很多以前一看见就挠头的题目，现在却变得非常简单，可能是长期的编程彻底地改变了我的思维方式吧。
但还是受不了Coding的诱惑，就像抽烟一样（虽然我不抽烟，但也知道它很上瘾），于是想找一些课题研究一下，便打开SIGGRAPH 2006的Paper list，搜索一些感兴趣的。第一眼就看到了< A Spatial Structure for Fast Poisson-Disk Sample Generation >这篇文章。感觉挺惊讶的，近来很少看到Sampling Pattern的研究，似乎已经没什么可以再深入下去的余地了，但它却提到了在O(NlogN)的复杂度内可以生成Poisson Disk Pattern。不可思议呢，以前我在做Sampling & Reconstruction的时候曾经尝试过研究Poisson Disk，但还是用一贯的思路，就是所谓Dart-Throwing。也有想要改进这种方法的念头，但最后还是没有深入下去，因为用Sobol Pattern生成的图像质量已经非常不错了，而且Sobol可以预计算序列，运行时几乎没有任何耗时。
稍微看一下，了解了大概的算法，觉得虽然有一定的难度，但要做一个这样的Implementation还是不太困难的。奋战一天下来，算是做出来了一个，然后下面分享一些我的想法。
二、Sampling & Reconstruction、Poisson Disk
这已经是我第二篇文章讨论Sampling & Reconstruction了，上一篇只是因为ChinaDV有一个疯子出来说什么反混淆与反走样，我详细解说一下什么是Sampling & Reconstruction。当然，这其实属于图形学、乃至信号处理里面最简单的概念。
Sampling主要分为Uniform Sampling、Adaptive Sampling和Stochastic Sampling。Uniform Sampling是直接在一个网格或者经过抖动的网格Pattern上进行采样，它受网格分辨率限制。当网格分辨率远低于Source的频率时，Uniform Sampling会出现非常大的误差，其直接结果就是不连续的锯齿，未经抖动的网格采样还可能会产生严重的水纹波。Adaptive Sampling是在Uniform Sampling的基础上进行改良的，当邻近两个采样点之间的差大于阀值时，再两点之间插入一个新的采样点，这样可以极大地避免不连续情况。但当首个Sample Pattern的分布和Source的频率相距太大时，虽然邻近两点的差值在阀值内， 有可能Source在两点之间有更多的连续变化存在，这样的重构会忽略掉采样点之间的高频变化。
Stochastic Sampling是完全的随机采样，它不受网格限制，使用一系列在采样区域内的散乱采样点进行采样，不根据网格划分采样点。Stochastic Sampling非常依赖Pattern的特性，Pattern的分布直接决定了采样的质量。Cook提出了使用Poisson Disk分布的Pattern会非常适合二维的采样，他指出人眼中的感光细胞也成Poisson Disk分布。所谓Poisson Disk Pattern即所有的采样点到其余所有的采样点的距离都大于一个阀值，可以认为未抖动的网格是Poisson Disk Pattern的其中一个特例。
这样Poisson Disk就要求任意两点之间的距离不小于阀值，比如10x10的区域内生成100个(以上)的采样点，阀值可以采用0.9（大于等于1将有可能使有的采样点不能放到区域内）。
传统生成Poisson Disk序列的方法为Dart-Throwing。可以把这个原始算法看作买六合彩。它不断“Throw”一个随机的采样点，然后和已有的采样点集合比较距离，若遇到小于阀值的就Discard，再重新“Throw”一个新的随机采样点；如果符合条件则“Dart”中了，添加到采样点集合里。就这样不断循环直到完全填满区域，或者生成的采样点“足够多”为止。如果采样区域非常大、或者采样点数目巨大，那么要计算完所有采样点的几率真的比中六合彩还要低得多。上了10,000个点的填充计算是基本上不可能完成的（根据我的经验……），所以在程序运行时生成Poisson Disk Pattern是非常不切实际的做法。
可是Poisson Disk分布的确太好，太适合各种图像重构的采样了，所以很多人会预计算一个足够大的Pattern，再把它Tile到采样区域里。
另外也提一下Sobol序列，Sobol是一个Quasi-Monte Carol序列，可以扩展到任意维度，它是一个稳定的、覆盖率非常好的随机序列，用在Stochastic Sampling中也可以得到非常高质量的采样效果。Sobol已经被广泛运用到各种图形图像系统中。
三、改进的Dart-Throwing
Dart-Throwing属于一种随机算法，无法计算它的大O，因为它甚至不可能结束，所以一直无法在程序中用作即时计算。但Poisson Disk的分布有非常独特的特点，可以让我们预先一些不可能中标的区域Discard掉，然后在100%中Dart的区域中Throw，这样可以保证我们的算法每次都是神投手。
哪些区域会100%中标？先看Poisson Disk的定义，即距离其它所有采样点的距离都在阀值外。这样在阀值距离内的所有范围都不可能中。如果阀值为r，已知采样点P，则以P为中心半径r~2r的圆环内是最佳的采样点出现区域。把随机数映射到这个范围内，可以得到符合条件的采样点。根据Paper，可以建立一个能进行布尔运算的贝壳形representation，每次把新的采样点和圆环（或者贝壳）相减，得到一个贝壳状的区域，那么新的采样点可以在贝壳区域内生成。
可是贝壳状的区域，而且涉及二维布尔运算，这会非常复杂，况且r~2r范围内生成的采样点，并不能保证覆盖率达到最大化，也有可能出现采样点过于稀疏而导致无法插入新点。
四、Boundary Sampling
当我们把r~2r的定义再进一步深入研究下去，可以得到实际上如果新建立的点在半径等于r的圆上，可以使覆盖率最大化，即浪费的空间最少（但这部等于所有采样点之间的距离都是r）。那么我们就把随机序列映射到一个贝壳区域内简化到直接把单个随机数映射到圆（或者弧线）上。
套用原来的思路，建立一种以圆弧为基础的representation，可以进行圆弧与圆之间的布尔运算，相比之下计算要简单很多，只用极坐标和简单的三角函数计算已经可以完成了。然后通过指定一个原始的采样点，我们可以一步一步地在新的点的圆弧上生成符合条件的新点。
在正方形区域内生成一个Poisson Disk Pattern，从点0开始。可以看到以点0为中心，r为半径的圆上任意一点插入的采样点都符合Poisson Disk分布，所以可以很简单地直接在圆上进行随机插值。
点1是新插入的采样点。分别以两点为圆心、r为半径的圆面相交的区域是不可能出现符合Poisson Disk分布的区域，所以可以直接把这部分的圆弧删除掉（蓝色），剩下（白色）的圆弧则是可以插入新的采样点的区域。
如此这般，继续在可能的弧线上插入新的采样点，这些点集就完全符合Poisson Disk分布，而且可以达到最优覆盖率。
当所有的圆弧都用完了以后，我们就得到了一个非常好的Poisson Disk Pattern了。
五、Implementation
关键在于可以进行布尔运算的圆弧的representation，我用了极坐标来表示圆弧，这样二维布尔运算变成了简单的加减运算。同时我建立了一张表来快速查询邻近的采样点，而空闲采样点集合则使用了std::set储存。其实上面已经解释得挺清楚的了，附上动态生成Pattern的Demo和代码。Demo打开后不断按 Space 就可以看到采样点的生成过程。
Demo & Source
