Transcript 画像処理・実習 第一回 - 東海大学 情報理工学部情報メディア学科

画像処理・実習
第十二回：
画像の直交変換（二次元フーリエ変換）
東海大学
情報理工学部情報メディア学科
濱本和彦
今回の内容
7. 画像の直交変換



7.1 直交変換
7.2 空間周波数とスペクトル表現
7.3 二次元フーリエ変換


離散フーリエ変換（DFT）
高速フーリエ変換（FFT）
実習


画像のフーリエ変換
フーリエ変換を利用したフィルタ処理
フーリエ変換って何？
フーリエ変換とは，時間（空間）の世界から，周波
数の世界へ変換する事です。
変換には，いろんな周波数の正弦波を用います。
だから変換基底がexp(jΘ)の形になるのです。
オイラーの公式：ｅ ｊΘ＝cosΘ＋ｊ sinΘ
元のデータが，いろんな周波数のcos波やsin波
にどれくらい似ているかを表します。
実部がcos波，虚部がsin波による変換結果とな
ります。
具体例１
波の成分が全くないので，
フーリエ変換
直流成分
（±が変化しない）
x
周波数は０になります。
実部として扱われます。
周波数
具体例2
sin波成分は虚部に
現れるため実部は０
です。sinとcosは直交
関係にあるからです。
周波数ｆ０のsin波なので，
実部
フーリエ変換
・・
・・
周波数
周波数ｆ０の位置に成分
が
得られます。
x
sin ω０X
ｆ
０
虚部
周波数
具体例3
周波数ｆ０のcos波なので，
フーリエ変換
・・
・・
x
周波数ｆ０の位置に成分
が
得られます。
ｆ
０
実部
周波数
cos波成分は実部に
現れるため虚部は０
です。
cos ω０X
虚部
周波数
具体例4
複雑な波形でも，それが
周期ｆ０の周期波形ならば，
フーリエ変換により
各周波数成分に分解
できます。
＋
ｆ
・・
０
・・
x
ｆ０の整数倍の周波数の
sin波とcos波の和で表現
できます。
＋
フーリエ変換
ｆ
０
2ｆ０ 3ｆ０ 4ｆ０
実部
2ｆ０ 3ｆ０ 4ｆ０
虚部
周波数
周波数
例えば，こんなことが出来ます
1
0
-1
信号成分
0
100
200
元の信号
（低周波）
1
フーリエ変換
0
ノイズ成分
1
-1
0
0
-1
5
10
受信された信号
0
5
周波数
10
ノイズ
（高周波）
1
逆フーリエ変換
0
-1
0
100
200
元の信号を復元できる
ノイズ成分
カット
フーリエ変換の対称性
フーリエ変換結果F(ω)には，次のような対称性があります。
F    F ( )
*
-ｆ０
ｆ
-ｆ０
ｆ
周波数
０
実部
実部はcos成分=偶関数
周波数
０
虚部
虚部はsin成分=奇関数
離散フーリエ変換
DFT : Discrete Fourier Transform
ディジタルデータに対するフーリエ変換。
コンピュータで計算→データは有限長。
得られた有限区間のデータを一周期と仮
定し，それが周期的に連続しているものと
してフーリエ変換を考える。
1 N 1
 2 
DFT : Ck   fi exp  j ki  (k  0,1,2,, N 1)
N i 0
N 

N 1
 2 
IDFT : fi   Ck exp j ki  (i  0,1,2,, N 1)
 N 
k 0
具体例
２周期分含まれているので
sin成分のみなので
実部の成分は0です
離散フーリエ変換
実部
X （連続の場合）
i （離散の場合）
周波数
２サンプル目の位置に
成分が得られます。
データ区間長 Ｔ
 2
sin k
 T

 2 
x   sin k i 

 N 
１／Ｔ : 空間周波数，k : Ｔに含まれる周期数
Ｎ : サンプリング数
１
２
１／Ｔ
虚部
周波数
DFTの対称性（周期性）
前の例と同じく，虚部の２サンプル目に成分がある場合
N
座標系（原点の位置）
に注意！
原点が中央の時
原点が左端の時
-2
N-2
2
同じデータになる
周波数
DFTの例
（0～N-1でのDFTです。原点は左端）
200
2
100
60000
1
0
0
0
40000
0
100
200
実部。直流分のみ。
100
20000
200
0
サンプリング数N=256,
周期の数k=16
sin波+直流分
100
0
100
200
パワースペクトル
0
実部2+虚部2
-100
0
100
200
虚部。16とN-16の位置に成分。
対称性に注意。
DFTの例
（0～N-1でのDFTです。原点は左端）
2
200
100
60000
1
0
0
0
100
200
サンプリング数N=256,
周期の数k=16
cos波+直流分
40000
0
100
200
実部。直流と，16および
N-16の位置に成分。
20000
0
0
100
200
パワースペクトル
100
0
-100
0
100
200
虚部。成分は無し。
二次元離散フーリエ変換
二次元信号（画像）f(m,n)，画素数M×Nの二次元DFTは
次の式で定義されます。
1 N 1 M 1
km
ln
F (k , l ) 
f
(
m
,
n
)
W
W

1
2
M  N n0 m0
ただし
 2 
 2 
W1  exp  j , W2  exp  j 
M
N


