Poisson Disk Sampling

Poisson Disk Sampling

2024年、9人の数学者からなるチームが、幾何学的ラングランズ予想に関する、実に1,000ページ近い証明を発表しました。これは純粋数学における偉大な成果であり、私は彼らが証明した命題のステートメントすら理解できないだろうと受け入れています。ましてや、その証明自体はなおさらです。

それとは全く対照的なのが、2007年にRobert Bridsonが発表した、1ページに満たない論文です。この論文は1,000件近い引用があり、完全に理解するのに10分もかかりません。コンピュータグラフィックスやシミュレーションで頻繁に発生する問題、すなわち「ランダムに配置するが、互いに近すぎないようにする」という問題に対するシンプルな解決策を提示しています。

例えば、森をプロシージャルに生成しようとして、木の配置方法が必要だとします。単純なランダムサンプリングの問題点は明らかです。

［画像: 木が重なっている例］

木が重なり合ってしまうのです。必要なのは、任意の2本の木の間に最小距離を設定できる能力です。この規則に従う木の分布は、「Poisson disk distribution（ポアソンディスク分布）」と呼ばれます。ナイーブな棄却サンプリング（rejection sampling）アプローチを試すことができます。これは、ランダムにダーツを投げ、他の点との最小距離内に入った点はすべて棄却するという方法ですが、より巧妙なデータ構造なしでは、各サンプルに対して衝突検出に線形時間がかかり、棄却率が急速に1に近づきます。Bridsonのアルゴリズムは、これを効率的に行う方法を提供します。

Bridson’s Algorithm

点の間の望ましい最小距離を r、d次元空間で作業していると仮定します。Bridsonのアルゴリズムは以下の通りです。

1. 空間を、辺の長さが r/√d のグリッドに分割します。これにより、各グリッドセルには最大で1つの点しか含められないことが保証されます。
2. active というリストを、空間から一様に選ばれたランダムな1つの点で初期化します。
3. active が空でない間、以下を繰り返します。
   1. active から要素 p を一様にランダムに選択します。
   2. p を中心とする、内半径 r、外半径 2r のアニュラス（環状領域）から、最大 k 回まで一様にサンプリングします。有効なPoisson diskサンプルが見つかった場合、グリッドを使用して効率的な衝突検出を行い、それを active に追加して、新しい p を選択します。k 回の試行で有効な点が見つからなかった場合、active から p を削除します。Bridsonは k=30 を推奨しています。

アニュラスを一様にサンプリングする最も簡単な方法は、単位ベクトル v⃗ ∈ℝᵈ と、区間 [1/2ᵈ, 1) から一様に選ばれた数 x を生成し、最終的なサンプルを 2rx¹/ᵈ ⋅ v⃗ とすることです。数年前に、これがなぜ機能するのかを説明する動画を作成しました。2次元では、単位ベクトルを選択することは、角度 θ ∈ [0, 2π) を選択することと同等です。高次元では、各成分が正規分布からサンプリングされたベクトルを正規化できます。

Improvements

Bridsonのアルゴリズムには、同じ数の点を生成するために必要なイテレーション数を劇的に削減する、2つの簡単な改良があります。1つ目は2次元で機能しますが、2つ目は高次元でも機能します。

まず2次元の改良から始めましょう。アルゴリズムが点 p を配置し、そのアニュラスをサンプリングして新しい点 q を取得する場合を考えます。p を q の親と呼びます。これらの点間の関係には貴重な情報が格納されています。必然的に q を中心とするアニュラスをサンプリングする際に、p に近すぎる点は考慮する必要のない角度の範囲全体が存在します。この範囲は、以下の図の破線で表されます。

［画像: 親子の関係による角度制限］

|p - q| = 1.50r

視覚的な直感は把握しやすいですが、それを数式に落とし込むのは面倒な三角法の演習です。詳細については省略し、証明なしで主張しますが、破線によって形成される円錐は角度 α を中心とし、その幅は 2β です。

\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*}

この公式で興味深いのは、β の計算式に現れる最小値です。これは、p と q の距離に応じて、アニュラスの内円または外円が円錐を制限する可能性があることを考慮しています。円錐とアニュラスの交差部分が完全に円内に含まれることを保証するために、最小値を選択する必要があります。p と q の間の距離が √3 ⋅ r を超えると、円錐の境界点が外円から内円にジャンプすることに注意してください。最小値の最初の項は外円との交差角度、2番目の項は内円との交差角度です。

これを実装するには、各点の親を格納するだけで十分です。その後、円錐の角度を計算し、次のサンプルの角度 θ を円錐の外の範囲で生成できます。

以下のグラフは、親の最適化ありとなしでBridsonのアルゴリズムによって生成された点の数を示しています。各データポイントは、辺の長さ ℓ=100、r=1 のグリッドで100回の試行の平均です。

［画像: 親の最適化による改善］

この改良は高次元にも一般化できる可能性がありますが、アニュラスの接触ベクトルを格納するためにより多くのスペースが必要になり、高次元ではアニュラスと球の交差の体積が比例して無視できるようになるため、効果は限定的になるでしょう。同様に、親に加えて点の子供を格納してアニュラスのセクションをさらに排除することを妨げるものはありませんが、これもより多くのストレージが必要になり、θ の選択がはるかに遅くなります。

次のサンプルの角度の選択方法を変更する代わりに、2番目の改良は次のサンプルの距離の選択方法を変更します。アニュラス内の各点から中心までの距離の分布を考えます。その累積分布関数（CDF）は、区間 [r, 2r] で xᵈ に比例します。その説明については、再び私の動画を参照してください。しかし、指数を d 以外の定数 c に変更するとどうなるでしょうか？そうすれば、点を中心に近づけたり遠ざけたりできます。c ≠ 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 が一様分布を与えることを思い出してください。

