Boltzmanです.
以前記事に書いた通り, ラベリング処理の実装が難航しています.
正直全然進捗がありませんが, やや行き詰ってきたところもあるので, ここで一度風呂場でアヒルに向かって話すつもりで記事にしようかと思います.
ちなみに”Rubber duck debugging”という言葉はchokopan氏に教わりました.

まず最も愚直な並列化はどのようなものだろうか.
私はこれにさえ考え至らなかったのだが, UMU氏に相談したところ, 「愚直なやり方でよければ」という枕詞とともに即答してくれた.
以下では二値化後に値が0で(黒)である領域についてラベリングを行うものとする.
1.まず, 各スレッドに対してピクセルをひとつづつ割り当てる
2.それぞれのスレッドで, 自分のピクセルの色を見て, 白であれば0を, 黒であればスレッド番号を新しいテクスチャに格納する
3.その後, 各ピクセルについて自分と, 自分の周囲8方向の隣接するピクセルと数字を比較し, その中で最も大きい値を新たに自分の値としてテクスチャに格納する
4.3.をたかだか有限回繰り返せばラベリングが完了する.
上の処理がたかだか有限回で終了する理由は, 同一の領域に含まれる任意の2ピクセルについて, 領域内でのその二点の最短距離は必ずテクスチャの面積より小さくなるためである. (もっと言えば, 始めの二値化画像において, 黒部分の面積をSとすると, 3.の操作はS回以下で構わない. )

しかしこの実装方法には現実的な問題点がある.
まず, この手法でラベリングを行った場合, 可能なラベル番号がテクスチャのピクセル数存在する.
現在私の実装では1024×720を生のまま使っているので, 可能なラベル番号は1024*720ある.
最終的には領域を面積順にソートするのであるが, この際に全てのラベル番号に対してソートするとなると,
各領域の面積を算出する時やソートのアルゴリズムの中でshared-memoryに1024×720の長さを持つ配列を確保する必要が生じる.
これはあまりに多いのでこのラベル番号から使われていない不要なものを捨てる必要がある.
この問題がこの手法の最大の問題点といえる.

考えるべきこととしては四つ
・冗長なラベル番号の削減
・3.の操作の実装方法
・3.の操作の回数の明確な上限
・そもそもほかにもっといい方法がないのか

まず, 不要なラベル番号の捨て方について考える.

はじめに, 実はラベル番号の数は初めから上で示した数の1/4でよいことがわかる.
初めから(2m, 2n), (2m+1, 2n), (2m, 2n+1), (2m+1, 2n+1)には同じ番号を振っても問題ない.
なぜなら, この4ピクセルにどのようなパターンで黒が割り振られても必ず隣接するためである.
必ず隣接するのであれば初めから同じラベル番号を割り振ってしまって差し支えない.

さらに, 1度3.の操作を行った後, 各スレッドにピクセルを4つずつ割り当て(始めに同じ番号を割り当てた4つ), 始めに割り当てた番号と比較する. この時, スレッド内の一つ以上の0でないピクセルについて, 3.を一度行った後に割り当てられている番号が始めに割り当てた番号と異なる番号だった場合には, そのスレッド番号は最終的な状態においてラベル番号として存在しないといえる.
ただし, 以下の場合については一度3.を行っただけではわからない.

リソースの消費を抑えるという意味では, この結果を受けて必要最小限にメモリを確保したいが, メモリの確保を命令するのはCPU側なので
CPU-GPU間で通信をしなければならず, 実行速度の大幅な低下が予想される.

ここからはその後の処理について検討する.

ソートのアルゴリズムにはbitonic sortを用いれば, ブロック一つ当たりのshared memoryの大きさは一定で構わないので, 問題は面積を求める際のsumである.
各ラベルの面積を求める作業は, 横軸をラベル, 縦軸を要素数とするヒストグラムを生成する作業に等しい.
このとき, 高速化を図るために, ブロック内で一度shared memoryにヒストグラムを生成し, そのあとにそれらをグローバルなメモリにまとめるという作業が必要になる.
この際, shared memoryに用意する必要があるのは, そのブロックの扱うテクスチャの範囲に残りうる最大のラベル数分だけで構わない.
理論的に残るであろうラベル番号の数は, 少なくとも(一つのブロックの中に存在しえる最大のラベル番号数)x(ブロックの総数)よりも小さい.
同一ブロック内に存在しえるラベル番号の数が最大となるときの領域は少し考えたところでは以下のようなものが思い浮かぶ.

これが正しければ, 例えば一つのブロックが32×32を扱う場合, shared memoryに用意する必要があるのはデータ用の16×16=256+インデックス用に256個分となる.

他の手法においては, 使用されているラベル番号をグローバルに記憶しておき, それを参照するという形をとっているが, この方式だと一つの場所にアクセスが集中する上に, 一度に一つのスレッドしかアクセスできないようにロックする必要があるため, 遅延が見込まれる.
また, この手法においてもあらかじめ仮ラベルの数を推測することはできないため, やはり領域を面積順にソートする際に今手法と同様の問題が生じることと思われる.

