Chapter1 4.6

CV勉強会発表資料
レベルセット法
4.6 – 4.12
Presented by takmin

概要
• 4.6 局所成長速度場と拡張成長速度場
• 4.7 Narrow Band
• 4.8 再初期化
• 4.9 差分法を用いたLevel Set Method
• 4.10 Geodesic Active Contour
• 4.11 Fast Marching Method
• 4.12 領域の統合を用いたLevel Set Method

4.6 局所成長速度場と拡張成長速度
場
4.4章のLevel Set Methodでは，成長速度Fijをグ
リッド（ピクセル）ごとに補助関数φij及び画素値
Iijから毎回算出していた．

局所成長速度場

Fi , j  k I ,i , j (a  b i , j ) (20)
濃度勾配が補助関数φ
大きいと小さの曲率
くなる

全てのピクセルごとにFを計算する
F=2 F=3 F=2 F=1 F=2 F=1 F=3 F=1 F=3 F=2
F=2 F=1 F=1 F=1

Φ=0

全てのピクセルごとにFを計算する
F=2 F=3 F=2 F=1 F=2 F=1 F=3 F=1 F=3 F=2
F=2 F=1 F=1 F=1 F=3 F=1 F=3 F=3 F=1 F=3

F=3 F=2 F=3 F=2 F=1 F=1 F=3 F=2 F=1 F=2 Φ=0

F=1 F=3 F=2 F=3 F=2 F=1 F=3 F=2 F=1 F=2
F=2 F=3 F=2 F=1 F=2 F=1 F=3 F=1 F=2 F=3

F=2 F=3 F=2 F=3 F=2 F=2 F=2 F=3 F=3 F=2

F=2 F=3 F=3 F=1 F=1 F=1 F=3 F=2 F=2 F=1

F=3 F=2 F=1 F=3 F=2 F=2 F=1 F=2 F=1 F=2

F=2 F=1 F=3 F=1 F=2 F=2 F=3 F=2 F=2 F=3
F=3 F=2 F=3 F=1 F=2 F=2 F=3 F=2 F=2 F=4

拡散成長速度場
Φ=0の領域でFを計算

F=1 F=2

F=1 F=3 F=2
F=2 F=1 F=2

F=3

F=2

注目画素の最近傍zero level探索

F=1 F=2

F=1 F=3 F=2
F=2 F=1 F=2

F=3

F=2

成長速度をコピー
F=3

F=1 F=2

F=1 F=3 F=2
F=2 F=1 F=2

F=3

F=2

最近傍zero level setの成長速度をコピー
F=3 F=3

F=1 F=2

F=1 F=3 F=2
F=2 F=1 F=2

F=3

F=2

最近傍zero level setの成功速度をコピー
F=3 F=3 F=2 F=1 F=1 F=1 F=3 F=1 F=1 F=2

F=1 F=2

F=1 F=3 F=2
F=2 F=1 F=2

F=3

F=2


F=3 F=3 F=2 F=1 F=1 F=1 F=3 F=1 F=1 F=2
F=3 F=3 F=2 F=1 F=1 F=1 F=3 F=1 F=1 F=2

F=3 F=3 F=2 F=1 F=1 F=1 F=3 F=2 F=1 F=2

F=3 F=2 F=2 F=1 F=2 F=1 F=3 F=2 F=1 F=2
F=2 F=3 F=2 F=1 F=2 F=1 F=3 F=2 F=2 F=2

F=2 F=3 F=2 F=1 F=2 F=2 F=3 F=2 F=2 F=2

F=2 F=3 F=3 F=1 F=2 F=2 F=3 F=2 F=2 F=2

F=2 F=2 F=3 F=1 F=2 F=2 F=3 F=2 F=2 F=2

F=2 F=2 F=3 F=1 F=2 F=2 F=3 F=2 F=2 F=2
F=2 F=2 F=3 F=1 F=2 F=2 F=3 F=2 F=2 F=2

