画像処理・実習第一回：

画像処理・実習
第十一回：
画像の直交変換（二次元フーリエ変換）
東海大学
情報理工学部情報メディア学科
濱本和彦
今回の内容
課題の解説
7. 画像の直交変換



7.1 直交変換
7.2 空間周波数とスペクトル表現
7.3 二次元フーリエ変換


離散フーリエ変換（DFT）
高速フーリエ変換（FFT）
実習


画像のフーリエ変換
フーリエ変換を利用したフィルタ処理
課題の回答例
輪郭抽出の基本：ノイズ除去→微分→２値化
Median()
Prewitt()
Fixed_threshold()
th : 128
輪郭線が途切れており，良い結果とは言えない
課題の回答例
Linear_Transform()
0-70 → 0-255
Median()
Prewitt(1.7)
Fixed_threshold()
th : 210
課題の回答例
ガウシアンで平滑化（一回）
Fixed_threshold()で２値化（しきい値42）
Laplacian()で微分処理
Thinningで細線化
画像の直交変換
変換とは？？
あるデータを他のデータに変換する作業は，一般的に行列で表現されます。
今，データ［ｆＮ］を，［ＣＮ］に変換する事を考えます。Ｎ＝４とします。
変換は次のように表現できます。
C0  0, 0 0,1
 C  

1
,
0
1,1
1
 
C2  2, 0 2,1
  
C3  3, 0 3,1
0 , 2
1, 2
2 , 2
3, 2
0 , 3   f 0 

1,3   f1 
 
2 , 3   f 2 
 
3,3   f 3 
これは，例えば，4画素の画像を何か処理して新しい画像を得ることを表します。
変換とは？？
この変換の式は，次のように表現できます。
CN    NN  f N 
それでは，［ＣＮ］から［ｆＮ］を求めるにはどうすればよいでしょう？
 f N   NN  CN 
1
ですね。
直交変換？？
 f N   NN  CN 
1
この逆変換の式が，つぎのように表現できるとき，
 f N   NN  CN 
T
つまり，
NN 
1
  NN 
T
が成り立つとき，この［ΦNN］による変換を，「直交変換」といいます。
直交？？
直交変換行列には，次の特徴があります。
0 (m  n)
 m , n  m,i  n,i 
i 0
1 (m  n)
N 1
このような特徴を持つベクトルφkを，正規直交基または基底と呼びます。
直交基底を用いた変換では，変換と逆変換が同じ行列で表現出来ます。
直交変換の代表的な例がフーリエ変換です。
演習7.1
1 j 12  0 1  j 12  0
1 j 12  3 1  j 12  3
 1[i ], [i ]  e
 e
  e
 e
2
2
2
2
1 1 1
1
1
   e0  e0  e0
2 2 4
4
4
1 1 1 1
    1
4 4 4 4
*
1
1
1
1
1
1 j 2  10 1  j 2  00
1 j 2  13 1  j 2  03
*
 1[i ],0 [i ]  e
 e
  e
 e
2
2
2
2
1 1 1 j 12  11 1 1 j 12  12 1 1 j 12  13 1
   e
  e
  e

2 2 2
2 2
2 2
2
1 1
1
1
    j     1    j   0
4 4
4
4
離散フーリエ変換
離散フーリエ変換では，基底は次のように表現されます。
1
 2 
k ,i 
exp j
ki 
N
 N 
変換は次のように表現されます。＊は複素共役を意味します。
N 1
Ck   f i   *k ,i
i 0
それでは，そもそもフーリエ変換とはいったい何なんでしょう？
フーリエ変換って何？
フーリエ変換とは，時間（空間）の世界から，周波
数の世界へ変換する事です。
変換には，いろんな周波数の正弦波を用います。
だから変換基底が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
ディジタルデータに対するフーリエ変換。
コンピュータで計算→データは有限長。
得られた有限区間のデータを一周期と仮
定し，それが周期的に連続しているものと
してフーリエ変換を考える。
2 

f i exp  j
ki  (k  0,1,2,, N  1)

N 

i 0
N 1
 2 
IDFT : f i   Ck exp j
ki  (i  0,1,2,, N  1)
 N 
k 0
1
DFT : Ck 
N
N 1
具体例
２周期分含まれているので
sin成分のみなので
実部の成分は0です
離散フーリエ変換
実部
X （連続の場合）
i （離散の場合）
周波数
２サンプル目の位置に
成分が得られます。
データ区間長Ｔ
 2
sin k
 T

 2
x   sin k

 N

i

１／Ｔ : 空間周波数，k : Ｔに含まれる周期数
Ｎ : サンプリング数
１
２
１／Ｔ
虚部
周波数
DFTの対称性（周期性）
前の例と同じく，虚部の２サンプル目に成分がある場合
N
座標系（原点の位置）
に注意！
原点が中央の時
原点が左端の時
-2
N-2
2
同じデータになる
周波数
DFTの例
（0～N-1でのDFTです。原点は左端）
200
2
100
60000
1
0
0
40000
0
100
200
実部。直流分のみ。
0
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
~
F ( k , n) 
M
1
F (k , l ) 
N
M 1

m 0
N 1
f (m, n)W1
~
ln
F
(
k
,
n
)
W

2
n 0
km
二次元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_processing11.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 W2
 km
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
  F (k , l )  j
2
m
 M
l 0 k 0
2
同様に
N 1 M 1
2

 km
ln
k  W1 W2

 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
n 2
ラプラシアンフィルタ関数
 2
A(k , l )   j
 M
  2
k  j
  N
2
  Ax (k )   Ay (l ) 
2
2

l

2
ラプラシアンフィルタ実習
ラプラシアンフィルタ関数
void f_laplacian(void) を完成させなさい。
lena.rawに対する処理結果を求めなさい。
ラプラシアン結果

Download Report