3. ステレオビジョンプログラムの配布

グラフカットを使用したステレオビジョンのプログラムを作りました. 共通のファイルとして scene1.row3.col1.ppmscene1.row3.col3.ppmmakefileを使います. scene1.row3.col1.ppmとscene1.row3.col3.ppmは記事2で紹介したTSUKUBAステレオ画像の内の真値データに対応する左右2枚の画像です. GPUとCUDAを使う場合は以下のmakeのターゲットをdis_cからdis_gに,dep_cからdep_gに,3d_cから3d_gに変更してください. またこの記事で配布するプログラムはすべて自由に使ってもらって構いません (the MIT license).

Graph Cut Program for Stereo Vision

graph_cut.cppgraph_cut.cuが グラフカットプログラムの本体です.graph_cut.cppはCPUによるグラフカット, graph_cut.cuはGPUによるグラフカットになっています. どちらも3D Grid Graphに対してPush Relabel + Global Update方式でグラフカットを行っています. Global Update部分はCPUではキューによるアルゴリズムを使い, GPUではシェアードメモリによるアルゴリズムを使っています. これらは全く異なるアルゴリズムです.一方Push Relabel部分はCPUとGPUで同じアルゴリズムを使用しています.

グラフカットはグローバル変数の

  • int SIZE_S;
  • int SIZE_H;
  • int SIZE_W;
  • int SIZE_HW;
  • int SIZE_SHW;
に値を設定してからgraph_cut()関数を呼び出すことで行います. graph_cut()関数からは呼び出し元のcost()関数がコールされてデータが準備され, グラフカットが実行されます. その後呼び出し元のcut()関数がコールされて結果が渡された後に, graph_cut()関数が終了します.

グラフカットの対象となる3D Grid Graphは次の図の様なステレオビジョン専用の10隣接枝グラフとなっています. ただし図は奥行き方向Sと横方向Wだけが見える平面図で,縦方向Hは見えていません. ソースノードSからシンクノードTへの方向を奥行き方向と言っています.

グラフの構造とgaze_line depth modelとの対応
この10隣接枝は
  • 0: (reverse=9) (init_val=data_cost) front
  • 1: (reverse=2) (init_val=penalty_h) up
  • 2: (reverse=1) (init_val=penalty_h) down
  • 3: (reverse=4) (init_val=penalty_w) left
  • 4: (reverse=3) (init_val=penalty_w) right
  • 5: (reverse=7) (init_val=zero) front right
  • 6: (reverse=8) (init_val=zero) front left
  • 7: (reverse=5) (init_val=inhibit_b) back left
  • 8: (reverse=6) (init_val=inhibit_b) back right
  • 9: (reverse=0) (init_val=inhibit_a) back
の0から9の方向番号で区別されています.枝には重みを与えますが, penalty_w,penalty_h,inhibit_a,inhibit_bはグラフ全体で共通です. これらはgraph_cut(int penalty_w, int penalty_h, int inhibit_a, int inhibit_b)関数の呼び出し時に引数で与えます. data_costは枝ごとにgraph_cut()関数から呼ばれるcost(int S, int H, int W)関数のリターン値が設定されます. inhibit_bはgaze_line-depth modelでのみ使われます.使わないときは重みを0とします. 重み0の枝は枝そのものが存在しない場合と等価になります. グラフカットの結果はgraph_cut()関数から呼ばれるcut(int S, int H, int W)関数により与えられます. inhibit_aは2Dのサイトにラベルを割り付ける問題を3Dのグラフカットで解くためのもので, ソースノードSの方向を向いた枝の重みです. このinhibit_aが小さいとグラフカットの結果をcut(int S, int H, int W)関数で表現できなくなります. cut(int S, int H, int W)はサイト(H, W)にラベルSを与えるという意味です. サイトラベリングとならなかった場合でもグラフカットは正しく行われていますが, 結果を正しく使うためにはプログラムの変更が必要です. またサイトラベリングにならなかったというエラーも表示されません. reverse値はいったん流したフローを流さなかったことにするために存在します. これがあるために手続きとしてはバックトラックを行いませんが,グリーディとならずに,最大フロー最小カットが得られます.

3D Grid Graphのサイズはグローバル変数の

  • int SIZE_S;
  • int SIZE_H;
  • int SIZE_W;
で与えますが,Sは奥行き方向,Hは縦方向,Wは横方向を表しています. SIZE_Wが384,SIZE_Hが288,SIZE_Sが28のとき, 横方向に384個,縦方向に288個,奥行き方向に28個のノードがあることになりますが, これに最も手前のソースノードSと最も奥のシンクノードTの特別な2個のノードを加えたものが全体のノードとなります. サイトラベリングとして見るときは,横方向に384個,縦方向に288個のサイトに, 0から28までの29個のラベルを割り当てることができます.サイトはどこで枝がカットされたかに対応するので, サイトは個数が1だけ大きくなります.cost(int S, int H, int W)関数の第1引数Sも0からSIZE_Sまでの値をとります. これに対して第2引数Hは0からSIZE_H-1まで, 第3引数Wは0からSIZE_W-1までの値をとります.

