古典的な論文を実装するための60行のコード：ポアソンディスクサンプリングを完了するのに0.7秒、Numpyよりも100倍高速

この記事はAI新メディアQuantum Bit（公開アカウントID：QbitAI）より許可を得て転載しています。転載の際は出典元にご連絡ください。

ランダムで均一な点で構成されるパターンは、植物や動物ではすでに一般的です。

ベイベリー、イチゴ、ライチ、ランブータンなどの果物の表面には粒状または毛状の構造があり、果物の表面にランダムかつ均一に分布しています。

同様の模様は、誰もが好んで食べる胃袋などの動物にも見られます。

同様に、コンピューターシミュレーションでは、空間内でポイントをランダムかつ均一に生成する必要があるシナリオが多数あります。

動物の毛を生成するときの毛穴の位置、マルチプレイヤーゲームでのプレイヤーの出生位置、森林を生成するときの木の位置など。

これらのシナリオに共通する特徴は、いずれも2 点間の距離が下限値以上である必要があることです(この下限値は事前に設定されており、これを変更することで生成される点間の間隔を制御できます) 。

ただし、完全にランダムに生成されたポイントを直接使用すると、一部の場所が「密集」し、一部の場所がまばらになるなど、非常に不均一な分布結果が得られる可能性が高くなります。

このようなポイントを使用して髪の位置生成をシミュレートすると、効果は非常に悪くなります。

したがって、ポイントを生成するプロセスに距離判定を追加して、要件を満たさないポイントを排除する必要があります。

以前は、 NumPyを使用してこのような効果を生成するには70 秒ほどかかることが多く、非常に非経済的でした。

現在、Taichi Graphics はTaichiをベースにした超高速アルゴリズムを実装しています。同じエフェクトが単一のCPUスレッドで実行され、そのようなパターンを生成するのにわずか0.7 秒しかかかりません。これは約100 倍の速さです。

彼らがどうやってそれをやったか見てみましょう。

Bridsonアルゴリズムを使用して実装

以前は、ダーツを投げるアルゴリズムが一般的でした（目隠しをした人がランダムにダーツを投げるようなもの） 。

毎回、エリア内のポイントがランダムに選択され、そのポイントとすでに取得されているすべてのポイントとの間の「競合」がチェックされます。

ポイントと既に取得されているポイント間の最小距離が指定された下限より小さい場合、ポイントは破棄されます。それ以外の場合は、適格なポイントとして既存のポイントのセットに追加されます。

十分なポイントが得られるまで、または N 回連続して失敗が発生するまで(N は正の整数) 、この操作を繰り返します。

しかし、このアルゴリズムは非常に非効率的です。

取得するポイントの数が増えるにつれて、競合の可能性がますます高くなり、新しいポイントを取得するのに必要な時間がますます長くなり、現在のポイントと既存のすべてのポイント間の距離を比較するたびに効率も低下します。

対照的に、 Bridson アルゴリズムはより効率的です。

このアルゴリズムの原理は、ロバート・ブリッドソンの2007年の論文「任意次元での高速ポアソンディスクサンプリング」[1]に由来しています（論文は非常に短く、A4用紙1ページのみです） 。タイトルと序文を除くと、実際のアルゴリズムの内容は短い段落だけです。

冒頭のアニメーションは、 400 x 400グリッド領域で実行される Bridson ディスク サンプリング アルゴリズムを示しています。アルゴリズムは、均等に分散された 100K 個のポイントを取得しようとしますが、実際には約 53.7K 個のポイントが生成されます。

このアニメーションは、単一の CPU スレッドで実行される Taichi を使用して生成され、コンパイル時間の計算を除いて 0.7 秒強かかります。同じコードをNumPyに変換すると約 70 秒かかります。 [2]

上のアニメーションからわかるように、ブリドソン アルゴリズムはキャベツの成長プロセスと非常によく似ています。つまり、シード ポイントから始めて、層ごとに新しいポイントを外側に追加していきます。

新しいポイントを追加するたびに、そのポイントは最も外側のポイントの周囲に配置され、最も外側のレイヤーを可能な限りカバーします。

存在するすべてのポイントまでの距離を毎回確認するのを避けるために、Taichi はグリッディングと呼ばれる手法を使用します。

空間全体がグリッドに分割されています。チェックが必要なポイントについては、そのポイントが配置されているグリッドを見つけ、そのポイントと隣接するグリッド内のポイント間の最小距離をチェックするだけです。