境界進行の比較
局所成長速度

拡張成長速度

4.7 Narrow Band
補助関数を画像全体に対して計算するのは無駄

ゼロ等高面の近傍のみ補助関数を更新すれば良い

φ
zero level
set
Narrow Band

δ

4.7 Narrow Band
ところで、、、、、、

図19って正しくは↓じゃないですか？？？？

φt

Step 1: 変数の準備
変数φi,j，Fi,j，Si,j ，Pk，n=0を準備する。
処理領域（narrow band）内か
を表すフラグ

φi-1,j-1 Fi-1,j-1 φi,j-1 Fi,j-1 φi+1,j-1 Fi+1,j-1
Si-1,j-1 Si,j-1 Si+1,j-1 φ=0
φi-1,j Fi-1,j φi,j Fi,j φi+1,j Fi+1,j
Si-1,j Si,j Si+1,j
φi-1,j+1 Fi-1,j+1 φi,j+1 Fi,j+1 φi+1,j+1 Fi+1,j+1
Si-1,j+1 Si,j+1 Si+1,j+1

Pk-1
Pk = (px,k, py,k)

Pk+1

Step 2: 符号付き距離場の構築
• 初期閉曲線C(p,0)=C0(p)を設定
• 補助関数値φijを初期閉曲線からの距離を元に設定
• 閉曲線の周辺に幅δのNarrow Bandを設定
• Narrow Band内のSi,jを1，それ以外のSi,jを0に設定

φ
C0(p)
φ=0 Si,j=1
Narrow Band

δ

Step 3: 成長速度の計算
• nに1を加える。
• Si,j=1の領域のみ成長速度Fi,jを計算

φ=0 Fi , j  k I ,i , j (a  b i , j )

Step 4: 補助関数の更新
• Si,j=1の領域のみ補助関数値φi,jを更新

i , j  i , j  tFi,j i , j
φ=0

φ

Step 5: Zero level setの検出
• Si,j=1の領域においてφ=0となる位置を検出し，閉
曲線C(p,t)を求める

C(p,t-1)

C(p,t-1)

Step 6: 再初期化
• nが適当な間隔，または閉曲線がNarrow Bandの端に近づいた場
合，再初期化（後述）を行い，閉曲線からの符号付き距離をφi,jに設
定する．
• 現在の閉曲線からの距離がδ以内にある領域のSi,jを1に，それ以
外を0にセットする．

δ

C(p,t-1)

Step 7: 繰り返し
• Step 3～Step 6をφi,jの変化量が閾値以下にな
るか，nがあらかじめ決められた繰り返し回数
を超えるまで繰り返す．

4.8 再初期化
Upwind Schemeで補助関数を逐次的に更新

積分誤差も積算される

定期的に現在のzero levelからの距離を元にφ全体を更新
再初期化

φ φ

再初期化のイメージ

4.8 再初期化
• 全てのグリッドで近傍zero level setを探索し，
距離を計算するのはコスト大
– Adalsteisonの方法
• First Marching Method（後述）
– 倉爪らの方法
• 再初期化と拡張成長速度場の計算を統一

4.8 再初期化
• Sussmanの方法
– Level Setの更新の間に以下の式を繰り返し解く
d
m  sgn(m )(1  m )
0
(37)
d

m  1 で収束傾きが一定
になる φ

4.8 再初期化
• Sussmanの方法（僕なりの解釈）
– Step 4 補助関数φの更新
i , j  i , j  tFi,j i , j
– 以下の式でφを更新
i , j  i , j  sgn(i0 j )(1  i , j )
,

傾きが一定
i , j  1 で収束になる φ

4.9 差分法を用いたLevel Set Method
• 補助関数を以下の式で更新していく方法

Φ t 1 1
A Φ t
(Φ t 1
 { })
t 1
i, j

行列Aを求めたい

