f n I g

情報デザイン専攻
Shin Yoshizawa: [email protected]
今日の授業内容
画像情報処理論及び演習II
www.riken.jp/brict/Yoshizawa/Lectures/index.html
www.riken.jp/brict/Yoshizawa/Lectures/Lec17.pdf
-フィルタ処理・エッジ強調線形フィルタとノイズ
吉澤 信
第5回講義
水曜日１限
教室6218
[email protected], 非常勤講師
1.
線形フィルタの基礎.
2.
画像復元の基礎.
3.
演習：線形フィルタ.
大妻女子大学 社会情報学部
今日の演習は第３回(Report0６)のレポートで出すので、
みなさん頑張ってくださいねーp(^^)q
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
ノイズ除去・画像復元・フィルタリングの背景


自然科学では観察・観測による画像解析が重要である.
測定・計測に基づくCAD/CAM/CAEも注目されている.
重要：線形フィルタ
 線形フィルタ（畳み込み和、Convolution）:
hy
hx
  f (i  m, j  n ) h ( m, n )
I (i , j ) 
m   h y n   hx
f 入力画像
I 出力画像
MRI
CT
j
構造光式表面
形状スキャナー
(0,0)
レーザー式表面
形状スキャナー
(i  m, j  n )
i
(i , j )
入力画像
h カーネル画像：フィルタ
hx
hx
hy
( m, n )
カーネル
画像
n
(0,0)
hy
m
カーネル画像(局所Window)サイズ：
共焦点レーザー顕微鏡
高速レンジスキャナー
(i, j ) フィルタを適用している画素の座標値
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
Normalized Convolution
重要：線形フィルタ
 線形フィルタ（畳み込み和、Convolution）:
I (i , j ) 
hy
hx
  f ( i  m, j  n ) h ( m, n )
m   h y n   hx
h
f
(2h y  1)  (2hx  1)
I
 fの重み付和
＝fとhの対応画
素値の積和
＝線形フィルタ.
 正規化された畳み込み和:
1
I (i , j ) 
W
W
hy
hy
hx
I new (x)  
Z (x, y) I (y)dy
 Z (x, y)dy
  f (i  m, j  n)h(m, n),
m   h y n   hx
hx
  h ( m, n )
m   h y n   hx
重みの和が１のときPartition of Unityを満たすと言
い、画像の輝度値が定数だった場合にフィルタ後の
輝度値も同じ定数である事を保証する(constant
reproductively)＝平均である. a *100  b *100  c *100
a *100  b *100  c *100
abc
 100
1
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
畳み込み：空間 VS 周波数

モザイクフィルタ
線形フィルタ（畳みこみ和）:
g (i, j ) 
w
ある範囲を
一色へ置き
換える.
w
  f (i  m, j  n)h(m, n)
n w m w
畳み込み積分の離散版.

畳み込み積分(convolution):
 
f1 ( x, y ) * f 2 ( x, y ) 
  f ( , ) f
1
2
( x   , y   )dd
©CG-ARTS協会
©CG-ARTS協会
  
畳み込みは周波数領
域では掛け算になる！
f * g  F 1[ F [ f ]F [ g ]]
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
Median Pyramid
メディアン(中央値)フィルタ

Block Median Pyramidal Transform: Cartoonish：
- ResizeとUpsamplingは３次補間などを使う.
- DownsamplingはMedianフィルタを使う.
- エッジをある程度保存.
V. Melink et al.
IEEE ICIP’99.
©CG-ARTS協会
 ある範囲の中央値に置換.
 異常値を検出.
 エッジを(ある程度)保存.
Resize+Median
+Upsampling
©CG-ARTS協会
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
平均化フィルタ
畳み込み：空間 VS 周波数２
ある範囲の
平均.
have ( x, y ) 
1
 x y
rect  , 
w2
 w w
1
1 x  2 かつ
rect  x, y   
0 その他
1
y  のとき
2
sin wu sin wv
wu
wv
平均化フィルタ
周波数領域
空間領域
©CG-ARTS協会
H ave (u , v) 
©CG-ARTS協会
2
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
Anisotropic(異方性)フィルタ
Gaussianフィルタ
 カーネル画像をデザインする事により特
定方向の重み付平均化を行える.
©CG-ARTS協会
Gauss関数の
重み付平均.
g ( x , y ) 
©CG-ARTS協会
Shin Yoshizawa: [email protected]
 x2  y2 

exp 
2
2 2 

1
2
Shin Yoshizawa: [email protected]
Gaussianフィルタの性質
 フーリエ変換：
F [ g ( x )]  F [
Gaussianフィルタの性質２
1
x2
1
2
exp( 2 )]( ) 
exp(  2 )
2
2
2 
2
 繰り返し適用: Gaussianフィルタの繰り返し適用は異なる
