Transcript ヘルムホルツステレオ法
ヘルムホルツステレオ法 CVPR2003サーベイ 宮崎大輔 CVLセミナー 2003年7月2日（水） 本発表の流れ 1. Helmholtzステレオ法とは？ 2. P. Tu, P. R. S. Mendonça, “Surface Reconstruction via Helmholtz Reciprocity with a Single Image Pair,” CVPR2003, Oral 3. T. E. Zickler, P. N. Belhumeur, D. J. Kriegman, “Toward a Stratification of Helmholtz Stereopsis,” CVPR2003, Oral Helmholtzステレオ法とは？ 歴史 • 1856年: von HelmholtzがHelmholtz reciprocityを発表 • 1900年: Rayleigh卿がより一般的な実験結果を提示 • 2002年: Snyderが一般化 Helmholtzステレオ法が開発された • 2001年: Magda, Kriegman, Zickler, Belhumeur [ICCV] • 2002年: Zickler, Belhumeur, Kriegman [ECCV, IJCV] そしてCVPR2003 • Tu, Mendonça: 一対の画像だけを使う手法 • Zickler, Belhumeur, Kriegman: 未校正下で行う手法 ステレオ法 • 多数の方向から物体を観測 • 光源固定 • window内での輝度の差な どを使って対応点を探す • カメラパラメータ（視線方向 など）とdisparity（視差）から 距離を求める • 距離(depth)が求まる 物体 照度差ステレオ法(photometric stereo) • 光源を変化させる→ある方向から光を照らした画像、 別の方向から光を照らした画像、・・・を撮影する • カメラ固定 • 法線(surface normal)が求まる Helmholtz Reciprocity Helmholtzステレオ法 • Stereo + Photometric stereo • 距離(depth)と 法線(normal) が求まる 計測結果の一例 計測装置の一例 取得画像の一例 計測結果の一例 Surface Reconstruction via Helmholtz Reciprocity with a Single Image Pair Peter Tu Paulo R. S. Mendonça GE Global Research Center, USA この論文の手法 • 以前のHelmholtzステレオ法を、一つのペアの reciprocal画像から形状を復元する • 動的計画法(dynamic programming)を使ってマッチ ングを行う – Ohta&Kanadeが1985年に動的計画法を使ったステレオ マッチング法を提案している 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) 画像２のepipolar線上のpixel epipolar 画 像 １ の pixel 線 上 の グリッド 動的計画法(dynamic programming) • コストが最小になる経路を探す →最適なマッチング 動的計画法の実装で有名な物にDijkstra法がある notation • • • • 表面上の点p 単位視線ベクトルv1,v2 単位法線ベクトルn Optical center c1,c2 C1 p v1 v2 n C2 Helmholtz reciprocityと画像の輝度 • BRDF (p, v1, v2 ) • Helmholtz reciprocity (p, v1, v2 ) (p, v2 , v1) • 輝度1 1 I1,2 (p) (p, v1, v 2 )n v 2 • 輝度2 || c2 p ||2 η: scale factor I 2,1 (p) (p, v1, v2 )n v1 1 || c1 p ||2 • 輝度2を輝度1に変換した輝度 2 n v || c p || 2 1 Iˆ1,2 I 2,1 n v1 || c2 p ||2 動的計画法のコスト a b 3D空間の点paが対応づけられている 3D空間の点pbが対応づけられている • Lambertianの場合 aからbにたどり着いたときのコスト = aのコスト + bでのwindow内の各pixelの輝度差 • Helmholtz reciprocityの場合 C(a, b) x1 x0 C( x)dx || pa pb || C( x) (I1,2 (p( x)) Iˆ1,2 (p( x)))2 p( x) xpa (1 x)pb ; 0 x 1 法線のepipolar平面への投影 • • • • 法線を２つの垂直ベクトルの和で表す n no ne neと視線（or光源）v1,v2はepipolar平面に含まれる noはepipolar平面に垂直である o e o e e v n v ( n n ) v n v n v n よって 1 1 1 1 1 v2 n v2 (no ne ) v2 no v2 ne v2 ne n no v2 ne Epipolar平面 v1 法線の近似 nˆ e (pa pb ) (v1 v2 ) || (pa pb ) (v1 v2 ) || pb pa v1 n̂ e v2 α: unknown scale 輝度 • 以上の事を使い、輝度を表す式に代入すると e 2 ˆ n v || c p ( x ) || 2 1 Iˆ1,2 (p( x)) I 2,1 (p( x)) e nˆ v1 || c2 p( x) ||2 • 右辺は、画像２、カメラパラメータ、近似法線から求ま る 飽和(saturation) • カメラのdynamic rangeを超えたらspecularとみなす →入射角と反射角が等しい • よって、コスト関数として以下を用いる C( x) (nˆ e v1 nˆ e v2 )2 結果 結果 Toward a Stratification of Helmholtz Stereopsis Todd E. Zickler Peter N. Belhumeur David J. Kriegman Yale University Columbia University University of California この論文の手法 • 以前のHelmholtzステレオ法を、キャリブレーションを 行わずに形状を復元する • 特異値分解(SVD; singular value decomposition) を使う – Tomasi&Kanadeが1992年にSVDによる因数分解法を提 案している 特異値分解(SVD)とは？ • SVDについては僕は初心者なので間違っている所 があるかもしれないので注意してください • SVDは行列Aを分解する • Uの列、Wの要素、Vの列（VTの行）を入れ替えても いいので、wjを大きい順に並べたとして、以下、議論 を進める • wjを特異値と呼ぶ 特異値分解 • N×Nの行列AのrankがN未満のとき、その行列を特 異という • 行列が特異のとき、Ax=bを解くことを考えよう • このとき、xは解が無数にあるか、解がないかである 特異値分解 • Ax=0を満たすxの張る空間を零空間(null space)と呼ぶ • 零空間が1次元の場合、零ベクトル(null vector)と呼ぶときも ある（普通、零ベクトルとは[0,…,0]のベクトルの事を言う） • Axの張る空間を値域(range)と呼ぶ w1 A= VT U wN 値域の正規直交基底 零空間の正規直交基底 特異値分解 • Ax=bの特殊解はSVDを使って以下で与えられる x V [diag(1/ wj )] (UT b) • Ax=bの解は特殊解に零空間の任意のベクトルを加 えた形で表される • Ax=bの解が無数にある場合、SVDは0に最も近い 特殊解を選ぶ • Ax=bの解がない場合（bが値域の外にある場合）、 SVDは最小２乗法の意味で最良の妥協解を求める。 すなわち、Axがbに最も近くなるような特殊解xをはき だす 特異値分解 Uncalibrated Helmholtz stereopsis • M個の位置にカメラ/光源を 配置する • カメラと光源を入れ替えてペ ア画像を取得 • 合計M(M-1)の画像を取得 • つまり、合計M(M-1)/2の reciprocal画像ペアを取得 Helmholtz reciprocity constraint Τ Τ ˆ s v ˆ s v j j i i eij ˆ 0 e n ji 2 2 | oi p | | o p | j • カメラ/光源位置 o1 oM • • • • • カメラがoi、光源ojのときの観測輝度 表面上の点 p 法線 n̂ 視線（pからoi） v̂i 光源の強さ si ハットは単位ベクトルを表す eij 正射影 • • • • • • • • カメラ/光源が十分遠くにあると仮定（2.5mなど） さきほどの制約： eijsiΤ e jisΤj nˆ 0 光源ベクトル： si sisˆi si は光源の強さで、 ŝ i はカメラ/光源方向（単位ベク トル） T w s 光源と法線の内積をwiとおく i i nˆ 制約は次のように書き直せる eij wi e ji wj 0 T T T T wiを要素とするベクトルw w s1 nˆ s2 nˆ sM nˆ 最終的に制約は次のように表せる Ew 0 具体例 M=4の例 e12 e21 e w e 31 13 1 e14 e41 w2 0 e23 e32 w3 e42 w4 e24 e e 34 43 E w ランク • M=4の場合rankは3になる e12 e21 e w e 31 13 1 e14 e41 w2 w 0 e23 e32 3 e42 w4 e24 e34 e43 • EのrankはM-1になる（らしい） • M-1未満のケースは考えない • １次元零空間は以下で表せる T ck w ck s1Tnˆ sT2 nˆ sTM nˆ , ck R Wの定義 • N×M行列Wを定義する（Nは点の数、Mはカメラ/光 源の数） • Wの行をさきほどの１次元零空間とする c1s1Tnˆ 1 Tˆ c2s1 n2 W T cN s1 nˆ N c1sT2 nˆ 1 c1sTM nˆ 1 c2sT2 nˆ 2 T ˆ cN s M n N • 図はWの列を表す Decomposition • W BSに分解したい • B : N×3行列(ckでスケーリングされた法線） • S : 3×M行列(光源の強さでスケーリングされた光源方向） (rank W = 3) T • SVDで分解 W UV 上位３つの特異値 • ~ ~ 12 B UΣ S Σ1 2VT 3×3の正則行列Qを使って正しい分解を行う ~ ~ 1 B BQ S QS だけを使って分解 • このQを求める必要がある – Uncalibrated stereo, Uncalibrated photometric stereoには様々な方 法があるが、この論文では以下で示す方法を使った – 詳しくは、[Shapiro-Zisserman-Brady1995], [Tomasi-Kanade1992], [Weinshall1993]を参照せよ アフィンカメラ はじめは正射影を仮定 x Pi Xk ti i k image point p P p T 1 T 2 projection matrix scene point 画像i 点k image offset 回転行列の上の２行をスケーリングしたもの （ただし、aspect比を1としピクセル歪みも無いとした場合） 視線/光源方向がPから求まる p1T pT2 sˆ T T | p1 p2 | Metric reconstruction • 先程のはMetric reconstructionではなく、ス ケール倍の曖昧性が残る • そこで、以下を満たす3×3 正則行列Qを求める piT1QQTpi1 piT2QQTpi 2 1 piT1QQTpi 2 0 • 本当の視線/光源方向 求まる ŝ i が • いまいち良く分からないが cos sin sin cos sin 2 cos2 1 sin cos sin cos 0 実装 • 4点以上（例えば20点）の対応点を全ての画像で与 えてやり、Tomasi-Kanadeの方法で視線/光源方向 ŝ i を求める 実装 • 全てのdisparity d∈[0, dmax]で探索し、最適なdを求 める • まず、点xで、各dにおける行列 Ex (d ) を計算する • dが最適であるかを評価する特徴としてrankE=M-1 という性質を利用する • rankを計算する方法は数多く提案されているが、こ の論文では単純に２つの最小の特異値を使って、以 下を評価値とする M 1 rx (d ) M ←ほとんど0 （rankE=M-1が成り立つ場合） 実装 • 点xoの周りのウィンドウWr（例えば9×9）でrx(d)を計算し、その和が最大と なるときのdをdisparityの推定値とする d arg max rx (d ) d xWr • wは wxo arg min || Exo (d )w ||2 , || w || 1 から求まる w • この解は、E=UΣVTとSVDしたとき、最小の特異値に相当するVの列にな る（１次元零空間） • 実際は誤差を防ぐためウィンドウWE（例えば3×3）で重み付きEwを計算 し、その和が最小となるwを求める wxo arg min rx || Exw ||2 , || w || 1 w xWE • この場合でも重み付きEをSVDしたときの１次元零空間としてwを与える事 ができる 実装 • 全てのwxからWを計算できる • 視線/光源方向 ŝ i より、W=BSと分解できる 結果 結果 結果 (c) Daisuke Miyazaki 2003 All rights reserved. http://www.cvl.iis.u-tokyo.ac.jp/