差分法の導出
φの更新式を以下のように定義する
   
   (38)     
   (39)
t 曲率  

「更新の大きさ」＝「φの曲率」×「φの傾きの大きさ」

曲率 φ
傾き

φi,j

差分法の導出
(38) ,(39)式を離散化する
   
   (38)     
   (39)
t  

 
n n 1  ( xn ,  yn ) 
    n 1 (40)
t  
n 1

 

これを最終的に，  f ( ) の形式に変形したい
n n 1

差分法の導出
 n   n 1  ( xn ,  yn ) 
    n 1
t  
n 1

 

力技で式を展開（省略）

in j  in 1
,

t
,j
  Ai  k , j l (in k , j l  in j )
k {0 , 1},l {1, 0}
 , (43)

1 2 in 1
,j
Ai  k , j l  2 (45)
h in 1, j l  in 1
k ,j

差分法の導出

  n
1 
 {1,0}  k {0,{1,0}
Ai  k , j l t i , j  Ai  k , j l tin k , j l in1
 ,j (49)
 k {0, 1},l  1},l

新しい補助関数値を用いた式古い補助
関数値

i, j を一列に並べたベクトル u   i, j 
n
n

n 1
Αu  u n

1 n1
u Α u
n

差分法の導出
n 1
Αu  u
n

  1n1   1n11 
  ,   , 
 0      
  n   n 1 
  1, H   1, H 
          
 A t  A t 1  A Ai , j 1t  Ai 1t  in j    in 1 
 i 1, j i , j 1  i k , j l t  ,   , j 
 k ,l

          
 n   n 1 
  W ,1   W ,1 
      
 0    n 
  W , H   W1 
n
    ,H 
行列Aがとてもデカイ！
要素数： W  H  2

差分法の流れ
1. 変数φi,jと変数Ai,jと繰り返し回数nをグリッド毎に
準備
2. n=0として初期閉曲線C(p,0)を設定し，Cからの
距離を元に符号付き補助関数値φi,jを求める
3. nに1を追加し，φからAi,jを計算する
4. 以下の式からφを更新
1 n1
u Α u
n

5. nの適当な間隔で，再初期化処理を行う
6. 3～5を収束するか指定回数を超えるまで繰り
返す

差分法の導出
曲率だけでなく，画像の濃度勾配を含めた場合

 k I (a  b )  (46)
t 濃度勾配曲率
の逆数

in j  in 1
,

t
,j
  Ai  k , j l (in k , j l  in j )  k I a in 1
k {0 , 1},l {1, 0}
 , ,j
(48)

kI b 2 in 1
,j
Ai  k , j l  2 (47)
h in 1, j l  in 1
k ,j

差分法の導出

  n
1 
 {1,0}  k {0,{1,0}
Ai  k , j l t i , j  Ai  k , j l tin k , j l  k I a in1 t  in1
 ,j ,j
 k {0, 1},l  1},l

i, j を一列に並べたベクトル u   i, j 
n
n

n1
Αu  B  u
n

u Αn 1
u n 1
B 

4.9.1 ADI法
•差分法は行列Aが巨大なため，計算量が大きい
•行列Aの大部分の要素がゼロ

近似解法で計算の高速化・安定化を行う

＜ADI法＞
x方向の解とy方向の解を逐次的に計算する

4.9.1 ADI法
x方向の更新
in j  in 1
 Ai 1, j l (in 1, j  in j )  Ai 1, j (in 1, j  in j )
, ,j
 ,  ,
t
通常の差分法  Ai , j 1 (in j 1  in j )  Ai , j 1 (in j 1  in j )
, , , ,

 A i  k , j l
k {0 , 1},l {1, 0}
(in k , j l  in j )
 ,
(43)

ADI法 i*, j  in 1
 Ai 1, j l (i*1, j  i*, j )  Ai 1, j (i*1, j  i*, j )