記事2で紹介したTSUKUBAステレオベンチマークの真値視差画像では 最大視差(max_disparity)が28,最小視差(min_disparity)が10ですから, SIZE_S = max_disparity - min_disparity = 18とするとラベルは0から18までの19個となります. これが10から28に対応します.このようにSIZE_Sはその値そのものはノードの数なのですが, ラベルとして見るときは最大ラベルを表し,最小ラベルは0なので,ラベルの数はSIZE_S + 1となります. プログラムを作っているときはいつもこのことに気を付けていなければなりません.

Pixel-Disparith Model

disparity_stereo.cppは pixel-disparity modelによって視差画像ファイルを生成するプログラムです.

make dis_c
とすると視差画像ファイルdisparity14.pngが生成されます.
disparity14.png
disparity_stereo.cppは100行程度の簡単なプログラムですのでぜひ読んでみてください.

pixel-disparity modelであることは以下のcost()関数

int cost(int D, int H, int W) {
  D += Min_depth;
  int RR, RG, RB;
  int LR, LG, LB;
  int LW;
  if (W+D > SIZE_W-1) LW = SIZE_W-1;
  else                LW = W+D;
  call_rgb(rframe, H, W, &RR, &RG, &RB);
  call_rgb(lframe, H, LW, &LR, &LG, &LB);
  int d_R = (LR>RR)? LR-RR: RR-LR;
  int d_G = (LG>RG)? LG-RG: RG-LG;
  int d_B = (LB>RB)? LB-RB: RB-LB;
  return d_R + d_G + d_B;
}
で,左画像(lframe)の方が視差(D)だけ右のピクセル(LW=W+D)を使っているところから分かります.

実行は

./プログラム名 右画像 左画像 最大視差 最小視差 penalty値 inhibit値 sacle値
となっています.最大視差と最小視差の範囲で最適解が探索されます. 範囲を広くするとそれだけ時間がかかりますし,結果が悪くなる可能性も増えます. 画像に写っている対象が,前後でどの範囲にあるかが分かっているときは, できるだけぴったりの最大視差と最小視差を指定してください. inhibit値は十分に大きければ結果に影響しません. penalty値は結果が良くなる値がありますが, その値は試行錯誤で決めるしかありません. 出力される視差ファイルはdisparity<penalty値>.pngとなります.

Gaze_line-Depth Model

depth_stereo.cppは gaze_line-depth modelによって視差画像ファイルを生成するプログラムです.

make dep_c
とすると視差画像ファイルpen_14_inh_1023.pngが生成されます.
pen_14_inh_1023.png
depth_stereo.cppも100行程度の簡単なプログラムですのでぜひ読んでみてください.

gaze_line-depth modelであることは以下のcost()関数

int cost(int D, int H, int W) {
  D += Min_depth;
  int RR, RG, RB;
  int LR, LG, LB;
  int RW;
  int LW;
  if ((W-D < 0) && (W+D > SIZE_W-1)) {
    printf("Error\n");
  }
  else if (W-D < 0) {
    RW = 0;
    LW = 2*W;
  }
  else if (W+D > SIZE_W-1) {
    LW = SIZE_W-1;
    RW = SIZE_W-1-2*(SIZE_W-1-W);
  }
  else {
    RW = W-D;
    LW = W+D;
  }
  call_rgb(rframe, H, RW, &RR, &RG, &RB);
  call_rgb(lframe, H, LW, &LR, &LG, &LB);
  int d_R = (LR>RR)? LR-RR: RR-LR;
  int d_G = (LG>RG)? LG-RG: RG-LG;
  int d_B = (LB>RB)? LB-RB: RB-LB;
  return d_R + d_G + d_B;
}
で,右画像(rframe)はdepth(D)だけ左のピクセル(RW=W-D), 左画像(lframe)はdepth(D)だけ右のピクセル(LW=W+D) を使っているところから分かります.

実行は