これは，一次元DFTを２回繰り返すことであり，次のように
計算されます。
1 M 1
~
F (k , n)   f (m, n)W1km
M m 0
1 N 1 ~
F (k , l )   F (k , n)W2ln
N n 0
二次元DFTの対称性（周期性）
座標系（原点の位置）
に注意！
原点が中央の時
原点が左上の時
周波数
N
M
二次元DFTの例
-N/2～N/2のDFT。原点は中央
周期の数:2
周期の数:25
周期の数:50
二次元DFTの例
-N/2～N/2のDFT。原点は中央
周期の数:2
周期の数:25
周期の数:50
二次元DFTの例
-N/2～N/2のDFT。原点は中央
+
二次元DFTの例
-N/2～N/2のDFT。原点は中央
高速フーリエ変換
FFT : Fast Fourier Transform
重複する計算を削減することにより高速化
N個のDFT→N/2個のDFTを2回でＯＫ
N/2個のDFT→N/4個のDFTを2回
2個のDFTの計算まで簡略化できる
データ数は2n個でなければならない
データがN個の場合


DFT : N2回の乗算が必要
FFT : Nlog2N回
高速フーリエ変換
FFT : Fast Fourier Transform
x(0)
x(4)
x(2)
x(6)
x(1)
x(5)
x(3)
x(7)
2点DFT
2点DFT
2点DFT
2点DFT
G(0)
G(1)
G(2)
G(3)
H(0)
H(1)
H(2)
H(3)
W 80
W
7
8
X (0)  G(0) +W80 H(0)
X (7)  G(3) + W87 H (3)
実習
image_processing12.cとimgfft.hを利用します。
void fft2(int flag)





二次元FFTの関数です。
原点（直流成分）は，中心[128][128]です。
flag : 0 フーリエ変換
flag : 1 逆フーリエ変換
フーリエ変換時は，パワースペクトルを出力します。
lena.rawとlenan20.raw（ノイズ有り）のパワース
ペクトルを比較しましょう。
逆変換することで元の画像が得られることを確
認しましょう。
Lenaのパワースペクトル
低周波部
にエネル
ギーが集中
高周波部にも均等にエネルギーが分布
低域通過フィルタ
高域を強制的に0としてカットします。
？
void lpf(void)を実行し結果を確認しなさい。
lpf()は，必ずfft2(0)の後に実行すること。
低域通過フィルタ結果
縦，横
とも20で
カット
縦，横
とも50で
カット
ほぼ移動平均平滑化に等しい
高域通過フィルタ
低域を強制的に0としてカットします。ただし，直
流分は残します。
？
void hpf(void)を実行し結果を確認しなさい。
hpf()は，必ずfft2(0)の後に実行すること。
高域通過フィルタ
縦方向，横方向とも均等に低域カットした場合の
画像を求めなさい。
縦方向の高域のみ（横方向は全てカット），横方
向の高域のみ（縦方向は全てカット），それぞれ
の画像を求めなさい。
直流分が無くなったらどのような画像になります
か？
ノイズがある場合にこの処理を施したらどのよう
な結果になるでしょうか？
高域通過フィルタ結果
均等にカットした場合
縦，横ともに5でカット
縦，横ともに10でカット
より低域が少なくなり，
エッジが目立つ画像と
なる。
高域通過フィルタ結果
縦方向高域のみ（横方向全カット） 横方向高域のみ（縦方向全カット）
縦方向は5でカット
横方向は5でカット
横線のエッジのみ抽出
縦線のエッジはほぼ消滅
縦線のエッジのみ抽出
横線のエッジはほぼ消滅
高域通過フィルタ結果
縦，横ともに10でカット
直流分あり
縦，横ともに10でカット
直流分なし
エッジだけの画像となる
高域通過フィルタ結果
ノイズがある場合
縦，横とも
10でカット
ごま塩状の雑音はほぼそのまま
残ってしまっている
周波数空間におけるX,Y方向微分
N 1 M 1
f (m, n)   F (k , l )W1kmW2 ln
l 0 k 0
 2   2 
  F (k , l ) exp j
km exp j ln 
 M
  N 
l 0 k 0
N 1 M 1
f (m, n) N 1 M 1
2
km
ln
  F (k , l ) j
k W1 W2
M
l 0 k 0
x方向（離散m方向）微分 m
これを逆変換すればX方向微分になる
f (m, n) N 1 M 1
2
km
ln
  F (k , l ) j l W1 W2
n
N
l 0 k 0
y方向（離散n方向）微分
同様に
周波数空間におけるX,Y方向微分
X方向フィルタ関数（Y方向も同様）
Ax(k)
-M/2
M/2
k
2
Ax (k )  j k
M
微分処理 実習
x方向微分の関数 void x_diffrentiation(void)
y方向微分の関数 void y_diffrentiation(void)
を完成させなさい。フィルタ関数が複素数であ
る事に注意しなさい。複素数と複素数の乗算に
なります。
lena.rawについて，それぞれの微分処理を
行った画像を出力しなさい。
微分処理結果
X方向微分
縦方向のエッジ抽出
Y方向微分
横方向のエッジ抽出
周波数空間におけるラプラシアン
f (m, n) N 1 M 1
2
km
ln
  F (k , l ) j
k W1 W2
m
M
l 0 k 0
 f (m, n)
 2  km ln
  F (k, l )  j k  W1 W2
2
m
 M 
l 0 k 0
2
同様に
N 1 M 1
2
 f (m, n)
 2  km ln
  F (k, l )  j l  W1 W2
2
n
 N 
l 0 k 0
2
N 1 M 1
∴ ラプラシアン＝
2
2 f (m, n) 2 f (m, n)
+
2
m
n2
ラプラシアンフィルタ関数
 2   2 
A(k , l )   j
k + j l
 M   N 
2
2
  Ax (k ) + Ay (l )
2
2
ラプラシアンフィルタ 実習
ラプラシアンフィルタ関数
void f_laplacian(void) を完成させなさい。
lena.rawに対する処理結果を求めなさい。
ラプラシアン結果