A LOOPY BELIEF PROPAGATION APPROACH TO THE SHAPE

A LO O P Y BELI EF P R O PAGATI ON A P P R O A CH TO
TH E S H A P E F R O M S H A D IN G P R O BLEM
Markus Louw, Fred Nicolls, Dee Bradshaw
In Proc. 2nd International Conference on Computer Vision
Theory and Applications, 2007
陰影からの形状復元へのループ確率伝搬アプローチ
紹介者:2DS09049E 小川真人 原研究室
1. 導入
 陰影からの形状復元(SFS: Shape From Shading )
 ループ確率伝搬法(LBP: Loopy Belief Propagation )
 論文概要
 Shape From Shading 陰影からの形状復元(SFS)
濃淡画像から立体形状を
生成する技術
画像例(http://perception.inrialpes.fr/)
濃淡画像
反射モデ
ルの仮定
既知の光
源方向
 Loopy Belief Propagation ループ確率伝搬法(LBP)
閉路のあるMRFの近似推論
有向グラフ
無向グラフ
ベイジアンネットワーク
(Bayesian network)
変数の依存関
係のグラフ化
因子グラフ
マルコフ確率場
(MRF : Markov Random Field)
メッセージ
伝搬
同時確率の最
大値を求める
 論文概要
• SFS問題のLBPを用いた計算を試みる
• 複雑な関数を用い、平滑性を維持する
• 多重解像度補間を用いる
2. ランバートライティングモデル
物体表面の拡散反射モデルを仮定

Rn( x)  I ( x)
反射関数=法線の関数

w
位置 x の画像強度

n (x)
θ
a
ab
(1)
a
b
b
明るい
ab
より明るい
3. ループ確率伝搬法(LBP)
確率の基本法則とグラフィカルモデル
加法定理 p( X )   p ( X , Y )
ベイズの定理
Y
乗法定理 p ( X , Y )  p (Y X ) p( X )
対称性 p(Y X ) 
p( X , Y )  p(Y , X )
p( X Y ) p(Y )
p( X )
b
a
a
c
p(a, b, c)  p(c a, b) p(b a) p(a)
b
c
p(a, b, c)  p(a) p(b a) p(c b)
グラフィカルモデル
• 「全確率変数上の同時分布が、一部の変数のみに
依存する因子の積としてどのように分解可能か」と
いう情報を表現する
a
b
c
p(a, b, c)  p(a) p(b a) p(c b)
p(a, c b) 

p(a, b, c)
p (b)
p(a) p(b a) p(c b)
p(b)
 p ( a b ) p (c b )
つまり、このグラフは変数間の条件
付き独立性を表現している
bが与えられた時、同時分布が積に分解される
 グラフ表現の色々
有向グラフ
c
b
a
p(a, b, c)  p(a) p(b a) p(c b)
無向グラフ
因子グラフ
fa
x1
x2
x3
x1
x2
x3
fb
fc
p( x1 , x2 , x3 ) 
1
 1, 2 ( x1 , x2 ) 2,3 ( x2 , x3 )
Z
p( x1 , x2 , x3 )  f a ( x1, x2 ) fb ( x2 , x3 )
p( x1 , x2 , x3 )  f a ( x1 , x2 ) fb ( x2 , x3 ) f c ( x2 , x3 )
確率伝搬法:グラフィカルモデル上の確率推論の計算法
周辺確率は xi を除くすべての変数に対して同時分布の和をとることで得られる
p( xi )   p( x)
x \ xi
x1
x2
xi 1
x3

xi
xK 1
xK

各ノードのとりうる状態の数をNとするなら、xの状態の数は N k であり計算量やメモリは
連鎖の長さ N に対して指数的に増加する
上のグラフで表わされる分布は関数の積に因数分解されることを見てきた
p( x ) 
1
 1, 2 ( x1 , x2 ) 2,3 ( x2 , x3 )  K 1, K ( xK 1 , xK )
Z
この条件付き独立性を利用すればはるかに効率的に計算できる!
「メッセージ」という量を導入し、その積で周辺分布を表す
p( x ) 
x1
x2
1
 1, 2 ( x1 , x2 ) 2,3 ( x2 , x3 )  K 1, K ( xK 1 , xK )
Z

p( xi )  
x \ xi
1
p( x ) 
Z
i 1i (ai )
i 1i (ai )
x3
xi 1
xi

xK 1
xK




 i 1,i ( xi 1 , xi )   2,3 ( x2 , x3 )  1, 2 ( x1 , x2 )  
 xi1
 x2
 x1
  


 

(
x
,
x
)


(
x
,
x
)
 i ,i 1 i i 1
 K 1, K K 1 K  
 xi1
 xK
 
i 1i ( xi )
メッセージ
i 1i ( xi )
1
p( xi )  i 1i ( xi ) i 1i ( xi )
Z
メッセージ
「メッセージ」はノードごとの逐次計算で得られる
N1 N 2
Ni
z1
zi
i 1
i 1i (ai )    a , z  j , j 1 ( z j , z j 1 )
i
z2
i
j 1
N1 N 2
12 (a2 )   a , z  1, 2 ( z1 , z2 )
2
z1
2
z2
N 2 N3
23 (a3 )   a , z 12 ( z2 ) 2,3 ( z2 , z3 )
3
z2
3
z3
Ni1 Ni
i 1i (ai )   a , z i 2i 1 ( zi 1 ) i 1,i ( zi 1 , zi )
i
zi1
i
zi
この方法をとることで、計算量はN K からNK2になる
もうちょっと直観的な解説
x1

xi 1
1
2
a2
N1

z1
i 1i (ai )
i 1i (ai )
x3
x2
z1
N
x1 の値に関係なく x2 がある状態を
とるときの確率を計算する
この関数がメッセージ
これを各々伝搬させていく
xi

xK 1
xK
因子グラフの場合のメッセージ伝搬
メッセージが2種類になるが本質的な違いは無い
 gx (x)
g
f
x
変数ノードから因子ノード
x f ( x)   ( x) g x ( x)
 x f (x)
(6)
g f
 y f ( y)
u\ x
y
x
f
 f x (x)
因子ノードから変数ノード
 f x ( x)   f (u)  y f ( y)
u\ x
y x
(7)
• 繰り返し後、各コーナー頂点ノードの高さの
事後分布は次式で与えられる
px ( x)   ( x)  g x ( x)
(8)
g
 gx (x)
g
x
グラフに閉路(ループ)がある場合、メッセージが収束するまで
伝搬を繰り返すことで近似する。(収束は保障されていない)
同時確率を最大化する場合
x
max
 arg max p( x )
x
x1
x2
xi 1
x3

xi
xK 1
xK

同時分布の最大値も効率的に計算することができる
1
max p( x )  max  max  1, 2 ( x1 , x2 )  N 1, N ( xN 1 , x N )
x
xN
Z x1
1






 maxmax 1, 2 ( x1 , x2 )  max N 1, N ( xN 1 , xN )  

  
xN
Z x1  x2 
このときの x の値を推定結果とする
4. SFSに対するLBPの定式化
• ピクセルの四隅に高さを表す変数ノードを置く
点(ノード)の高さがわかれば立体を復元できる
高さの確率をモデル化し条件下で最も尤もらしい値を推定
• (X+1)(Y+1)個の確率変数の同時分布と考える
A11
A12
X
N
0
Y
Pr{A  a}  p(a)
高さの階調数を N とすると
a が取りうる場合の数は N (X+1)(Y+1)
A  ( A11 , A12 ,, A( X 1)(Y 1) )T
a  (a11 , a12 ,, a( X 1)(Y 1) )T
どのような依存関係をモデル化するか?
• 変数ノードの作る三角形に対してエネルギー
関数(因子)を与える
L
θ
n
158/255
各頂点3対は一意の平面を形成し、平面の光源方向に対する角度は確率をその3
点の配置に割り当てることを可能とする。
L
これらのグラフを数式で表わすと…
θ
1. 頂点ノードX
2. 画像データY
3. 既知の距離データZ
xi
これらの条件付き確率は次のようになる
xk
PX Y , Z  
n
xj
158/255
t
exp



ijk ( xi , x j , xk , yijk ) exp i ( xi , zi )
i , j ,k:k  j i
i


 ( xi , x j , xk , yijk )  yijk  nijk ( xi , x j , xk )  L
t
ijk
頂点ノードi,j,kを囲む
画像領域の平均強度
(9)
(10)
照度方程式(式5)の
右辺
静止/動光源を情報を含む場合、関数を変形させる


 ( xi , x j , xk , yijk )  yijk  nijk ( xi , x j , xk )  L
t
ijk
(10)
N



t
 ijk ( xi , x j , xk , yijk )   yijks  nijk ( xi , x j , xk )  Ls
s 1

yijk
Nは画像枚数。 ベクトル
はsで番号づけられた各強度画像の平均強度を

格納していることを示す。各画像でs回繰り返す。 L は画像sでの光源方向
s
(13)
まとめるとこのようなモデルになる
変数ノード:(h+1)(w+1) 個
因子ノード:4wh 個
• 平滑化を合わせた場合は8whになる(後述)
4.1 境界条件と距離データ
• LBPでは制約条件を簡単に扱うことができる
xi ( X  h)  exp h  u( xi ) 
u ( xi ) は既知の距離データ。
(14)
境界条件の設定、点に距離データを与え
た場合のどちらも同じ方法で扱う。
4.2 平滑化
4.2.1 変形サイズ三角形での平滑化
• 高周波誤差を含む解に収束する場合があるの
で、これを軽減するために大きな三角形因子を
追加する。
C
D
三角形面を大きくする関数を加えることで
低周波成分に対する制約が増えるので
結果的に高周波誤差が減る
(関数の形式は式(10)と同じ)
4.2.2 同一線上の3点での平滑化
A
物体表面の滑らかさ
を保つための関数
B
 h( x j )  h( xi ) h( xk )  h( x j ) 
 /  )
 ( xi , x j , xk )  exp(

d ( j, k ) 
 d (i, j )
2
s
ijk
ノード i,j 間の距離
(16)
平滑化パラメータ
ノード間の勾配量が似ている方が尤もらしいとする関数
平滑化項を合わせると条件付き確率は次のように変形される
PX Y , Z  
t
exp



ijk ( xi , x j , xk , yijk ) exp i ( xi , zi )
i , j ,k:k  j i
PX Y , Z  
(9)
i
t
exp



ijk ( xi , x j , xk , yijk ) 
i , j ,k:k  j i
s
 ijk
( xi , x j , xk ) exp i ( xi , zi )
i
(17)
4.3 多重解像
• 高さの近似精度を上げるため、一定の繰り返しの後
に高さ解像度を増やす。
1. for(i=1;i<N;i++)
2.
for(j=1;j<M;j++)
3.
各コーナー頂点ノードに対して現在の高さ推定値を計算
4.
}
px ( x)   ( x)  g x ( x) (8)
g
5.
各コーナー頂点ノードの高さ解像度を上げる
6. }
新しい高さ解像度の計算
1. 新しい高さ解像度を計算:hnew=khold
2. n←1; a←1; b←1
3. while(n<numlabels)
4.
if(LB<Hx+a・hnew<UB)
5.
vertex height(n) ←Hx+a・hnew
6.
n←n+1; a←a+1
7.
}
8.
if(n<numlabels)
9.
if(LB<Hx-b・hnew<UB)
10.
vertex height(n) ←Hx-b・hnew
11.
n←n+1; b←b+1
12.
}
13.
}
14. }
メッセージ
x の状態
(高さ)
0
2
4
6
8
10
補間
0
2
4
6
8
10
徐々に高さの刻み幅を狭くして
推定精度を上げていく
Hxは位置xでの高さのMAP推定値。 hnew は現在のLBP解像度に対する高さ解像度。
k は各解像度サイクルでの前の解像度 hold に適用する圧縮率。LBとUBは下限と上限
解像度上昇に伴うメッセージの更新
頂点ノードの新しい高さを反映するために、頂点ノードからファク
ターノードへのメッセージを更新する。
それにはオリジナル値の点の間の線形補間を用いる。