./プログラム名 右画像 左画像 最大視差 最小視差 penalty値 inhibit値 倍率 scale値
となっています.最大視差と最小視差の範囲で最適解が探索されます. 範囲を広くするとそれだけ時間がかかりますし,結果が悪くなる可能性も増えます. 画像に写っている対象が,前後でどの範囲にあるかが分かっているときは, できるだけぴったりの最大視差と最小視差を指定してください. inhibit値は十分に大きければ結果に影響しません. penalty値は結果が良くなる値がありますが, その値は試行錯誤で決めるしかありません. 出力される視差ファイルはpen_<penalty値>_inh_<inhibit値>.pngとなります. 倍率は画像をその倍率分拡大して処理を行い,結果としてサブピクセル処理を行うものです. しかし倍率2でも8倍のグラフサイズとなりますし,画像の拡大では情報が追加される訳ではないので, 2を超える倍率指定は意味が無いようです. 記事2で注視線の数をピクセルの数と同じ(倍率1の時)にすると, 奥行き方向の分解能が半分になると指摘しましたが, この倍率を2にして初めて奥行き方向の分解能がpixel-disparity modelと同じになります. pen_14_inh_1023.pngは倍率1で作成したものです.

3D Stereo Vision Program

3d_stereo.cppは gaze_line-depth modelによるステレオビジョン処理の結果を3D表示するプログラムです.

make 3d_c
とすると3D表示プログラムが起動されます. 1文字キーコマンドを受け付けます.1文字キーコマンドの一覧は?キーで表示されます. 基本的な使い方は,プログラムが起動したら,","でグラフカットを実行し, それを"h","j","k","l","f","n"で視点を変えて眺めてください.終了は"q"です. figureモードではステレオカメラの左右の画像やカメラ中心,注目視線(gaze_lineではありません)が表示されます. 以下は1文字キーコマンドの一覧です.
1文字キーコマンド一覧
文字機能
qquit
?help
etoggle figure mode
mleft picture
,stereo picture (do graph_cut)
.right picture
sdepth(disparity) to front
tdepth(disparity) to back
heye move left
jeye move down
keye move up
leye move right
feye move far
neye move near
!decrease base line
@increase base line
#decrease focal length
$ncrease focal length
%decrease pixel size
^increase pixel size
zdecrease penalty
xincrease penalty
Zdecrease inhibit
Xincrease inhibit

プログラムの起動は

./プログラム名 右画像 左画像 最大視差 最小視差 penalty値 inhibit値 倍率
となっています.最大視差と最小視差の範囲で最適解が探索されます. 範囲を広くするとそれだけ時間がかかりますし,結果が悪くなる可能性も増えます. 画像に写っている対象が,前後でどの範囲にあるかが分かっているときは, できるだけぴったりの最大視差と最小視差を指定してください. inhibit値は十分に大きければ結果に影響しません. penalty値は結果が良くなる値がありますが, その値は試行錯誤で決めるしかありません. penalty値とinhibit値はプログラムの動作中に1文字キーコマンドで変更できますので, 変更してから再度","でグラフカットを実行すると違う結果を見ることができます. 倍率は画像をその倍率分拡大して処理を行い,結果としてサブピクセル処理を行うものです. 2倍の倍率指定でもずいぶん表示の密度が高まったように感じられます.

背景知識

ステレオマッチングは左右の画像の対応する点のずれ, すなわち視差,を検出して,この視差から対象の空間配置を推定する技術です. 視差から空間位置を推定するのは三角測量の計算をすれば良いだけですので, どうやって対応する点を見つけるかが問題の本質であると考えられてきました. ですので左右のステレオ画像から空間位置を推定する問題はステレオマッチングと呼ばれています. 教科書の章立てや項目名もステレオビジョンではなくてステレオマッチングとなっていることが普通です. 2つの画像から対応する点を見つける2D(画像)処理の後に,空間の配置を決定する3D処理を行う, すなわち別々の処理として扱ってきたということです. 対応を見つけようとするステレオマッチングという用語や視差という概念に縛られすぎてきた, ステレオビジョンを考え直して見ようというのがこのHPの趣旨ですが, しばらくはこれまでの考え方を説明します.

対応を見つけるために,特徴のある点(特徴点,コーナーと呼ばれる)の特徴量の類似性や配置の対応検出, 矩形領域の画像の類似性などが使われてきました. これらは

  • 特徴ベースのステレオマッチング
  • 領域ベースのステレオマッチング
  • ピクセルベースのステレオマッチング
と分類されていました.ピクセルベースはその他の意味です. ところが,最近ではすべてのピクセルに対して視差を検出しようとする 稠密ステレオ(Dense Stereo)が一般的になってきて, ピクセルベースが当たり前になって来たようです. それでこの分類はもう使われなくなっています.

稠密ステレオでは,右画像の全てのピクセルに視差を割り当てる問題を解くことになります. 視差を割り当てるのがなぜ右画像なのかというと,右方向を正とすることが自然なので,視差を正とすると, 右画像のあるピクセルに投影された空間上の点は左画像では少し右のピクセルに投影されるのですから, 左画像でどれだけ右に現れるかを視差と呼んで, 「右画像のピクセルに右方向の視差,すなわち正の視差を割り当てる」のが自然というわけです.