次に, 3.の操作をどのように実装するかについて考察する.
以下に三つのやり方を用意した.
①自分と, 自分に隣接する8つのピクセルについて番号を比較し, 最大のものを新たに自分の番号として採用する.

上の方法を2つに分割する.
②自分と, 自分に隣接する4つのピクセルについて番号を比較し, 最大のものを新たに自分の番号として採用する.

次は上の二つとは少し異なるやり方
③全体を1マス2×2のグリッドに区切って(区切り方は4通りある), そのグリッド内の各ピクセルについて, グリッド内の最大の番号を新たに採用する.

①に比べて, ②③では3.の操作をする回数が増えると考えられるが, メモリのアクセス集中が緩和できる分, 実行速度は上昇すると考えられる.
アクセスの重複は①で9, ②で5, ③で0回である.

上の3つの方法に関して, 具体的に3.の操作の実行回数がどれくらいになるのか実際に計算してみた.
プログラムはGPUを使わず, Cで記述しているため, 実行速度に対する検討はできない. 矢上の計算機サーバーにはC++/Fortranのコンパイル環境しかないため, Cで実装して回数のみを見ることにした.
二値化画像をランダムに生成してその画像を①~③の方法でラベリングするという作業を繰り返し1万回行い, ラベリング終了までの平均回数と最大回数を見る.

32×32画像の場合
①：平均 37.6 最大 64
②：平均 75.3 最大 129
③：平均 78.1 最大 170

128×128画像の場合
①：平均 157 最大 177
②：平均 315 最大 354
③：平均 299 最大 414

1024×1024画像の場合
①：平均 1257 最大 1286
②：平均 2515 最大 2572
③：平均 2321 最大 2468

実装の仕方が適当だったので1024×1024に関しては計算に50時間以上かかった.

どの手法を用いてもオーダーは線形になることがわかる.
最大と平均のいずれにおいても, 画像の幅が大きくなると③(Method3)は②(Method2)よりも少ない3.の試行回数で済むことがわかる.
メモリ参照の衝突分を考慮しても, ラベリングに要する時間は③(Method3)のほうが②(Method2)よりも短くなることが予想される.
①(Method1)と②(Method2)を比較すると, ②(Method2)では①(Method1)の調度倍程度の回数を要することがわかる.
マスクを分割した際, 二つのマスクについて同じ回数処理を行うよう設計したが, おそらく前半では上半分のマスクの試行が, 後半においては下半分のマスクの試行が, それぞれほとんど無駄になっているために, 倍の回数がかかっているのではないかと思われる.
現在の実装では, 上半分のマスクとした半分のマスクを交互に適用しているが, これを, 始め下半分のマスクから始めて, ラベルに変化がなくなるまでこれを続け, 次に上半分のマスクを繰り返し適用して, こちらも変化がなくなるまで続け, 最終的にどちらのマスクでも変化がなくなるまで続ける, という形に変えれば, もっと回数は削減できるのではないかと考えられる.

さらに, ①の場合において3.の操作の必要回数の明確な上限について考察する.

3.の操作の必要回数は, グラフの問題に帰結すると考えられる.
すなわち,

上図のようなグラフAを考える.
このグラフは, 横w個, 縦h個の頂点を持ち, 頂点の総数はw*hである.
ある頂点は, 自分の周囲8方向の頂点と必ず辺で接続されている. それ以外の頂点とは辺を持たない.
このグラフから任意の頂点とその頂点に接続する辺を除いた部分連結グラフGを考える.
(任意の頂点とその頂点に接続する辺を除いた場合非連結グラフができる可能性もあるが, そのような場合は考えない)
Gの任意の2頂点u, vを考える. (u,v∈V(G))
u, vを始点/終点とする単純道のうち, 最も頂点数の少ないもの(最短経路)に含まれる頂点数をnとして, n-1をu-vの「距離」と呼ぶことにする
このような「距離」のうち, 最大のものはいくらか.

グラフの最短経路・最長経路問題というとNP完全がちらつくが, この問題は果たして解けるのだろうか…
上の測定の結果から推定すると, この最大「距離」もグラフの辺の長さに比例すると類推できる.

正確な最大値が算出できなくても, 処理前のラベルと処理後のラベルで変化があるかどうかを確認すれば終了しているかどうかは判定できるが, その場合動作が安定しない可能性がある.
また, 処理の繰り返しはCPU側から命令する必要があるので, (たかだか1bitで事足りるが)マイフレームCPU-GPU間の通信が発生する.
繰り返しをホスト側のコードで実装する必要があるのは, ブロック間の同期命令を(少なくともHLSLでは)デバイス側に記述できないためである.

最後により良い方法の検討について
正直私は理論はからっきしなのでよくわからない.
知人から局所的にラベリングをしてからそれをマージする方法をすすめられた.
確かにほかの領域分割の手法ではマージをする方法がある.
しかしこの手法は領域分割のみに言及しており, 各領域に固有な番号を振る方法については言及がなく, 実際そこまでの実装は難しいように思う.

考察にもならない考察はこのあたりでいったん区切り, とりあえずなんでもいいから実装してみようと思います. 環境はPyOpenCLにすると思います.

Posted on: 2016年8月19日, by : Lait-au-Cafe