sigmaのGaussianフィルタ.
Gaussianのフーリエ変換はGaussian：ただしsigmaは反比例.
周波数
領域
g ( x ) * ( g ( x ) * f ( x )) 
1
g
2
2
( x) * f ( x)
g ( x ) * ( g  ( x ) * f ( x ))  g  ( x ) * ( g ( x ) * f ( x ))  ( g  ( x ) * g ( x )) * f ( x )
暈け大
暈け小
大きな
sigma
小さな
sigma
暈け大
2
2 2
1
 )] * f ( x )
exp(  2 ) exp(
2
2
2
2

1
 F 1[ exp(
( 2   2 ))] * f ( x )
2
2
x2
1

exp(
) * f ( x)
2( 2   2 )
2 (  2   2 )

暈け小
Shin Yoshizawa: [email protected]
1
g
2
 2  2
( x) * f ( x)
小さいsigmaでの繰り返し適用
大きいsigmaで１回適用
Shin Yoshizawa: [email protected]
Gaussianフィルタの性質３
x2  y 2
 分離:
g ( x, y ) 
exp(
)
2 2
2 2
2
高次元のGaussianは
1
x
1
y2
低次元Gaussianの積.  ( 2  exp( 2 2 ))( 2  exp( 2 2 ))
 g ( x ) g ( y )
1
 異方性フィルタ:
…
 F 1[
＝
空間
領域
 F 1[ F [ g  ( x )]F [ g ( x )]] * f ( x )
Gaussianスケールスペース



g ( x, y ) 
x 2  y 2 | (0,0)  ( x, y ) |2 | x |2 , x  ( x, y )T
a  x 
1
xT Ax
a
exp(
)
x Ax  ( x, y ) 11 12    g (x) 
2
2
2 2

a
a
y
 21 22  
T
a  x 
a
 1 0

xT Ax  ( x, y ) 11 12    A  
a
a
y
 0 1
 21 22  
 xT Ax  x 2  y 2
異なるσでスムージングすると、様々なレベルでの画像.
σを“スケール”として、「異なるスケールの画像」.
Gaussianフィルタを繰り返し適用:
1
2
exp(
2
x2  y2
)
2 2
世界は、スケールによって
異なる構造をもつ、という考え方.
 σで偏微分する事で
スケールの極値や
変化率を図る.

3
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
エッジ強調フィルタ(空間領域)
画像復元
 周波数領域でのエッジ強調フィルタと同様に空間領域で
も、元画像＋k(元画像-平滑化画像)でエッジ強調画像を
作成可能.
H h emph (u, v )  1  kH high (u, v )  1  k (1  H low (u, v ))
-
＋ k(
=
元画像
元画像
エッジ強調画像
平滑化に用いた
sigmaスケールでの
エッジ強度画像
)
劣化(暈けやノイズ)
画像の復元
＝ノイズ除去.
平滑化画像
(Gaussian)
エッジ画像
(高周波の
バンド画像)
絶対値
＋反転
©CG-ARTS協会
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
Noise ?

Noise ?
例えば加算ノイズモデル(Additive Noise Model):
True
Signal:
f
元信号
I  g* f n
Degradation
Function: g
e.g. PSF
Restoraed
復元信号 Signal: Iˆ
©wikipedia
Observed
Signal: I
得られる信号
+
観察・観測・測定
装置による暈け Noise:
Photon (Shot) Noise
Salt and Pepper Noise
(Impulse)
n
復元
フィルタ
Shin Yoshizawa: [email protected]
Additive Gaussian Noise
Brownian Noise
Periodic Noise
Multiplicative Speckle Noise
©S. Dangeti, 2003.
Shin Yoshizawa: [email protected]
Impulse Noise (Salt and Pepper)
©G. Baker, 2005.
Adaptive
Median Filtering
R. C. Gonzalez
R. E. Woods
Digital Image Processing
pp. 332-333, 2008
Adaptive Median Filter
 アルゴリズム:各画素で、外れ値を判別＋半径を増加.
r : 局所Windowの半径
A) Stage A:
rmax : 半径の最大値
1. : A1  I med  I min
I ( x, y ) : 輝度値
A

I

I
2. : 2
med
max
3. : if A1  0  A2  0 goto Stage B.
4. : r  r  1 , if r  rmax goto A-1, else output I med
B) Stage B:
I
の出力が外れ値だった場合.
1. : B1  I ( x, y )  I min
2. : B2  I ( x, y )  I max
3. if B1  0  B2  0 output I ( x, y ) else output I med
med
局所Window内の
I med : 中央値
I min : 最小値
I max : 最大値
I min
I med
I max
4
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
Deconvolution(逆畳み込み)
加算性白色ガウスノイズ
I  f n
*
白色性：パワースペクトル(フーリエ変換の自己相関)が一定
+
Input
f  F 1[
* Degradation: Gaussian
PSF: point spread function
I  f *g
Deconvoluved
Output
=
ガウス分布
f
n
I
F[I ]
]
F[g]
©G. Baker, 2005.
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
Additive Noiseの確率密度関数
 輝度値のヒストグラムの可視化で、
