はじめに

これはレイトレ合宿アドベント・カレンダー(2018)の記事です。ボロノイ図の生成の仕方、またそれを用いたブルーノイズの特性を持つ点列の生成方法、点描などの応用についてざっと紹介します。

下準備

はじめに基本的な、ボロノイ図、パワー図、ボロノイ図の生成法、ブルー・ノイズを紹介します。ここで紹介するもののほとんどは２次元に特化しており、離散化された場を対象としています。

ボロノイ図

Wikipediaには「ボロノイ図（ボロノイず、英語: Voronoi diagram）は、ある距離空間上の任意の位置に配置された複数個の点（母点）に対して、同一距離空間上の他の点がどの母点に近いかによって領域分けされた図のことである。」とあります。それぞれの領域をセルと呼びます。

CGではフレークの形状や、Worleyノイズなどでよく出てくるので知っている人がほとんどでしょうか。Inigo Quilez氏のサイトに行けば様々な応用例を見ることができます。http://www.iquilezles.org/www/articles/voronoilines/voronoilines.htm

自然界でもボロノイ図のようなパターンはいたるところで現れます。

Can you guess what this sea creature is? Some call it the lava lamp of the sea! (Footage via @EVNautilus) pic.twitter.com/vfaS4AjUCz
— Azula (@azula) 18 August 2018

パワー図

ボロノイ図では与えられた点$p$とセル$i$の重心$s_i$との距離を以下で測りますが、

l(\mathbf{p}, \mathbf{s_i})=\sqrt{(\mathbf{p}-\mathbf{s_i})^2}

距離の計算には任意の関数を使うことができます。特に以下を用いたものをパワー図といい、$r_i$による重みづけを可能としながら、境界が直線になるような領域の分割を行うことができます。

l^2(\mathbf{p}, \mathbf{s_i})-r_i^2

詳しく知りたい方には杉原厚吉氏の「なわばりの数理モデル」という書籍をお勧めします。

単純な生成方法

さて、ボロノイ図はどのように描くのでしょうか。定義に倣うと、それは最近傍点を見つけることに帰着します。最も単純な方法は全ての点を保持するkd-treeやBVHなどのデータ構造を利用し、単純に最近傍点を求めるというものでしょうか。うまく作れば、これでも十分なスピードが出せます。

Jump Flooding Algorithm

2次元に限定した場合、ボロノイ図の生成にはJump Flooding Algorithm(以下JFA)という効率の良いアルゴリズムを使うことができます。出来上がったボロノイ図は近似ではあるものの、その亜種を使用した場合、誤差は非常に小さく抑えられるため魅力的なものとなっています。

JFAという名前に憶する必要はありません。少し変わった方法で塗りつぶしを行うアルゴリズムです。通常の塗りつぶしでは、隣接するピクセルに値を入れていきますが、JFAでは、下の図にあるように、大きなステップで塗りつぶしを始めます。はじめのラウンドでは画像サイズの半分の幅を使い、ラウンドが進むにすれ、ステップサイズを半分にしていき、ステップサイズが1になれば終了です。

画像の右端を見ればわかりますが、toroidal domain(ドーナツのように上下左右がつながっている領域)でもうまくボロノイ図を計算することができます。toroidal domainで計算したものはテクスチャとして繰り返し貼ってもつなぎ目がありません。

アルゴリズムの説明ではscatterを用いましたが、scatterを用いると、並列化した際に、スレッドが競合するため、atomicな演算を使用しなければなりません。なので、並列化をするにはgatherを使用して実装します。ボロノイ・セルの数が少ない場合、シングルスレッドでscatterを利用した方が効率が良いですが、後半のラウンドにおいては、並列化できるためgatherを用いた方が効率が良くなります。CPUで実装する場合はハイブリッドがよいでしょう。

単純な実装では解像度が２のべき乗に限定されるという制限が付きます。FFTのようですね。

ブルー・ノイズ

ランダムに点を打った場合は、ある場所では近傍の点どうしがくっついていたり、ある場所では少し離れていたりします。大雑把にいってしまうと、ブルー・ノイズと呼ばれる特性を持つ点の集合は、ランダムに配置されているものの、周りとの点との距離が比較的均等であるといった特徴を持ち、その生成にはボロノイ図を利用するものが多々あります。

