グラフカットとは

グラフとはノード（点）をエッジ（線）でつないだ図形のことです． ノードの位置やエッジの長さの詳細は無視して相対的な繋がり具合だけを問題にします． エッジに向きを与えたグラフを有向グラフと言います． エッジには重みと呼ばれる数値を割り当てることができます． ここでは有向エッジに負でない整数を割り当てる場合を考え， この数値がその方向への最大流量を意味するものとします． ノードの中に2つの特別なノードS(湧き出し点)とT(吸い込み点)を置き， グラフの中をSからTに向かって水を流すことを考えます． このとき「最大どれだけの水を流せるのか」を見つける問題を 最大フロー問題と言います．もうそれ以上SからTへ水を流せないと言う時にはどこかのエッジが 最大流量に達していて，その最大流量に達したエッジの一群が「一枚の壁」を作って SからTまでの経路を塞いでいるはずです．すなわちこの壁によってグラフが S側とT側に分離されていることになります．

一方，そのエッジを切断することで，グラフをS側とT側に分離することができる， S側からT側への向きを持ったエッジ群を考え，これをカットと言います． カットに含まれるエッジの重みを合計した値をカット値と言います． カットはいくつも存在しますが，その中でカット値が最小のものを 最小カットと呼びます． 最小カットでグラフを分割すること， または最小カットを見つけることをグラフカットと言います．

もうそれ以上水を流せない状況ではこの最小カットのところが先の「一枚の壁」になっているはずですから， 最大フローと最小カット値は一致することになります． これを最大フロー最小カット定理と言います． できるだけ水を流そうとすると最小カットが見つかるということになります． 後で詳しく説明しますがグラフカットのアルゴリズムはどれも最大フローを流すことで最小カットを見つけています．

組み合わせ最適化問題の評価値をグラフのカット値として表現できる場合には グラフカットにより最適解を求めることができます． 他の手法では最適化解（最適を目指して頑張った解）しか求まらないのに対し， 本当の最適解が見つかります．

残余フロー

まずどのアルゴリズムにも共通する残余フローという考え方について説明します． あるエッジにどれだけの水量を流せば全体としての最大フローが得られるのかは あらかじめ分かっているわけではありません．とりあえず流してみるしかありません． とりあえず流すことで当然， あとどれだけそのエッジに水を流せるかの値は最大流量よりも小さくなります． この「あとどれだけ流せるか」の値を残余フローと言います． ここまでは当たり前に感じますが，重要なのは， とりあえず流した量と同じだけの量を反対向きのエッジの残余フローに加えることです． 反対向きに流すということは最初にとりあえず流した量を減らす（キャンセルする）ということなので， とりあえずで流しすぎて失敗するということがなくなるわけです． バックトラック（キャンセル）の機能を持つこの残余フローのおかげで 一見グリーディ（欲張り法）に見えるアルゴリズムで最大解が求まります．

グラフカットアルゴリズム

グラフカットアルゴリズムには2つの基本的な考え方が存在します． 1つ目は増大パス法（augmenting path method）， 2つ目はプッシュリラベル法（push relabel method）です．

増大パス法は残余フローが残っているSからTまでの経路を見つけては流すという処理を そのような経路が見つからなくなるまで行うというものです． 経路の見つけ方や流し方に様々なバリエーションがあります． 記事1で紹介したBKアルゴリズム(Boykov-Kolmogorov algorithm) も増大パス法に分類されます．

プッシュリラベル法では実際の水の流れにヒントを得た過剰フロー（over flow）というアイデアを使います． 浸食などを考えない時，実際の水の流れがどのように決まって行くかというと， 水は高い方から低い方へ流れるのですが， その高さ（水位）は以下に説明する仕組みにより自動的に決まって行きます． 上流から流れて来る量よりも下流へ流れ出す量が少なければ， そこに水が溜まって，その場所の水位は上昇することになります． そしてこの水位の上昇は上流からの流入を減らし，下流への流出を増大させます． 逆に上流から流れて来る量よりも下流へ流れ出す量が多ければ， 溜まっている水が減って，その場所の水位は下降することになります． そしてこの水位の下降は上流からの流入を増やし，下流への流出を減少させます． この水位の変化により流れが調整されて最終的に定常的な流れに到達します． 水位の差が大きくなると流量が増え， 水位の差が小さくなると流量が減るという関係が流れを調整する訳ですから， 最大流に達して水位の差が増大しても流量が増えなくなると， 水位の差が拡大して，自然の水流では滝となります． ちょうどEFT(Field Effect Transistor)のピンチオフ点のように． そしてこの滝が最小カットに対応することになります． そして滝の下流ではどの場所も流れて来た量が流れ出しているだけでなく， 流れを増やす余裕を持っています．