加算ノイズの分布(確率密度関
数)がある程度わかる．
Wiener Filter
I  f n

ノイズをｎ、観測信号をI、真の信号をfとして加算ノイズ：
I (t )  f ( t )  n ( t )
を考える. 最小二乗的に最適なフィルタをhとして、復元信号r
とｆの差の二乗を最小化する. r (t )  h (t ) * ( f (t )  n (t ))
2

E(t)   r (t )  f (t ) dt  min
©R. C. Gonzalez and R. E. Woods, 2008.


ｆとnは無相関なので





 f (t )n(t )dt  0





2
2 2
2
2 2
2
2
2 2
2 2
2
2
2
2
2 2
 | h( f  n)  f | dt   (h f  2h fn  h n  2 f h  2 fhn  f )dt   (h f  h n  2 f h  f )dt   ( f (1  h)  h n )dt

f2
E (t )  0  2 f 2 (1  h )  2hn 2  0   f 2  ( f 2  n 2 )h  0  h  2
h
f  n2
Gaussian
Rayleigh
Gamma
Exponential
Uniform
Impulse
Shin Yoshizawa: [email protected]
1
F [ f ]2

F [ f ]2  F [n ]2 1  F [n ]2 / F [ f ]2
通常は，fとnは不明なので，
適当な定数で近似.
Shin Yoshizawa: [email protected]
Wiener Filter2

F [h(t )] 
Wiener Filter2
ノイズをn、観測信号をI、真の信号をf、PSFをgとすると、
©J. Robles, cssvc.ecsu.edu
©CG-ARTS協会
I ( u, v )  g ( u, v ) * f (u, v )  n (u, v )

復元画像と原画像の誤差を最小にするような逆フィルター
ｈ(u,v)はPSFが無い場合と同様にして、
F [h(t )] 
F [ f ]2 F [ g ]
1
F [ g ]2

2
2
2
2
F [ f ] F [ g ]  F [n ]
F [ g ] F [ g ]  F [n ]2 / F [ f ]2
F [h(t )] 
1
F [ g ]2
F [ g ] F [ g ]2  
©wikipedia
通常は，fとnは不明なので，
適当な定数で近似: 単純な
逆畳み込みと比べて、F[g]が
ゼロに近くても発散しない.
5
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
演習: 線形フィルタ
www.riken.jp/brict/Yoshizawa/Lectures/index.html
www.riken.jp/brict/Yoshizawa/Lectures/Lec17.pdf
www.riken.jp/brict/Yoshizawa/Lectures/Ex10.zip
演習１7－１


平均化フィルタの作成: Ex10.zip内のAverageFilter.cxxの
中にあるコメントに従って平均化フィルタを作成しましょう.
カーネル画像(局所Window)のサイズをr=１，５，１０で実行
してみましょう.
カーネル画像(局所Window)のサイズ：(2r+1)x(2r+1)
演習17-1:平均化フィルタの作成.
演習17-2: 正規化Gaussianフィルタの作成.
演習17-3: エッジ強調フィルタの作成.
11x11: r=5
21x21: r=10
61x61: r=30
今日からはReport0６の内容です。
Report04・０５が出来ていない人も、演習
17-1だけは先にやりましょう！
Shin Yoshizawa: [email protected]
Shin Yoshizawa: [email protected]
演習１7－2


正規化Gaussianフィルタの作成: Ex10.zip内の
GaussianFilter.cxxの中にあるコメントに従って正規化
Gaussianフィルタを作成しましょう.
パラメータをr=5，sigma=2.5、及びr=10,sigma=5.0で実行し
てみましょう.
ガウス関数の標準偏差パラメータ: sigma
21x21: r=10, sigma=5.0
41x41: r=20, sigma=10.0 121x121: r=60, sigma=30.0
演習１7－3


エッジ強調フィルタ(＝元画像+k(元画像-平滑化画像))の
作成: Ex10.zip内のEdgeEnhancementFilter.cxxの中にある
コメントに従ってエッジ強調フィルタを作成しましょう.
パラメータr=10, sigma=5.0でk=1.5,3.0,5.0の三種類につい
エッジ強調の強度: k
て実行してみてください.
k=1.5
21x21: r=10, sigma=5.0
k=3.0
k=5.0
6