,j
（x方向） t
 Ai , j 1 (in 1  in 1 )  Ai , j 1 (in 1  in 1 )
,j
1
,j ,j
1
,j
(44)
y方向は前時刻n-1の値を利用する

4.9.1 ADI法
i*, j  in 1
,j

t (52)
 Ai , j 1 (in 1  in 1 )  Ai , j 1 (in 1  in 1 )
,j
1
,j ,j
1
,j

 i*1, j 
 * 
a
i 1, j ai , j ai 1, j  i , j   bi , j (53)
 * 
 i 1, j 
更新する補助関数値
ai 1, j   Ai 1, j t
（x方向のみ）
ai , j  1  Ai 1, j t  Ai 1, j t
ai 1, j   Ai 1, j t

j ,j ,j  ,j  ,j
1

bin, 1  in 1  Ai 1, j t in 1  in 1  Ai , j 1t in 1  in 1
1
,j 

4.9.1 ADI法
 *
i, j を一列に並べたベクトル u
*
 
  *
i, j
n 1
A xu  b*

   0
 
Ax   ai 1, j ai , j ai 1, j 
0   
 

1 n 1
u Α b
*
x

まずはx方向のみ計算！

4.9.1 ADI法
僕の解釈
 *
i, j を一列に並べたベクトル u*  i*, j   このテキストの表現が間違って
いる？
n 1
A xu  b
* 本当はこういうこと？

   0  1*,1  1*,H   b1,1  b1, H 
    
  i , j      bij
*
 ai 1, j ai , j ai 1, j  
0
    W ,1  W , H   bW ,1  bW , H 

* *
  
2
行列Axの要素数= W 計算列数= H
行列Ayの要素数= H 2 計算列数= W
計： W 2  H  H 2 W W  H 2 よりは尐ない

4.9.1 ADI法
y方向の更新
in j  in 1
 Ai 1, j l (in 1, j  in j )  Ai 1, j (in 1, j  in j )
, ,j
 ,  ,
t
通常の差分法  Ai , j 1 (in j 1  in j )  Ai , j 1 (in j 1  in j )
, , , ,

 A i  k , j l
k {0 , 1},l {1, 0}
(in k , j l  in j )
 ,
(43)

ADI法 in j  in 1
, ,j
(y方向) t
 Ai , j 1 (in j 1  in j )  Ai , j 1 (in j 1  in j )
, , , ,
(44)
x方向は先程の逐次更新の値を利用する

4.9.1 ADI法
 n
i, j を一列に並べたベクトル u
n
 
 n
i, j

A yu  b
n *

 * を使って計算
1 *
u Α b
n
y

更新完了！

4.9.2 AOS法
＜ADI法＞
x方向の計算を行ってから，その結果を元にy方
向を計算

＜ASO法＞
x方向とy方向を逐次的に別々に計算して，平均
を取る

4.9.2 AOS法
x方向
i*, j  in 1
,j

t
A xu*  u n1
 Ai , j 1 (in 1  in 1 )  Ai , j 1 (in 1  in 1 )
,j
1
,j ,j
1
,j

y方向
i*,*j  in 1
 Ai 1, j l (in 1j  in 1 )  Ai 1, j (in 1j  in 1 )
,j
1, ,j 1, ,j
t A y u**  u n 1
 Ai , j 1 (i*,*j 1  i*,*j )  Ai , j 1 (i*,*j 1  i*,*j )

u*  u**
un 
2

4.10 Geodesic Active Contour
以下の式で定義された成長速度を持つ曲線を
Geometric Active Contourと言う

Ct  F ( )N (11) Fi , j  k I ,i , j (a  b i , j ) (20)
閉曲線の速度場法線ベク
濃度勾配曲率
時間変化トル
が大きい
ほど小

Ct  N の時，以下のLを最小にする

L   C (q) dq   (q) dq ds
'

Geometric Active Contour
Ct  N の時，以下のLを最小にする
L   C (q) dq   ds
'