アルゴリズム

さて、ここからブルーノイズの性質を持つ点列の色々な生成法を紹介していきたいと思います。

Lloydの弛緩法

ブルー・ノイズと呼ばれる特性を持つ点列は様々な方法で生成できますが、代表的なものにLloydの弛緩法を使うものがあります。この弛緩法には先ほど説明したボロノイ図を使います。

アルゴリズムは非常に簡単で、任意の数の点を最初に生成し(重なっていなければどのような方法でも構いません)、それらに対してボロノイ図を作ります。そのあと、各点を自身が含まれるボロノイ・セルの重心に移動させます。移動した点に対してさらにボロノイ図を生成し、同様の処理を繰り返すと、次第に点は移動しなくなり、収束していきます。出来上がったボロノイ図は重心ボロノイ分割と呼ばれ、しばしばCVT(Centroidal Voronoi Tessellation)と省略されます。

この手法は、ある濃度の分布に従って点列を分布させるという用途には向きません。収束するまでに非常に長い時間がかかる、または不可能であるからです。(重心の計算の際に、陰に密度関数をぼかしてしまうという欠点があります。) 具体的に言うと、明るい場所には少しの点を、暗い場所にはたくさんの点を打つ、といった点描には使えません。

CCVT

Capacity Constrained Voronoi Tessellationの略称です。
各ボロノイセルが決められた容量（2次元ではここでは面積と捉えてもらえばイメージしやすい）を持つように制約をかけて最適化する方法です。Lloydの弛緩法の、生成された点列に規則性が見られるという欠点も克服でき、質のいいブルー・ノイズ特性を持つ点列が生成できます。以下は2次元で生成する場合の基本的なアルゴリズムです。

1. 画像内に必要な数の点を生成し、ボロノイ図を生成します。
2. 全てのピクセルにどのボロノイ・セルに属するかという情報を持たせます。(ランダムでも構いません。また容量の調整は属するピクセルの数を調整することで行えます。その場合ステップ1ではボロノイ図の生成をする必要はありません。)
3. 適当な2点を選び、2点の所属するセルを入れ替えて、「それぞれの点とそれが所属するセルの重心との距離」の和が減少するなら、所属するセルの入れ替えを行います。
4. 入れ替えた場合、各セルの重心は移動するので、Lloydの弛緩法と同様にしてセルの重心位置を更新します。
5. 3.に戻ります。

元の論文のアルゴリズムの動作は遅いです。ただ、いくつか高速化手法が提案されていて、代表的なものには次にあげる二つの文献があります。

Fast Capacity Constrained Voronoi Tessellation

いくつか最適化のテクニックが挙げられています。

処理の対象としているセルを$\mathbf{s_i}$と$\mathbf{s_j}$としましょう。セルに含まれる平均のピクセル数が$\frac{m}{n}$個であった場合、すべてのペアを試した時の計算量は$O(\frac{m}{n}log(\frac{m}{n}))$となります。このコストがCCVTの処理時間の大部分を占めているので、これを減らすことを考えます。

ここで$\mathbf{s_i}$はi番目のボロノイ・セルの重心位置、$\mathbf{p}$はいま見ているピクセルとすると、ボロノイ・セルをアサインしなおした場合のエネルギーの変化分$\Delta E$は次の式で表すことができます。

\Delta E=|\mathbf{p}-\mathbf{s_i}|^2-|\mathbf{p}-\mathbf{s_j}|^2

処理を簡単にするために、いったん全てのピクセルを$\mathbf{s_i}$に所属させましょう。そこからいくつかを$\mathbf{s_j}$に戻すわけですが、最適な数をみつけるのは$\Delta E$の中央値を見つけることと等価になります。ここでquickselectのような効率の良いアルゴリズムを使った場合、計算量は$O(\frac{m}{n})$まで減らすことができます。

また、バウンディング・サークルが重ならないセル同士の場合、ピクセルが所属するセルの交換は行わないことでも高速化することができます。さらにCCVTのステップ3をみてもわかるように、並列化は比較的容易であるため、並列化によって大幅な高速化を達成しています。それ以外にもローレベルな最適化について触れられていますがあまり重要ではないでしょう。

