差分法・変分法と平滑化・エッジ

(1)

吉澤信

[email protected], 非常勤講師大妻女子大学社会情報学部

画像情報処理論及び演習II

第6回講義水曜日１限

教室6218

情報デザイン専攻

-フィルタ処理・エッジ強調-

差分法・変分法と平滑化・エッジ

Shin Yoshizawa: [email protected]

今日の授業内容

1. 勾配とエッジの基礎：差分法

.

2.

Laplacianと拡散方程式の基礎:変分法.

3. 演習：エッジ強度抽出と拡散方程式

.

www.riken.jp/brict/Yoshizawa/Lectures/index.html www.riken.jp/brict/Yoshizawa/Lectures/Lec19.pdf

今日の演習は第５回のレポートで出すので、

みなさん頑張ってくださいねー

p(^^)q 第４回レポートは今日〆切です！

明日以降も提出を受け付けます、ただし減点.

後期二回目レポート「Report05」

1. 第５回レポート(後期二回目)：フィルタ処理の内容.

2.

Lec18.pdf, Lec19.pdf, Lec20.pdfに問題のヒントあり.

3. １２月１２日〆切.

www.riken.jp/brict/Yoshizawa/Lectures/index.html www.riken.jp/brict/Yoshizawa/Lectures/Report05.doc

復習：勾配(Gradient)



勾配(Gradient): スカラー場の各点で変化が最大の方向と変化率を大きさに持つベクトル場.



勾配作用素:

).