new
x f

( x)  L 
old
x f
( x), h
new
old
,h

(15)
Lは線形補間で、hnew, holdは新旧の解像度の高さラベル。2点a,c
の間の点bの線形補間は以下で与えられる。
あるいは
f (b)  f (a)  (b  a)  f (c)  f (a)  (c  a)
f (b)  f (c)  (c  b)  f (c)  f (a)  (c  a)
5. 実験結果
実画像と人工画像に対しテストを行った
b
a. 人工画像:異なる滑らかな3D形状から生成
(ランバートモデルでレンダリング)
b. 実画像:ビールの泡(処理前に平滑化)、壺
5.1 人工画像の結果
真値
入力1枚
入力2枚
入力3枚
高さラベル26個
200回の繰り返し中に圧縮率0.8の多重解像7回
泡(20)
5.2 実画像の結果
泡(10)
1,2列目は表面の点を与えない泡
3列目は表面の点を与えない壺
4列目は一つの点の高さを与えた
(同一線上平滑化パラメータσ=1)
壺(5)
点を与えた壺(7)
6. 議論・考察
• 主な問題は局所最小解に落ち込む傾向があ
ること。再構成時にひずみとなる。
• 境界条件や距離データが既知の場合でも観
測画像に対応付けられる表面がいくつも存在
する。
• 克服する方法には、さらに同一線上の3点でのエネル
ギー項や、高さの高周波誤差を除去するための変形
スケールの三角形平滑化の追加が含まれる。
• 頂点ノードの高さのラベル数を増やすことで信頼性を
上げることができるが計算時間を大幅に増加させる。
• 小さい画像(50×50以下)に対してはより信頼性を上
げられる。
• 画像サイズが大きくなると、より局所最小解に陥りや
すくなる。
7. 結論
提案アルゴリズムの特徴
• 頂点の高さ解像度を調整できる
• 他のSFS手法と違い、境界条件と点の距離
データの両方を組み込むことができる。
• 異なる反射モデルを簡単に適用できる
• 深さの増加に伴い、計算時間とメモリ量が増
加する。
ご清聴ありがとうございました