Capacity Constrained Voronoi Tessellation Revisited

最適化の際に使用されるコスト関数の変化が以下のように「線」に単純化されることを利用して高速化を図ります。

\Delta E=|\mathbf{p}-\mathbf{s_i}|^2-|\mathbf{p}-\mathbf{s_j}|^2=ax+by+c

$a,b,c$は省略。一般に容量制約のついたセルはボロノイ・セルではなくパワー・セルとなります。パワー・セルの境界は直線となるので上の式変形はいわば当然の結果です。

またCCVTのステップ3でセルの入れ替えが起こるのは実際隣同士であるため、隣接するセルをトラックすることで、対象となる点を減らすことができます。他にも初期の点列の生成方法など細かいことが書かれていますが、大切なのはこの二点でしょう。

BNOT

Blue Noise Through Optimal Transportの略省です。Optimal Transportを用いて最適化を行いますが、詳細を理解して、実装する時間がなかったので、そのうち加筆したいと思います。実際は次に紹介する方法が非常に簡単で、大変効率が良いのでそちらをお勧めします。(次に紹介するものを先に実装してしまい取り組むモチベーションが下がったというのが本音でしょうか。)

Weighted Linde-Buzo-Gray Stippling

Siggraph Asia 2017で発表された点描のためのアルゴリズムです。非常に実装が簡単で、結果もかなり良好です。動画に点描を適応させても、ちらつきが抑えられるなどといったメリットもあります。アルゴリズムは単純で以下の通りです。

1. 適当な数の点を生成します(1つでもよい)。
2. 点一つに注目し、その点が含まれるボロノイ・セルのキャパシティが、対象となる濃度関数とくらべ、小さければ削除、大きければ2つに分割します。この処理をすべての点に対し行います。
3. 新しい点列に対しLloydの弛緩法を適用します。
4. 2.に戻ります。

栄養のある所では細胞が分裂し、足りないところでは細胞が死んでしまう、と例えると理にかなったとても分かりやすいアルゴリズムであるといえるでしょう。著者に質問したところマルチクラス・サンプリングへの拡張も容易であるようです。

コード・サンプルはこちら。

Blue-noise Optimized Point Sets based on Procrustes Analysis

こちらもSiggraph Asia 2017で発表されたアルゴリズムです。Procrustes Analysisを用いたもの、という内容になっていますが、ドロネー三角形を描写した後、それぞれの三角形が正三角形へと近づくように最適化を行っていきます。こちらも直感的に理解しやすいシンプルな方法です。

三角形の頂点がゴム紐で正三角形の頂点から引っ張られているイメージで、正三角形の大きさあるいは、ゴム紐の強さを変えることで適応的なサンプリングもおそらく可能となるでしょう。マルチクラス・サンプリングへの適応は自明ではありません。

応用

Sampling

レンダリングにおけるサンプルの生成は、言うまでもないですね。インポータンス・サンプリングにも使えます。

Stippling

既に述べた内容ではありますが、濃度に応じてサンプルを生成する方法は点描に使用することができます。

optimized over the weekend🏴‍☠️ pic.twitter.com/tCtnZJELy1
— shinji🤔ogaki (@redqueenrender) 13 November 2017

Noise

Gabor Noiseなどカーネルをペタペタと重ねていって作るタイプのノイズでは、カーネルの配置を任意に決めることができます。その配置にこれまで紹介した方法を使うことができます。

Instancing

インスタンシングを使ってオブジェクト例えば木などを配置するのに便利です。マルチクラスをサポートしたサンプリング方法を用いれば、いろいろな種類の樹木を均等に混ぜて配置することが出来るので、自然な林や森を作ることが可能です。

Space Colonization

(image from http://gandhigames.co.uk/fractal-trees-includes-unity-implementation/)
また、葉脈などの模様を作るの方法の一つにSpace Colonizationと呼ばれるアルゴリズムがありますが、そこでもボロノイ図が活躍します。植物に関する話題はまた別の記事で紹介したいと思います。

さいごに

いろいろなアルゴリズムを紹介しましたが、ここで取り上げたものは、サンプリング、プロシージャルなシェーダやモデリングに欠かせないものなので、知っておいて損はないのではと思います。

ボロノイと点列の生成のはなし