過剰フローとは流出量よりもどれだけ流入量が多いかの値です． 過剰フローがあるノードは水位が上昇している場所に対応します． 過剰フローを隣のノードに渡すということは，そのノードへの流出量を増やして 自分は流入量と流出量をバランスさせることを表します． 水は高い方から低い方にしか流れませんので， 過剰フローは高さが低い方にしか渡せません． もちろん残余フローが残っている方向にしか渡せません． 過剰フローの一部を渡すこともできます． 過剰フローを渡せないときは自分の高さを高くします． SからTに向かって過剰フローを渡していくと 滝の下流には過剰フローが無いようにできます． 滝は残余フローがゼロとなったエッジ群で，もうそこを過剰フローは通過できません． 過剰フローを渡す操作をpush, 過剰フローを渡せなかったとき自分の高さを高くする操作をrelabelと言います． このpushとrelabelにより滝とその下流に過剰フローが存在しない状況を作り出せば 滝の位置が確定してグラフカットが完了します．

高低差が大きいほど多く流れると言うのではなく， 高低差があれば流れるとしたところが実際の水の流れとは異なります． これは滝を見つけるのが目的であるからです． ですので高さはどちらが高いかが判定できれば十分なので整数で表し， 1だけ大きければ十分であるとします． 過剰フローを持つ全てのノードに対して， pushかrelabel（詳細はすぐ後で説明します）を1回適用することを単位操作として， これを繰り返す場合には， 流れを増やせるところの高低差は0または1となります（証明は難しいので省略）． 従って，Tの高さを0として高低差1で全ノードを繋ぐと， 最も高いノードの高さは全ノード数-1となります． 高低差が0のところもありますので，これが流れを増やせる範囲の考えうる最高の高さです． そこで全ノード数を流れを増やせない滝の上側の高さに割り当てます． 単位操作を開始する前の高さの初期値はS以外は0，Sは全ノード数とします． すなわちSは滝の上側にあるものとします． relabelは「周りのノードの高さを調べて一番低いものより1だけ大きく」とする 文献と「単に自分の高さを1だけ大きくする」とする文献があります． 「周りのノードの高さを調べて一番低いものより1だけ大きく」とする方が証明がやり易そうですが， 「単に自分の高さを1だけ大きくする」でも正しいことが証明できるようです． 実際のプログラムでは「単に自分の高さを1だけ大きくする」の方が速くなります．

pushとrelabelだけでもグラフカットは完了するのですが， 速度を向上させるためにglobal relabelingやgap heuristicという操作を組み合わせて使用します． global relabelingはTから幅優先探索を行ってその時点のすべてのノードの高さを一気に決めてしまう操作です． あるノードから残余フローが残っているエッジを辿ってTへ辿り着けるとき， その最短のエッジ数をそのノードの高さとします．これを高速に行うアルゴリズムが幅優先探索です． Tへ辿り着けないノードは滝の上側にあるわけですから， 以後の操作から除くことができます． 50単位操作程度毎にglobal relabelingを行います．

単位操作を繰り返しているとき， 流れが到達していないところの高低差は0， 流れが存在してその流れを増やせるところの高低差は1です． 従ってこの両方の範囲で高さとノード数のヒストグラムを作ると， ある高さ以上でノード数が0となるのは当然として， それよりも低い高さには少なくとも1個のノードが存在します． 従って全体を対象とする高さとノード数のヒストグラムの中に， ノード数が0でない高さに挟まれてノード数が0の高さ（ギャップ）が出現した場合には， 高さがギャップより大きいノード群は流れを増やせる流れから切り離されている， すなわち滝の上側にあることが分かります． 滝の上側のノードは以後の操作から除くことができます． これをgap heuristicと言います．

並列処理向きグラフカットアルゴリズム

プッシュリラベル法では， 過剰フローを渡されたノードと過剰フローを渡せなかったノードが， 次の単位操作の主役となるわけですから， CPUではキューを使って主役を管理すれば効率を上げることができます． ところがGPUやFPGAのような並列処理向きのデバイスを使う場合は キューでは並列性が生かせません． どうしても全ノードに対して主役かどうかの判定を行うことになります． そうであればglobal relabelingとpushをパイプライン的に行うことで， この全員参加を有効に使えるのではと考え思い付いたのが， Wave Front Fetchアルゴリズムです．研究ノートによれば2014/10/16に思い付き， 2014/11/04に命名したようです．