この距離が指定された下限よりも大きい限り、それ以上のポイントをチェックする必要はありません。この手法はグラフィックスや物理シミュレーションで非常によく使用されます。

このサンプリング プロセスは、スレッドが新しいポイントに「侵入」すると、他のすべてのスレッドの距離の判断が変わるため、並列化が困難です。したがって、Taichi はこのアルゴリズムを実行するために単一の CPU スレッドのみを使用します。

 ti . init ( アーキテクチャ= ti . cpu )

上記のコードでは、プログラムを CPU 上で実行するために arch=ti.cpu が指定されています。

シングルスレッド + CPU なので、純粋なPythonで書けばよいのではないかと思われるかもしれません。心配しないでください。計算部分は ti.kernel 関数に配置されます。この関数は Python 仮想マシンでは実行されませんが、Taichi によってコンパイルおよび実行されるため、純粋な Python 実装よりも何倍も高速になります。

Bridson アルゴリズムの具体的な実装を紹介する前に、この Taichi プログラムに何行のコードが含まれているかを推測してみませんか?

Taichiのインストールとインポート

まず、より多くの機能とさまざまなプラットフォームでのより安定したサポートを備えた最新の Taichi リリース バージョンを使用することをお勧めします。この記事の執筆時点での最新バージョンは1.0.3です。

 pip インストールtaichi == 1.0.3

次に、コードの先頭に次のように記述します。

 taichi をti としてインポートする
taichi.math をtm としてインポートします

これにより、Taichi と Taichi の数学モジュールがインポートされます。一般的に使用される数学関数に加えて、math モジュールは非常に便利なベクトル演算も提供します。

準備

ポアソン サンプリング アルゴリズムでは、サンプリング ポイント間の距離には下限 r があります。

全体の領域はN×N個の同じサイズの正方形で構成されたグリッド領域であると仮定します。この場合、各小さな正方形の対角線の長さはr 、つまりグリッドの辺の長さはr/√2になります。

したがって、任意の小さな正方形には最大で 1 つの点が含まれます。次の図に示すように:

これは、先ほど説明したグリッド化方法です。つまり、任意の点 p について、その点が位置する正方形を D とすると、p からの距離が r 以下である任意の点は、D を中心とし、5×5 の正方形で構成される正方形領域内に位置している必要があります。

距離を確認するときは、このサブ領域についてのみ計算する必要があります。

得られたサンプリングポイントを記録するために、 1次元配列サンプルとN×Nの2次元配列グリッドを使用します。

1. samples は現在サンプリングされたすべてのポイントの座標を格納し、各要素は 2 次元座標 (x, y) です。
2. grid[i, j]は、配列サンプル内の(i, j)番目のグリッド内のサンプリングポイントのインデックスを格納する整数です。 grid[i, j] = -1は、このグリッドにサンプリングポイントがないことを意味します。

したがって、初期設定は次のように記述できます。

 グリッド_n = 400
 res = ( グリッドn 、 グリッドn )
 dx = 1 / 解像度[ 0 ]
 inv_dx = 解像度[ 0 ]
 半径= dx * ti . sqrt ( 2 )
 希望サンプル数= 100000
 グリッド= ti . フィールド( dtype = int 、 形状= res )
 サンプル= ti.Vector.field ( 2 , float , shape = desired_samples )

ここで、グリッド サイズは 400x400 に設定され、グリッドが占める平面領域は [0,1]×[0,1] なので、グリッドのステップ サイズは dx = 1/400 になります。最小サンプリング間隔は各小さな正方形の対角線の長さ、つまり半径 = sqrt(2)*dx です。

サンプリング ポイントの目標数を desired_examples=100000 に設定しましたが、これは 400 x 400 のグリッドに 160000 個の小さな正方形が含まれているため推定値です。各正方形には最大で 1 つのポイントがあることを考慮すると、距離制約を満たすポイントの最大数は間違いなく 160000 未満です。

最初はグリッド内にポイントがないので、グリッド内のすべての値を-1に設定する必要があります。

 グリッド. 塗りつぶし( - 1 )

新しいポイントを生成する方法

新しいランダム ポイントを追加するときは、常に既存のポイントの近くの位置をランダムに選択し、それを既知のポイントと比較して、最小距離制約を満たしているかどうかを確認します。 満たされている場合は、既存のポイントに追加し、満たされていない場合は破棄して新しいポイントを選択します。

以下の点に留意することが重要です。