輪郭線上速閉曲線の
度場の合計長さ C(p,t-1)

曲率の合計

輝度勾配について考えない場
合，閉曲線の長さを最小にす C(p,t)
る方向に動く

Geometric Active Contourの問題点：
•対象によっては曲線の成長が不安定
•ノイズやギャップ（切れ）に敏感すぎる

Geodesic Active Contour
LR   g C (q)C ' (q) dq   g C ( s)ds (67)

境界の特徴

このLRを最小にするアプローチ

LR   g C (q)C ' (q) dq   g C ( s)ds (67)

この最小化を最急降下法で行うと（証明は参考文献にて）

Ct  g ( I )N  g ( I )  NN
 g ( I )  g ( I )  NN
  F ' ( )N
と解ける．
ちなみに， F ' ( ) は F ( ) の微分という意味ではない．（お陰で長時間悩んだ）

補助関数の更新式は，以下の通り

t   F ' ( ) 
 g ( I )  g ( I )  N 
 g ( I )   g ( I )  
境界特徴境界特徴補助関数
の勾配の勾配
補助関数
補助関数
の勾配の
の曲率
大きさ

まだ直感的にわからない・・・。

Geometric Active ContourとGeodesic Active Contourを比較して，
理解を深める．

ここで以下の式の，a=0，b=1とすると

Fi , j  k I ,i , j (a  b i , j ) Fi , j  kI ,i , j i , j

g(I)は境界に関する関数なので，kと対応させる．
1
g ( I )  k I   
1  G  I
濃度勾配の大きさ
F  g (I )

Geometric Active Contour Geodesic Active Contour

t   F ( )  t   F ( ) 
'

 k I    g ( I )   g ( I )  
 g ( I ) 
Geometric
との違い

濃度勾配の空間的な変化
を見ている

ボケやノイズ，ギャップに
強くなる

Gradient vector fieldを用いたGradient vector flow：
• エッジの位置情報だけでは動的輪郭の移動方向が正しく判別出来
ない
• あらかじめエッジの大まかな位置を表すベクトル場を作成しておき，
Geodesic Active Contourをこのベクトル場に基づいて移動させる．

4.11 Fast Marching Method
閉曲線が，常に一方向にしか動かないと仮定するこ
とで高速化を図る方法

常にF>0の例

常にF>0で，一次元の例:

x  F  T (x) dT ( x)
F 1
移動距離成長到達時間 dx
速度

多次元へ一般化

Eikonal方程式

T ( p ) F  1

Eikonal方程式 T ( p ) F  1
2 2
 T   T 
T  
 y    y 
  
   

maxD 
1
x
i, j i, j ,  2
, 
T , D T ,0  max Di jyT , Di jyT ,0
x 2
 1 / Fi , j
または

maxD 
1
x
i, j 
2
  2

T ,0  min Di jxT ,0  max Di jyT ,0  min Di jyT ,0
, , ,   2
 1 / Fi , j

x
Ti 1, j  Ti , j y
Ti , j 1  Ti , j
D T
i, j D Ti, j
h h
Ti , j  Ti 1, j Ti , j  Ti , j 1
Di xT 
,j Di jyT 
,
h h

Step 1: 初期化，Heapの構築
初期境界に属するグリッドをknownのリストに追加し，到達
時間Tを0にセットする

Heap形式のリスト

known(T=0)

knownの四近傍グリッドのうち，knownでないグリッドをtrial
とする．

trial Heap形式のリスト

known(T=0)

trialグリッドの時間をTi,j=1/Fi,jから求め，リストに登録する．

(4,2) T=1
known(T=0) (5,3) T=2
(4,4) T=4
(3,3) T=3

これらのグリッドを到達時間に関する昇順でソート

(4,2) T=1
known(T=0) (5,3) T=2
(4,4) T=4
(3,3) T=3

これらのグリッドを到達時間に関する昇順でソート

