關于 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