1. 既存のポイントの周囲の領域がすでに塗りつぶされている場合、後で新しいポイントを追加するときにその近傍を考慮する必要はありません。これを記録するために下付き文字のヘッドを使用できます。サンプル配列内の添え字 < ヘッドのポイントはすでに入力されているため、これ以上考慮する必要がないことに同意します。考慮する必要があるのは、添え字 >= ヘッドのポイントのみです。最初はヘッド = 0 です。
2. samples は長さ 100K の配列です。これは、実際にこれほど多くのポイントを取得できるという意味ではありませんが、具体的な数を事前に決定することはできないため、現在取得されているポイントの数を記録するために下付き文字の末尾も使用する必要があります。最初のポイントとして領域の中心を選択するため、最初は tail = 1 になります。もちろん、この初期点の位置は任意です。
3. 前述したように、点 p が既存の点との最小距離制約を満たしているかどうかを確認する場合、すべての点を走査して確認する必要はありません。 p が位置する正方形を中心とした5×5 の正方形で構成される正方形領域を確認するだけです。

ポイントが既存のポイントと競合するかどうかを確認するロジックを別の関数に記述します。

 @ ti.func
 def check_collision ( p 、 インデックス):
 x , y = インデックス
衝突= False
 iが範囲内( max ( 0 , x - 2 ), min ( grid_n , x + 3 ) ) の場合:
 jが範囲内( max ( 0 , y - 2 ), min ( grid_n , y + 3 ) ) の場合:
 グリッド[ i , j ] != - 1の場合:
 q = サンプル[ グリッド[ i , j ]]
 ( q - p ) .norm () < 半径- 1e-6 の場合:
 衝突= 真
リターン衝突

ここで、p はチェックする点の座標であり、index=(x, y) は p が配置されている正方形の下付き文字です。

x-2 <= i <= x+2 かつ y-2 <= j <= y+2 を満たすすべての添字 (i, j) を反復処理し、正方形 (i, j) 内にすでに点があるかどうか、つまり grid[i, j] が -1 に等しいかどうかを確認します。そうであれば、p までの距離が半径より小さいかどうかを確認し、対応する判定を返します。

準備が完了したら、正式なサイクルを開始できます。

 @ ti.kernel
 def poisson_disk_sample ( 望ましいサンプル数: int ) - > int :
 サンプル[ 0 ] = tm.vec2 ( 0.5 )
 グリッド[ int ( grid_n * 0.5 ), int ( grid_n * 0.5 )] = 0
 頭、 尾= 0、1
 head < tail かつhead < desired_samples の場合:
 source_x = サンプル[ ヘッド]
 頭+= 1

 _が範囲( 100 ) 内にある場合:
 シータ= ti . ランダム( ) * 2 * tm.pi
 オフセット= tm.vec2 ( tm.cos ( theta ) , tm.sin ( theta ) ) * ( 1 + ti.random () ) * 半径
new_x = source_x + オフセット
新しいインデックス= int ( 新しい x * inv_dx )

 0 <= new_x [ 0 ] < 1 かつ0 <= new_x [ 1 ] < 1 の場合:
 衝突= check_collision ( 新しい x 、 新しいインデックス)
 衝突せず、 末尾がdesired_samples 未満の場合:
 サンプル[ 末尾] = new_x
 グリッド[ new_index ] = 末尾
尾+= 1
 戻る尾

まず、エリアの中心、つまり座標 (0.5, 0.5) の点を初期点として選択し、それを「シード」として使用して、ランダムな点をエリア全体に徐々に広げていきます。

次の while ループはアルゴリズムの本体です。このループはシリアルに実行され、1 つのスレッドのみを占有します。

考慮する必要がある最初のポイントsamples[head]を見つけるたびに、それを中心とする半径[radius, 2*raidus]の円内の100個のポイントをランダムに選択し、これらの100個のポイントが[0,1]×[0,1]の範囲を超えていないかどうか、および既存のポイントと競合していないかどうかを確認します。

すべての条件が満たされている場合、それは適格なポイントです。その座標とグリッドの添え字をサンプルとグリッドに更新し、既存のポイントの末尾の数を 1 増やします。 100 ポイントすべてがチェックされた後、既存のポイントのセットにさらにポイントを追加できます。

半径 [radius, 2*raidus] の円内でサンプリングすると、最小距離制約を満たしながら、既存のポイントからそれほど離れていないポイントを取得できることに注意してください。

100 個のポイントをすべてチェックした後、samples[head] の周囲に新しいポイントを配置するための空白領域がないと想定できるため、head を 1 増やして、次の samples[head] の近くの領域を再チェックします。