(4,2) T=1
known(T=0) (5,3) T=2
(3,3) T=3
(4,4) T=4

それ以外のグリッドをfarとし，到達時間を無限大とする

far (T=∞)

(4,2) T=1
known(T=0) (5,3) T=2
(3,3) T=3
(4,4) T=4

Step 2: Heap先頭の取り出し
Heapの先頭に置かれた，trialリストの中で一番到達時間T
が短いグリッドを選択し，リストから除外して，knownとする
far (T=∞)

(4,2) T=1
known(T=0) (5,3) T=2
(3,3) T=3
(4,4) T=4


known(T=1) Heap形式のリスト

known(T=0) (5,3) T=2
(3,3) T=3
(4,4) T=4

Heapから除外する際，同時にHeapを再構築する．

(5,3) T=2
known(T=0) (3,3) T=3
(4,4) T=4

Step 3: 取り出されたグリッド近傍の選
択
現在選択されたグリッドの４近傍のうち，farに属しているグ
リッドをトライアルのリストに追加する

(5,3) T=2
known(T=0) (3,3) T=3
(4,4) T=4

Step 3: 取り出されたグリッド近傍の選
択
現在選択されたグリッドの４近傍のうち，farに属しているグ
リッドをtrialにする

(5,3) T=2
known(T=0) (3,3) T=3
(4,4) T=4

Step 4: Heapへの追加と再構築
現在選択されたグリッドの４近傍のうち，trialに属しているグ
リッドの到達時間Tを計算して仮登録し，Heapを再構築する．

(5,3) T=2
known(T=0) (3,3) T=3
(4,4) T=4
(4,1) T=2
(5,2) T=1.5
(3,2) T=4

Step 4: Heapへの追加と再構築
現在選択されたグリッドの４近傍のうち，trialに属しているグ
リッドの到達時間Tを計算して仮登録し，Heapを再構築する．

(5,2) T=1.5
known(T=0) (5,3) T=2
(4,1) T=2
(3,3) T=3
(4,4) T=4
(3,2) T=4

Step 5: 繰り返し/終了処理
trialに属しているグリッドが存在すればStep 2へ戻る．それ
以外の時は処理を終了する．

(5,2) T=1.5
known(T=0) (5,3) T=2
(4,1) T=2
(3,3) T=3
(4,4) T=4
(3,2) T=4

4.12 領域の競合を用いたLevel Set
Method

濃度勾配が明確でない場合

領域の競合を用いる

Method
領域分割のためのエネルギー関数：
C
E (C , u )    ( I  u ) dxdy   
2
ds    u dxdy
 s  C

Mumford-Shahの枠組み問題と呼ばれる

なんのことかわからん！

Method
分り易く画像全体を２つの領域 ΩinとΩout に分ける場合を考
える
E (C , in ,  out ) 境界線の長さ

C
   ( I  in ) dxdy    ( I   out ) dxdy   
2 2
ds
 in  in s
領域Ω inでの輝度値の
バラツキ

I: 濃淡画像の輝度値
μ in ,μ out: 領域内外の輝度値の平均値 Ωin

このエネルギー関数の最小化！
Ωout

Method
E (C ,  in ,  out ) をストークスの定理を用いて時間tで偏微分

 
Ct   ( I  in ) 2  ( I  out ) 2 N   k N (82)

このCtでエネルギー関数が最小

Ct  F ( )N t   F ( ) 

なので，

t
 
  ( I  in ) 2  ( I  out ) 2    k  (83)

Method

  ( I   )  ( I   )    
in
2
out
2
k (83)
t

レベルセットの補助関数更新時に，この領域分割の項を足
す？（テキストに記述なし）

Chapter1 4.6

Recommended

Recommended

More Related Content

What's hot

What's hot (20)

Viewers also liked

Viewers also liked (8)

More from Takuya Minagawa

More from Takuya Minagawa (20)

Chapter1 4.6