［画像: 距離分布のCDFとサンプル］

c = 2.00

スライダーで F₀(x) がゼロ除算で未定義になるにもかかわらず c=0 を設定できることに気づいたかもしれません。これを修正するために、F₀(x) = lim_{c→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 ⋅ 2ˣ であるべきです。ここで x は区間 [0, 1) の一様乱数です。それ以外の場合は、2ry¹/ᶜ を使用します。ここで y は 1 と 1/2ᶜ の間の均一乱数です。区間の境界は、c が正か負かに応じて入れ替わります。

以下のヒートマップは、c がBridsonのアルゴリズムで生成される点の数に与える影響を示しています。以前と同じ実験設定を使用し、辺の長さ ℓ=100、r=1 のグリッドで100回の試行を行いました。

［画像: cの値による改善］

これは、c を非常に負の値に設定するか、c が負の無限大に近づく極限を取るべきであることを示唆しているようです。これにより、各点が親から正確に r の距離を持つようになります。これは、生成される点の数を最大化し、より密な配置を作成しますが、分布が「ランダムに」感じられるという犠牲を払います。c = -∞ の極端なケースでは、点の長い文字列や、制限された距離が届かないギャップなどの多くのアーティファクトが発生します。後から点の生成方法のツリーを再構築することも難しくありません。

したがって、点の密度を最大化することと、ランダム性を維持することのバランスを取る必要があります。もし 15 ≤ k ≤ 40 のいずれかの値を固定するなら、親の最適化を使用する場合、c = -1.4 - 17/√k と設定するのが最良であると私は見つけました。この式は、生成される点の数が、一様かつ最大ポアソンディスクサンプラーからの期待される出力とほぼ一致するように経験的に導出されました（これについては後述）。15から40までの各 k について、半径 r/2 の円を各点の周りに描いたときに、それらの円が全体の面積の54.7%を占めるようにする c の値を見つけるために二分探索を行いました。このパーセンテージは、ランダム逐次吸着モデルにおける円盤の飽和カバレッジです。もちろん、これは2次元にのみ適用され、c は高次元では異なるように調整する必要があります。

Stippling

これまでは r を一定に保ってきましたが、その必要はありません。点の間の最小距離を関数 r : ℝᵈ → ℝ として動的に設定できます。したがって、点 p を配置した後、内半径 r(p) のアニュラスをサンプリングします。これの楽しい応用は、画像内の各ピクセルの明るさとして r を定義して、スタンプリング効果を生成することです。左の図は白黒で行われ、右の図は3つのPoisson diskサンプルのセットを組み合わせて、各カラーチャンネルに使用しています。

［画像: スタンプリングの例］

Birdsonのアルゴリズムは本質的に逐次的なものですが、並列実行するように設計された他のアルゴリズムもあり、大幅なパフォーマンス向上をもたらします。私が最も気に入っているのはPixelPieで、これは完全にGPU上で実行されます。このアルゴリズムを使用して、Poisson disk samplingを使用したリアルタイムビデオスタンプラーであるPoisson Camを作成しました。これは、シェーダープログラミング、Rust、ストリーム圧縮アルゴリズム、そしてもちろんPixelPieアルゴリズム自体を学んだため、私が最も楽しかったプロジェクトの1つでした。

Maximality and Uniformity

2022年、Scott A. Mitchellは、2次元でPoisson diskサンプルを生成するための美しく新しい方法を紹介する、非常にクールな論文を発表しました。この論文は、私の知る限り、発表以来ほとんど注目されていません。この最後のセクションは、それに捧げられています。

アルゴリズムを説明する前に、それがBridsonのアルゴリズムよりも優れている3つの点を説明したいと思います。

1. Maximality（最大性）：アルゴリズムが終了した後、Poisson diskプロパティに違反することなく、これ以上点を配置できないことが保証されます。
2. Uniformity（一様性）：アルゴリズムは、すべての最大Poisson diskサンプルのセット全体にわたる一様分布をサンプリングします。
3. Determinism（決定性）：アルゴリズムは棄却サンプリングに依存しないため、点の配置に失敗することはありません。

最大性と一様性は、階層的ダーツ投げなど、多くのアルゴリズムで達成されていますが、Mitchellのアルゴリズムは、何らかの棄却サンプリングなしでこれを達成した最初のアルゴリズムです。パフォーマンスも非常に優れています。厳密なベンチマークは行っていませんが、Mitchellの実装は、私が実装したBridsonのアルゴリズム（親の最適化を使用）と同程度に高速に実行され、ほぼ同じ数の点を生成することがわかりました（k=20、c=-5.2 の場合）。私の知る限り、このアルゴリズムの唯一の欠点は、大幅に複雑であることです。以下は、大幅に簡略化された概要です。

1. 空間を、辺の長さが r/√d のグリッドに分割します。
2. もう1点配置する余地がある間、以下を繰り返します。
   1. グリッドから、残りの面積で重み付けされたセル c をランダムに選択します。
   2. c を、互いに素な三角形とチョーク（chock）に分解します。1
   3. 面積で重み付けされたランダムな三角形またはチョーク t を選択します。
   4. t から一様に点 p をサンプリングします。それを最終セットに追加します。
   5. グリッドから p を中心とする半径 r の円を切り出します。

元の論文では、アルゴリズムの詳細な説明が非常によく行われているため、私ができる最も役立つ貢献は、次の視覚化です。

［画像: Mitchellのアルゴリズムの視覚化］

Footnotes

1. チョーク（chock）とは、円、放射状の線、および接線によって境界付けられる3辺の形状です。↩