ループは、すべてのポイントの周囲のスペースが埋められたとき、つまり head = tail になったとき、または desired_samples ポイントが取得されたとき、つまり tail = desired_samples になったときに終了します。このとき、サンプル内の添字が 0 から tail-1 までの範囲にある点はすべて存在する点です。

アニメーション効果を表示

わずか数行のコードで、サンプリング プロセス全体をアニメーション化できます。

 num_samples = ポアソンディスクサンプル( 必要なサンプル)
 gui = ti . GUI ( "ポアソンディスクサンプリング" 、 res = 800 、 background_color = 0xFFFFFF )
 カウント= 0
 速度= 300
 guiの実行中:
 gui.circles ( samples.to_numpy ()[: min ( count * speed , num_samples )] ,
 色= 0x000000 、
 半径= 1.5 )
 カウント+= 1
 GUI . 表示()

ここでは、生成された 300 ポイントごとに 1 フレームを描画するようにアニメーションの速度を制御します。

プログラムの主要なポイントをすべて説明したので、各部分をまとめてみましょう。

 taichi をti としてインポートする
taichi.math をtm としてインポートします
ti . init ( アーキテクチャ= ti . cpu )

 グリッド_n = 400
 res = ( グリッドn 、 グリッドn )
 dx = 1 / 解像度[ 0 ]
 inv_dx = 解像度[ 0 ]
 半径= dx * ti . sqrt ( 2 )
 希望サンプル数= 100000
 グリッド= ti . フィールド( dtype = int 、 形状= res )
 サンプル= ti.Vector.field ( 2 , float , shape = desired_samples )

 グリッド. 塗りつぶし( - 1 )

 @ ti.func
 def check_collision ( p 、 インデックス):
 x , y = インデックス
衝突= False
 iが範囲内( max ( 0 , x - 2 ), min ( grid_n , x + 3 ) ) の場合:
 jが範囲内( max ( 0 , y - 2 ), min ( grid_n , y + 3 ) ) の場合:
 グリッド[ i , j ] != - 1の場合:
 q = サンプル[ グリッド[ i , j ]]
 ( q - p ) .norm () < 半径- 1e-6 の場合:
 衝突= 真
リターン衝突

@ ti.kernel
 def poisson_disk_sample ( 望ましいサンプル数: int ) - > int :
 サンプル[ 0 ] = tm.vec2 ( 0.5 )
 グリッド[ int ( grid_n * 0.5 ), int ( grid_n * 0.5 )] = 0
 頭、 尾= 0、1
 head < tail かつhead < desired_samples の場合:
 source_x = サンプル[ ヘッド]
 頭+= 1

 _が範囲( 100 ) 内にある場合:
 シータ= ti . ランダム() * 2 * tm . pi
 オフセット= tm.vec2 ( tm.cos ( theta ) , tm.sin ( theta ) ) * ( 1 + ti.random () ) * 半径
new_x = source_x + オフセット
新しいインデックス= int ( 新しい x * inv_dx )

 0 <= new_x [ 0 ] < 1 かつ0 <= new_x [ 1 ] < 1 の場合:
 衝突= check_collision ( 新しい x 、 新しいインデックス)
 衝突せず、 末尾がdesired_samples 未満の場合:
 サンプル[ 末尾] = new_x
 グリッド[ new_index ] = 末尾
尾+= 1
 戻る尾

num_samples = ポアソンディスクサンプル( 必要なサンプル)
 gui = ti . GUI ( "ポアソンディスクサンプリング" 、 res = 800 、 background_color = 0xFFFFFF )
 カウント= 0
 速度= 300
 guiの実行中:
 gui.circles ( samples.to_numpy ()[: min ( count * speed , num_samples )] ,
 色= 0x000000 、
 半径= 1.5 )
 カウント+= 1
 GUI . 表示()

コードの総行数: 60 。

もう一つ

具体的には、このコードは次の 2 つの操作を実装します。

1. 完全なポアソンサンプリングアニメーションは60 行のコードで実装されています。
2. 400 x 400グリッドで53,000 ポイントが収集されましたが、所要時間は1 秒未満でした。

関連するコードは、記事の最後にある元のリンクにあります。

厳密に言えば、この記事で実装されているアルゴリズムは、 Bridson の論文に記載されているアルゴリズムとは少し異なります(より単純です)が、効果は似ています。

違いがわかりますか？ (ヒント: 元の論文のステップ 2 の「アクティブ リスト」の扱いに注意してください)

プロジェクトアドレス:

https://github.com/taichi-dev/poisson-sampling-homework