基本的な考え方は次のようになります． 流れの到達点であるTは新しい時刻を， T以外のノードは自分の時刻を， 自分に向かう残余フローで繋がっているノードへ送り出します． この結果T以外のノードはいくつかのノードから時刻を受け取ることになりますが， 最も新しい時刻で自分の時刻を更新するとともに， この時刻を送ってきたノードに過剰フローを渡します． この処理を全ノードが実行することを単位操作と言うことにして， 滝の下で過剰フローがなくなるまでこの単位操作を行います． 単位操作毎に様々な時刻の波が一斉にT側からS側に伝わっていく様子や， 異なる時刻の波がぶつかって古い時刻の波が消滅する様子や， 新しい波の先頭が過剰フローを取り込む様子や， 過剰フローを受け取ることで残余フローが無くなってもう波を伝えられなくなる様子 がイメージできると思います． Wave Front Fetchアルゴリズムという名前はこのイメージを表したものです． このアルゴリズムは過剰フローという考え方を使っていますのでプッシュリラベル法に分類されますが， 役割の証明が難しいrelabelを使わない分だけ簡単となっています．

Wave Front Fetchアルゴリズムの詳細

Wave Front FetchアルゴリズムをNVIDIA社のCUDAで記述すると

• tim_read(S, H, W) read the time of node(S, H, W)
• tim_write(S, H, W, V) set V to the time of node(S, H, W)
• ovf_read(S, H, W) read the overflow of node(S, H, W)
• ovf_add(S, H, W, V) add V to the overflow of node(S, H, W)
• edg_read(S, H, W, D) read the flow of edge(S, H, W, D)
• edg_add(S, H, W, D, V) add V to the flow of edge(S, H, W, D)

__global__ void wave_init(void) {
  ///////////////////////////////
  int total_id = blockDim.x * blockIdx.x + threadIdx.x;
  int S  = total_id / size_hw;
  int sa = total_id % size_hw;
  int H  = sa       / size_w;
  int W  = sa       % size_w;
  ///////////////////////////////
  tim_write(S, H, W, 0);
}

__device__ void push(int S, int H, int W, int nS, int nH, int nW, int D, int R) {
  int tim = tim_read(S, H, W);
  int Tim = tim_read(nS, nH, nW);
  if (Tim <= tim) return;
  int edg = edg_read(S, H, W, D);
  if (edg == 0) return;
  int ovf = ovf_read(S, H, W);
  if (ovf > 0) {
    bool q = edg > ovf;
    int pp = (q)? ovf: edg;
    ovf_add(nS, nH, nW, pp);
    ovf_add(S, H, W, -pp);
    edg_add(nS, nH, nW, R, pp);
    edg_add(S, H, W, D, -pp);
    if (q) tim_write(S, H, W, Tim);
  }
  else tim_write(S, H, W, Tim);
}

__global__ void wave(int tim) {
  ///////////////////////////////
  int total_id = blockDim.x * blockIdx.x + threadIdx.x;
  int S  = total_id / size_hw;
  int sa = total_id % size_hw;
  int H  = sa       / size_w;
  int W  = sa       % size_w;
  ///////////////////////////////
  if (ovf_read(S, H, W) < 0) {
    tim_write(S, H, W, tim);
    return;
  }
  if  (S!=size_s-1)                 push(S, H, W, S+1, H,   W,   0, 9);
  if ((S!=size_s-1)&&(W!=size_w-1)) push(S, H, W, S+1, H,   W+1, 5, 7);
  if ((S!=size_s-1)&&(W!=0))        push(S, H, W, S+1, H,   W-1, 6, 8);
  if                 (W!=size_w-1)  push(S, H, W, S,   H,   W+1, 4, 3);
  if                 (W!=0)         push(S, H, W, S,   H,   W-1, 3, 4);
  if                 (H!=size_h-1)  push(S, H, W, S,   H+1, W,   2, 1);
  if                 (H!=0)         push(S, H, W, S,   H-1, W,   1, 2);
  if ((S!=0)       &&(W!=size_w-1)) push(S, H, W, S-1, H,   W+1, 8, 6);
  if ((S!=0)       &&(W!=0))        push(S, H, W, S-1, H,   W-1, 7, 5);
  if  (S!=0)                        push(S, H, W, S-1, H,   W,   9, 0);
}

void wave_cut(int loop) {
  wave_init<<< SIZE_SHW/256, 256 >>>();
  for (int time = 1; time <= loop; time++) {
    wave<<< SIZE_SHW/256, 256 >>>(time);
  }
}