, ( ) ) , , ( ) , ( (

) , ( ) , (

y x

I y I

y x I x

y x I

y I x y I x I I

 



 



 







) , ( x  y



 



 エッジの大きさ＝

勾配の大きさ：

2 2

y

x

I

I  





勾配ベクトルの表記：

©wikipedia.

©www.mathworks.co.jp.

１階微分は接線、傾き

.

重要：画像のエッジ



画像の勾配: 画像を高さ関数と考えたときの勾配ベクトル場、画像のエッジ部分で大きい勾配ベクトルをもつ画像.

入力I(x,y) エッジ強度画像 I(x,y)

©wikipedia.



勾配ベクトルの方向: 画像

エッジと垂直な方向.

) , (x y

I x方向微分Ixy方向微分Iy

©wikipedia.

) / arctan( I

_y

I

_x

 

Ix

Iy

x

I

 ^{勾配ベクトルの大}_{きさ=エッジ強度:}

2 2

y

x I

I

I  



復習：Laplacian、Laplace-Poisson方程式



ラプラス作用素(Laplacian): 滑らかさを記述.

 Laplace方程式：自然科学の多くの分野で重要.



Poisson方程式: Laplace方程式の右辺が関数.

2 2 2 2 2

y

x 

 



 















 0

 I g I 



2 2 2 2

y I x I I



 



 



©J. Sun et al.

SIGGRAPH 2004.

Source画像

Target画像

Poisson 方程式を

) 解く！

, (x y h h

(2)

微分(導関数)の近似：差分法



) ( ) lim ( ) (

0

x f x f x

x

f   







微分の定義：



テーラー展開：



 ) ( ) ) (

(   



 f x f x

x x f



１階微分の前進1次差分近似、後退1次差分近似：

. , )

! ( ) ) ( )(

! ( ) 1 )(

)! ( 1 ( ... 1

...

) )(

(

! '' 2 ) 1 )(

( ' ) ( ) (

0 ) 1 (

) 1 (

2

b c a a n x

a a f b c n f a b a n f

a b a f a b a f a f b f

n n n n n n

n       

 















^





1 ) 1 (

2 ( )

)!

1 ( ... 1 ) (

! '' 2 ) 1 ( ' ) ( )

( ^ ^

 





 fⁿ xhⁿ

h n x f h x f x f h x f



 ) ( ) ) (

( f x f x

x x f



 



2次以降の項で打ち切ると…

) ) ( ( ) ) ( ( ' ) ( ) ( ' ) ( )

( ² Oh²

h x f h x x f f h O h x f x f h x

f         

誤差はhの２乗に比例:O(h²)







h x ) (a h f 

a

xx

f'()|

h ) (x f f(a)

a

１階微分は接線

微分(導関数)の近似：差分法2



高階微分の近似はより多くの評価点が必要：例えば

2階微分の前進1次差分近似：

) ) ( ( ) ( ) ( ) 2 )) ( ) ( ( ) ( ( ) ( ' ) (

'' 2 2 ²

2 Oh

h x f h x f h

h x f h x h f h O

x f h x f dx x d dxf x d

f   

 



 

 

 



１階の前進１次差分近似

) ) ( ( ) ( 2 ) 2 ) ( (

'' ₂ Oh²

h x f h x f h x x f

f       誤差はhの２乗に比例:O(h²)

x ) 2 (a h f 

a

xx

f''()|

h ) (x f f(a)

a

2階微分はラプラシアン ) (a h f 



１次精度の差分近似はｎ階の微分をn+1の評価点で近似する.



同様にn次精度の差分近似は１階の微分をn+1の評価点で近似する（次のスライド）.

微分(導関数)の近似：差分法3



高次の近似もより多くの評価点が必要：例えば１階微分の前進2次差分近似：

) ( ) (

! '' 2 1 ) ( ) ) ( ( ' ) ( ) (

! '' 2 ) 1 ( ' ) ( )

( ² ³ f xh Oh³

h x f h x x f f h O h x f h x f x f h x

f           

) ) ( ( ) ( ) ( ) 2 )) ( ) ( ( ) ( ( ) ( ' ) (

'' 2 2 ²

2 Oh

h x f h x f h

h x f h x h f h O

x f h x f dx x d dxf x d

f   

 



 

 

 



) ( )

! ( 2 ) 1 ) ( ) ( 2 ) 2 ( (

! 2 1 ) ( ) ) ( (

' Oh²h Oh³

h x f h x f h x f h

x f h x x f

f          

) 2 (

) ( 3 ) ( 4 ) 2 ) (

(

' O h

³

h

x f h x f h x x f

f       

代入

注目！ｈの２乗の誤差がhの３乗の誤差になる！

１階の前進１次差分近似

２階の前進１次差分近似

微分(導関数)の近似：差分法４

) ) ( ( ) ( 2 ) 2 ) ( (

'' ₂ Oh²

h x f h x f h x x f

f      

誤差はhの3乗に比例: O(h³)



中心差分を使うと、評価点の数は同じで、より高精度になる：例えば2階微分の中心2次差分近似：

x ) (a h f 

a

xx

f''()|

h ) (x f f(ah)

a

2階微分はラプラシアン )

(a f

) ) ( ( ) ( 2 ) ) ( (

'' ₂ Oh³

h

h x f x f h x x f

f      

) , ) (

(

) ( 2 ) ( ) ( 2 ) ( ) 2 ( ''

2 ) ( '' ) ( ) ( ) ( ) ( ) (

/ ) ( ) ( 2 '' ) 1 ( ' / ) ( / ) (

/ ) ( ) ( 2 '' ) 1 ( ' / ) ( / ) ( ) ( ) (

! '' 2 ) 1 ( ' ) ( ) (

) ( ) (

! '' 2 ) 1 ( ' ) ( ) (

3 3

3 2

b a b O

a ab

b x af x f b a a x x bf f

x f b a ab

b x af a x bf

b b O b x f x f b x f b b x f

a a O a x f x f a x f a a x f b O b x f b x f x f b x f

a O a x f a x f x f a x f

 







 



 

 





 



















 























２階の微分の前進１次差分：

←微小距離が異なるときの中心差分近似

h b a 



GradientとLaplacianの離散化(差分近似)

2 2 2 2

y I x I I



 



 

) 

,

( y

I x I I



 



) ) ( ( ) ( 2 ) ) (

( ³

1 2 Oh

h

h x f x f h x f f

n

i i

i i i i

i    









x

} ,..., , {x₁ x₂ x_n

 x

x f x f x f f

n 2 2 2 2 2 2 1 2

... 

 

 











１次精度前進：

 ２次精度中心：画像ではこれが基本.

)) , ( ) 1 , ( ), , ( ) , 1 (

( I x y I x y I x y I x y

I     



) 1 , ( 2 ) 2 , ( ) , ( 2 ) , 1 ( 2 ) , 2

(        



 I I x y I x y I x y I x y I x y 2 )

) 1 , ( ) 1 , , ( 2

) , 1 ( ) , 1

( (      



 I I x y I x y I x y I x y ) , ( 4 ) 1 , ( ) 1 , ( ) , 1 ( ) , 1

( x y I x y I x y I x y I x y

I

I         



数学的にはh<<1だが、画像などはh=1を良く使う.

GradientとLaplacianの離散化(差分近似)2



下記１階微分の中心差分近似の分母になぜ２が出てくるのかは、２階微分の中心差分と同様に計算するとわかる.

2 )

) 1 , ( ) 1 , , ( 2

) , 1 ( ) , 1

( (      



 I I x y I x y I x y I x y

h h x f h x x f f

b a b O

a ab

b x f x f a x f a x f x b f

ab x f b a b a

x f a b b a

b x f a a x f b

b b O x f b x f b x f b b x f

a a O x f a x f a x f a a x f b O b x f b x f x f b x f

a O a x f a x f x f a x f

2 ) ( ) ) ( ( '

) , ) (

(

)) ( ) ( ( )) ( ) ( ) ( ( '

) ( ' ) ( ) ( ) ( ) ( ) (

/ ) ( ) ( 2 '' / 1 ) ( ' / ) ( / ) (

/ ) ( ) ( 2 '' / 1 ) ( ' / ) ( / ) ( ) ( ) (

! '' 2 ) 1 ( ' ) ( ) (

) ( ) (

! '' 2 ) 1 ( ' ) ( ) (

3 2 3

2

2 2

2 2 2 2

2 2

2 3 2

2

2 3 2

2

3 2



 



 







 



 

 



 



















 





















h b a  



(3)

画像では3x3の作用素(オペレータ)

 微分フィルタはエッジを検出できるが，ノイズに対しても敏感

.

 ノイズを抑えながらエッジ抽出

:

 微分と平滑化の組み合わせ.

 横と縦の組み合わせ.

©CG-ARTS協会

エッジ強度画像を白をエッジとするか、

黒をエッジとするかは、表現の違い.

前進１次後退１次中心２次前進１次後退１次中心２次

微分と平滑化オペレータの合成

プリューウィットオペレータ:

- x方向微分：横に微分＋縦に平滑化. - y方向微分: 縦に微分＋横に平滑化.

©CG-ARTS協会

＊:畳み込み

重要：ソーベルオペレータ

最もよく使われている一階の偏微分オペレータ:

- 平滑化を中央に重み付

.

©CG-ARTS協会

勾配の細線化による線検出

勾配強度画像を細線化してもエッジ検出が可能.

2 2

y

x

I

I  



２階微分オペレータ(Laplacian)

２次の中心差分で

Laplace

オペレータの４連結での近似

:

2 2 2 2 2

y

x 

 



 















©CG-ARTS協会

©CG-ARTS協会中心２次

) , ( 4 ) 1 , ( ) 1 , ( ) , 1 ( ) , 1

( x y I x y I x y I x y I x y

I

I         



２階微分オペレータ(Laplacian)２

８連結では？

中心２次８連結

)) 1 , 1 ( ) 1 , 1 ( ) 1 , 1 ( ) 1 , 1 ( ( 5 . 0

) , ( 6 ) 1 , ( ) 1 , ( ) , 1 ( ) , 1 (

























y x I y x I y x I y x I

y x I y x I y x I y x I y x I I

1 1 1 1

5 . 0

) ) ( ( ) ( 2 ) ) (

(

³

1 2

O h

h

h x f x f h x f f

n

i i

i i i i

i

    



 



x

5 . 0

5 .

0 0 . 5

 6

2 1 1

5 . 1 0

2 

₂



 h

h

斜めのhは上下左右のhは1、

(4)

重要：平滑化とエッジ



平滑化(滑らかにする事)の意味は？



前回の講義でやったノイズと何の関係があるの？



前々回の周波数分解との関係は？

平滑化・Low Path→ノイズの除去・エッジの削除



変分法(Variational Calculus): 極小、極大を汎関数で停留条件を満たす様に求め、対応する偏微分方程式を導出.



汎関数(functional): 関数の関数.



停留条件:第一変分がゼロ(関数での１階微分がゼロ).



復習：極小、極大:

変分法

0 min

) , , , , ( )) , (

(     



E dxdy

u u u y x F y

x u

E

_x _y









 I dxdy

y x I

E

²

2 )) 1 , (

例(エッジ強度の積分)：

(

極大値

極小値 Euler-Lagrange方程式

) (x f

0 ) ( 'x  f

0 ) ( ' x  f

0 ) ( ''x  f 0

) ( ''x 

f f''(x)0

積分の領域：の境界曲線ではとする.

変分法２

) , ( ) , ( ) , ,

( x y u x y x y

U    

0 min

) , , , , ( )) , (

(     



E dxdy

u u u y x F y

x u

E

_x _y





 





 F x yUU U dxdy F x yu u u dxdy

y x U

E( ( , )) ( , , , x, y) ( , , , x x, y y)

摂動(微小変化)した比較関数で置き換え、その偏微分をゼロとする:





 





 



 



 



 



 



 



 

 

 

 

u dxdy F u

F u F

U dxdy U

F U U

F U U F

dxdy U U U y x F u

E E

y y x x

y x

) (

) , , , , ( )

(

0 0 0





 

 



 

←全微分

0 ) , ( x y 

 

 

↑第一変分

変分法３







 



 



 



 



 



u dxdy F y u

F dxdy x

u F y u F x

u dxdy F u F

y x y

x y y x x

) (

)) ( ) ( (

) (



 









 



 



 



 



y y

x x

u F y u F y u

F

u F x u

F x u

F



) (



準備:





 

 

 



 



 



 ) ( ) 0

( dx

u dy F u dxdy F

u F y u F

x x y x y



 ←グリーンの定理：面積

分を線積分に変換.

積分の定義域：の境界曲線

 

では

 ( x , y )  0

なので.

グリーンの定理→



©wikipedia



  



 



 

 

 



 dxdy

u F y u F dxdy x

u F u F

y x y

y x x

) (

)

(   

変分法４







 



 



 



 



 



 

u dxdy F y u

F x u F

u dxdy F u

F u E F

y x

y y x x

)) ( ) ( (

) (





↑がゼロになるためには、下記偏微分方程式を満たす:

0 ) ( )

( 



 



 



y

x

u

F y u

F x u F

 ２変数で1階微分の汎関数に対する公式：

min )

, , , , ( )) , (

(   



dxdy u u u y x F y

x u

E

_x _y

Euler-Lagrange 方程式

変分法５

0 ) ( )

( 



 



 



y

x

u

F y u

F x u F

エッジ強度の積分エネルギーの場合は:

2 min

1 2









dxdy

I

( )

2 1 2

1

2 2 2

y

x

I

I I

F    

y y x x

I I I F I F I

F 



 



 



 0 , ,

y I I x I I I y x I F y I F x I F

y x y

x





 







 

 



 

 



 



 



2 2 2 2

) ( ) (

つまり、Laplace方程式の解がエネルギーを最小化する:

0 ) ,

( 

 I x y

(5)

変分法６

0 ) ( )

( 



 



 



y

x

u

F y u

F x u F

例えば、Poisson方程式は:

min ))

, ( 2 2 (

1



 2 



dxdy y x If

I

( 2 )

2 1 I

2

I

2

If

F 

_x



_y



y y x x

I I I F I f F I

F 



 



 



 , ,

I f y f

I x f I yI xI I F y I F x I F

y x y x







 









 

 



 

 



 



 



2 2 2 2

) ( ) (

つまり、Poisson方程式の解がエネルギーを最小化する:

) , ( ) ,

( x y f x y

I 



重要：変分法によるエネルギー最小化



Laplace方程式はディリクレ・エネルギーを最小化する事で導かれ、その解は調和写像と呼ばれる.

 ディリクレ・エネルギー: 勾配の大きさを積分→エッジの大

きさの和＝凹凸具合＝ノイズの大きさ.

2 min

1  

²





dxdy I

) ,

( ) , ( ) , (

) , ( 0 ) , (

















y x y x g y x I

y x y

x I





) , ( x y g



変分法によりLaplace方程式がディリクレ・エネルギーの Euler-Lagrange方程式として導かれる↓

定義域境界:

定義域:

境界条件:

凹凸具合が最小化＝滑らかな(調和な)解＝平均化.

重要：拡散方程式



ディリクレエネルギーの最小化過程は拡散方程式(熱伝導方程式)として記述出来る:



時間の変数を加えて関数(画像)を拡張しその時間方向への接線(時間変数での一階微分)が拡散のスピードになる.

) , , ) (

, ,

( I x y t

t t y x

I  









拡散過程

t

は、時間の極限でLaplace方程式を満たす＝解は調和関数となる (

 I  0

).

例：クーラーを止めたら、温度は一定(又は周りの部屋の温度の平均)になる.

LaplacianとDoG、LoG

).

, ( ) , ( ) ,

( x y g x y g

₂

x y

DoG

_



_



_

©wikipedia

2 ) 2 exp(

) 1 ,

( 2

2 2

2 

 

y y x

x

g   

) , (

* ) , ( ) , ,

( x y g x y I x y

L  

_



 



 



 



 

 



 





  k

y x L k y x L k

y x L k y x L y x y L x

L k

) , , ( ) , , ( ) , , ( ) , , lim ( ) , , ) ( , ,

( 1

) , , ) ( , ,

( I x yt

t t y x

I 



) , (

* ) , (

* )) , ( ) , ( (

) , , ( ) 1 ( ) , (

* ) 1

( ²

2

y x I DoG y x I y x g y x g

y x L k y

x I LoG k

k 



  



















LoG: Laplacian of Gaussian.



DoG: Difference of Gaussian.

) , (xy g LoG 

差分による拡散方程式の離散化 ) , , ) (

, ,

( I x y t

t t y x

I  



) , ( )

, ( ) ,

1

( i j I i j I i j

I

ⁿ^



ⁿ

  

ⁿ

) , ( 4 ) 1 , ( ) 1 , ( ) , 1 ( ) , 1

(x y I x y Ixy Ix y I xy

I

I        



)) 1 , 1 ( ) 1 , 1 ( ) 1 , 1 ( ) 1 , 1 ( ( 5 . 0

) , ( 6 ) 1 , ( ) 1 , ( ) , 1 ( ) , 1 (

























y x I y x I y x I y x I

y x I y x I y x I y x I y x I I



時間は、前進１次差分近似：epsilonは微小時間.

4連結

８連結



空間は中心2次差分近似:

＝  

差分の陽解法と陰解法

) , ( )

, ( ) ,

1

( i j I i j I i j

I

ⁿ^



ⁿ

  

ⁿ



陽解法は…



陰解法は…

) , ( )

, ( ) ,

(

¹

1

i j I i j I i j

I

ⁿ^



ⁿ

  

ⁿ^

注目！



連立方程式となる：

b x b

x ¹

1 1

1

( , ) ( , ) ( , ) ( ) ( , )

) (



























A A

j i I I

j i I j i I j i I

I 

ⁿ ⁿ ⁿ



ⁿ

A

) , ( 4 ) 1 , ( ) 1 , ( ) , 1 ( ) , 1

(x y I x y Ixy Ix y I xy

I

I        

 なら

行列の対角成分は、非対角成分は、

０か-1の疎な行列になる.

1   4

(6)

差分の陽解法と陰解法2



画像の各画素にIDをk=(i×sx+j)と与えると…

k

1

k k1 k2

sx k

sx k

1

sx k

2

sx k

1

sx k

1

sx k

2

sx k

1

sx k

sx k2

1 2



 sx k

2 2



 sx k 1

2



 sx k )

, ( 4 ) 1 , ( ) 1 , (

) , 1 ( ) , 1 (

y x I y x I y x I

y x I y x I I













は、

) ( 4 ) ( ) 1 (

) ( ) 1 (

k I sx k I k I

sx k I k I









 となる、ここで、 × の行列を考える.

sx sy

差分の陽解法と陰解法3

) ( 4 ) ( ) 1 (

) ( ) 1 (

k I sx k I k I

sx k I k I I





























4 4 0 1 0 1 4 1 1 0 0

1 4 1 0 0

1 4 1 0

0 1 4 1

1 0 1 0 1 0 0 0 1 4













1行にはゼロでない要素が５つで対角成分が-4、それ以外は1: ラプラス行列とも呼ばれる.

sxsy

sxsy k

sx k 

 1 sx k k 

 1 k

)

(  

 I  A

A

行列の対角成分

は、

非対角成分は、０か-1の疎な行列になる.

4 1  

 

フィルタの繰り返し適用

)

1 ( n

n Filtering I I ^ 

}

;

);

( ){

iteration (

tmp I

I Filtering tmp

for



 陰解法では絶対安定だが、陽解法は微小時間epsilonが大きい値は解が不安定になる(クーラン条件). 演習では簡単の為陽解法を使う.

復習：エッジ強調フィルタ(空間領域)



周波数領域でのエッジ強調フィルタと同様に空間領域でも、元画像＋k(元画像-平滑化画像)でエッジ強調画像を作成可能.

)) , ( 1 ( 1 ) , ( 1

) ,

( u v kH u v k H u v

H

_h__emph

 

_high

  

_low

元画像平滑化画像 (Gaussian)

= ＋

エッジ強調画像

-

元画像

k( )

エッジ画像 (高周波のバンド画像) 絶対値

平滑化に用いた＋反転 sigmaスケールでのエッジ強度画像

[a] から [c] を引くと

Laplacianオペレータによるエッジ強調

©CG-ARTS協会



Laplacian

をたすと平滑化なので１から引いたオペレータは鮮鋭化(エッジ強調).

4連結 8連結

この8連結は差分法的には精度が悪い近似.

Laplacianオペレータによるエッジ強調2

©CG-ARTS協会

多重解像度解析などの周波数フィルタと原理は同じ.

(7)

Laplacianオペレータによるエッジ強調3



Multiresolutional Meshとして３D形状への応用もある.

数理モデリング・解析の基本

問題・要求エネルギー

(

積分方程式

)

エネルギーを考える(設計). 変分法.

偏微分方程式

数値解析線形化

連立方程式差分法等.

結果数値解法.

ギザギザを滑らかにしたい：ノイズを細かいエッジと考えると…

2 min

1 2









dxdy I

エッジの大きさの和を最小化.

) , , ) ( , ,

( I x y t

t t y x

I 



数値解法で実際に計算する.

) , ( ) , ( ) ,

1(i j I i j I i j

Iⁿ^  ⁿ 



 ⁿ

) , ( 4 ) 1 , ( ) 1 , ( ) , 1 ( ) , 1

(x y Ix y Ixy Ixy Ixy

I

I        



演習: 平滑化・エッジ

www.riken.jp/brict/Yoshizawa/Lectures/index.html www.riken.jp/brict/Yoshizawa/Lectures/Lec19.pdf

www.riken.jp/brict/Yoshizawa/Lectures/Ex11.zip

↑はReport05の内容です。

今日〆切のReport04が出来ていない人は質問してください！

演習19-1:エッジ強度画像の作成.

演習19-2: 拡散方程式による平滑化フィルタの作成.

演習１9－１

エッジ強度画像の作成: Ex11.zip内のEdgeMagnitude.cxx の中にあるコメントに従ってソーベル作用素を使った勾配強度画像を生成するプログラムを作成しましょう.

演習１9－2

拡散方程式による平滑化フィルタの作成: Ex11.zip内の Diffusion.cxxの中にあるコメントに従って拡散方程式による平滑化フィルタを作成しましょう.

微小時間epsilon=0.25で繰り返し10,20,30,40,50,100で実行してみましょう.

繰り返し5回繰り返し10回繰り返し50回繰り返し100回

繰り返し5回 10回 20回 30回 40回 50回 100回

来週の予定

特徴保存フィルタ: Nonlinear Diffusion, Bilateral &

Non-Local Meansフィルタ.