右画像の全てのピクセルに視差を割り当てる様な問題はサイトラベリング問題と呼ばれます. サイトは場所のことで場所に何らかのラベルを割り当てる問題というわけです. 今の場合はピクセルが場所で,視差がラベルになります. 割り当ての数はラベル数のサイト数乗 ( assign_num = label_num site_num ) という膨大な数となります. この中から1つの割り当てを選ぶ訳ですが,普通は割り当てに対して単一の数値を与える関数を定義して, その数値が最小となる割り当てを選択します. たった一つの数値で良いか悪いかを判断して良いのかと感じますが, これは物理学からのアナロジーで, エネルギーという単一の実数値を与える関数(その様な関数をハミルトニアンという)から その系のすべてが分かるという定式化を真似たものです. ですので,割り当てに与えられたこの単一の数値をエネルギー, エネルギーを与える関数をエネルギー関数, そして最小のものを選ぶことをエネルギー最小化と呼んでいます. そして割り当てを,物理的状態を表す用語である,配置と呼ぶこともあります.

サイトラベリング問題のエネルギー関数は f を割り当てとして E ( f ) = Edata ( f ) + Esmooth ( f ) と表現されます.簡単のためサイトはピクセル,ラベルは視差として説明します. Edata ( f ) は画像データからピクセル毎に視差の割り当てが適切かどうかを決める項です. ピクセルだけを見て対応がはっきりと分かるようなステレオ画像であれば, この項だけで十分なのですが,実際は同じような明るさと色のピクセルばかりですから, ピクセル間の関係を見る第2の項 Esmooth ( f ) が必要になります.この第2項をどのようにするかは非常に重要な論点で様々なものが提案されています. まず隣接する2つのピクセルを問題にすることが最初の第一歩となります. たとえば視差は隣接するピクセル間で大抵はそんなに大きく変化しないだろうというものです. すなわち第2項はsmooth性に関するものだということで Esmooth という名前が使われます.

2つのサイトの関係を問題にする場合を1階と言い, 3つ以上のサイトの関係を問題にする場合を高階と言います. 1階であってもどのような関数にするかは大問題です. 特に視差は,大部分でスムーズに変化するものの, 対象の境界部分では大きく変化するのが当然ですから, 単にスムーズにするだけではいけない, 厄介なものだと考えられてきました. それで少しの変化では変化に応じてエネルギーを高くするものの, ある程度を超えるとそれ以上はエネルギーを大きくしない, 打ち切り比例 ( truncated linear ) と呼ばれる関数が良く使われています. なおある配置のエネルギーを高くするということはそのような配置を起こりにくくすることに対応します. 一方,gaze_line-depth modelでのdepthは, 視差のように対象の境界で大きく変化することがありませんので, 単純比例 ( linear ) の関数を使うことができます. 実は禁止条件があるのでdepthは隣接間で最大1しか変化しません. 禁止条件も入れて凸の関数にすることができます.

またエネルギーを最小にする配置を簡単に見つけられるかどうかはエネルギー関数の形に強く依存します. 1階のエネルギー関数が凸と言われる性質を持つときには, グラフカットを使って一気に最小配置を求めることができます. 打ち切り比例 ( truncated linear ) は凸ではないので, alpha expansionやalpha-beta swapなどの近似アルゴリズムが使われます.

この記事で紹介している最初のdisparity_stereo.cppは, pixel-disparity modelをlinearなエネルギー関数を使って一気に最小配置を求めています. このlinearなエネルギー関数とはサイト対に対して Esmooth ( f ) = penalty * | i - j | と表現されるものです.ここで ij は隣接するpixelのdisparity値です.

2番目のdepth_stereo.cppと3番目の3d_stereo.cppは, gaze_line depth modelを凸なエネルギー関数を使って一気に最小配置を求めています. この凸なエネルギー関数とはサイト対に対して Esmooth ( f ) = penalty * | i - j | + inhibit * ( | i - j | - 1 ) * T ( | i - j | > 1) と表現されるものです.ここで ij は隣接するgaze_lineのdepth値です. T ( bool ) はboolが真のとき1,そうでないとき0となるデルタ関数です. この式を見ると分かるように,inhibitを大きくすると隣接するgaze_lineのdepth値は最大1しか変化しません. このinhibitが記事2で指摘したgaze_line-depth modelの重要な性質を実現している訳です. これにより大部分でスムーズに変化するものの, 対象の境界部分では大きく変化するという厄介な視差の性質を考慮する必要がなくなり, グラフカットによる完全な最小化が可能な凸のエネルギー関数を使えるようになったのです.