関数push()を見ると， まず自分の時刻と隣の時刻を比べて隣が新しくないときは意味がありませんので，returnしています． 時刻は値が大きいほど新しい時刻であることを表します． 次に，隣がいくら新しくても隣へ向かう残余フローがなければ意味がありませんので，returnしています． 残余フローが残っていて過剰フローも持っているというときには， 過剰フローの全体か一部を隣に渡すことができます． しかし残余フローを超えて渡すことはできませんし， 過剰フローを超えて渡すこともできません． そこで残余フローと過剰フローのどちらが小さいかを調べて， その量の過剰フローを隣に渡し，隣へ向かうエッジの残余フローを減らし， 逆向きの残余フローを増やしておきます． 過剰フローを隣に渡すということは隣の過剰フローを増やし， 自分の過剰フローを減らすということです． 過剰フローを渡した結果，残余フローが無くなってしまった場合は， 経路が切れてしまったということですから， 間違った情報が伝わって行かないよう自分の時刻の更新は行いません． 残余フローが残っている場合と過剰フローが最初から無かった場合は自分の時刻の更新を行います．

関数wave()を見ると， 最初の部分は自分が一斉に起動されたスレッドの内のどれに該当するのかを判定する部分です． その後の過剰フローが負かどうかを判定している部分は， Tへの残余フローを持つかどうかを判定しており， Tへの残余フローを持つ場合は，Tへの経路があるということなので自分の時刻を新しい時刻に更新しています． このあとのpush()のなかで，Tへの残余フローを持つノードへ過剰フローが渡されると，過剰フローの負の度合いが減らされます． その結果もし過剰フローが正となってしまった場合は，Tへの残余フローがなくなって， 他に渡すべき過剰フローを持つ普通のノードになります．このように過剰フローに負の値も使うと SとTを除いたノードだけでデータ構造を作ることができます．

Wave Front Fetchプログラムの配布

配布データとプログラムは迷惑かもしれませんが

• 3_rec_L.png
• 3_rec_R.png
• rectify.cpp
• makefile
• graph_cut.cu_wave_push_only
• graph_cut.cu_push_relabel
• graph_cut.cu_wave_push
• graph_cut.cu_with_shfl
• graph_cut.cu_without_shfl

これらのgraph_cut.cu...はgraph_cut.cuという名前のコピーを作って

make do
>,fffjjjq

起動形式は

./rectify right_image left_image max_disparity min_disparity penalty inhibit divisor option LOOP

結論

配布プログラムを実行してみると分かるのですが，inhibitをちゃんと入れた10エッジ版で終了判定まで行うと， push relabel+global relabelingが最も速いようです． それではWave Front Fetch方式の利点は何かというと，コードが単純なことであると思います． push relabel方式は，push relabelの部分は単純なのですが，global relabelingがないと全く速度がでません． Wave Front Fetch方式ではこれだけでほぼ同程度の速度が得られます． Wave Front Fetch方式の終了判定をどうするかは悩ましい問題ですが， 一つの考え方として，グラフカットが完了していなくてもある回数で打ち切るという方法があります． 最初のかなり早い段階でほとんどのフローは流れてしまっており， ほとんどの時間を最後のほんの少しのフローを流すために使っています． ほとんどの場合，途中で打ち切ってもステレオビジョンの表示はほとんど劣化しません． これをどの程度に選ぶかを明らかにするという研究の方向もあり得るのではと思います．

ノード当たり32bitとした6エッジ版がいまのところ最速ですが， GTX1080やGTX1080Tiなどの最高性能のGPUを使っても， リアルタイムで使うにはもう一歩という状況です． 余裕を持ってリアルタイム処理を行うには， ピクセル数を減らしたり，グラフカットを階層的に行うなどの工夫が必要です． 将来自動車の自動運転などに応用するということであれば， GPUの低消費電力化がどれほど進むかが大問題です． やはりFPGA化やASIC化が必要かもしれません． その時にWave Front Fetch方式の単純性が役立つように感じます．

inhibitを無くすことの意味

32bit6エッジ版ではinhibitを無くしてしまっているのですが，その意味について整理しておきます． pixel-disparity modelでは隣のピクセルであっても大きく変化する厄介な視差(disparity)を扱うために， 視差の変化に比例(linear)するpenaltyではなく，打ち切り比例(truncated linear)のpenaltyが使われています． これは「視差の変化が大きいところでpenaltyの効果を減らす」ものです． gaze_line-depth modelでは奥行き(depth)の変化は 視差と違って最大でも1ですから，変化に比例するpenaltyどころが，2以上の変化を許さないinhibitも与えています． これは「奥行きの変化が大きいところでpenaltyの効果を増やす」ものです． このようにpixel-disparity modelとgaze_line-depth modelでは変化が大きいところの考え方が逆になっています． inhibitを与えない場合は変化に比例するpenaltyだけを与えることになりますが， 「奥行きの変化が大きいところでpenaltyの効果を増やす」ことに違いはありません．